Barbara Porro

Thamilla Zamoum

Cédric Mallien

Benjamin C C Hume

Christian R Voolstra

Eric Röttinger

Paola Furla

Didier Forcioli
email: forcioli@unice.fr

Horizontal acquisition of Symbiodiniaceae in the Anemonia viridis (Cnidaria, Anthozoa) species complex

Keywords: Symbiosis, Clonality, Cryptic species, RAD-Seq, microsatellite, ITS2 ORCID

All metazoans are in fact holobionts, resulting from the association of several organisms, and organismal adaptation is then due to the composite response of this association to the environment. Deciphering the mechanisms of symbiont acquisition in a holobiont is therefore essential to understanding the extent of its adaptive capacities. In cnidarians, some species acquire their photosynthetic symbionts directly from their parents (vertical transmission) but may also acquire symbionts from the environment (horizontal acquisition) at the adult stage. The Mediterranean snakelocks sea anemone, Anemonia viridis (Forskål, 1775), passes down symbionts from one generation to the next by vertical transmission, but the capacity for such horizontal acquisition is still unexplored. To unravel the flexibility of the association between the different host lineages identified in A. viridis and its Symbiodiniaceae, we genotyped both the animal hosts and their symbiont communities in members of host clones in five different locations in the North Western Mediterranean Sea.

The composition of within-host symbiont populations was more dependent on the geographical origin of the hosts than their membership to a given lineage or even to a given clone. Additionally, similarities in host symbiont communities were greater among genets (i.e. among different clones) than among ramets (i.e. among members of the same given clonal genotype). Taken together, our results demonstrate that A. viridis may form associations with a range of symbiotic dinoflagellates and suggest a capacity for horizontal acquisition. A mixed-mode transmission strategy in A. viridis, as we posit here, may help explain the large phenotypic plasticity that characterises this anemone.

INTRODUCTION

From their earliest development stages, all metazoans live in close association with different microorganisms (symbionts) forming a so-called holobiont (Margulis & Fester, 1991). The evolution of these holobionts is influenced by the transmission mode of symbionts. Two modes of symbiont acquisition have been described: from the parents or from the environment, called vertical transmission and horizontal acquisition, respectively. Vertical transmission is traditionally more associated with mutualism whereas horizontal acquisition corresponds more to parasitism [START_REF] Moran | Lifestyle evolution in symbiotic bacteria: Insights from genomics[END_REF]. Vertical transmission indeed ensures the maintenance of beneficial associations at the embryonic stage [START_REF] Zilber-Rosenberg | Role of microorganisms in the evolution of animals and plants: The hologenome theory of evolution[END_REF], but with less flexibility in the host-symbiont association. A strong correlation between the host and symbionts phylogenies can therefore be expected (for a review on the topic, see [START_REF] Baumann | Biology of Bacteriocyte-Associated Endosymbionts of Plant Sap-Sucking Insects[END_REF]. However, this binary vision of a clearly split transmission mode between mutualism and parasitism is not so obvious and many examples of mutualistic symbiosis maintained by horizontal acquisition have been observed and described [START_REF] Henry | Horizontally Transmitted Symbionts and Host Colonization of Ecological Niches[END_REF][START_REF] Kaufman | Bacterial Symbionts Colonize the Accessory Nidamental Gland of the Squid Loligo opalescens via Horizontal Transmission[END_REF][START_REF] Nelsen | Dissociation and horizontal transmission of codispersing lichen symbionts in the genus Lepraria (Lecanorales: Stereocaulaceae)[END_REF][START_REF] Rowan | Intraspecific diversity and ecological zonation in coralalgal symbiosis[END_REF]. In such case, the horizontal acquisition, which could be hazardous because of the uncertainty in partner occurrence, is nevertheless an efficient way to obtain new or locally better adapted symbionts.

Host symbiont acquisition is not necessarily exclusive by means of vertical or horizontal mechanisms: a situation named by [START_REF] Ebert | The Epidemiology and Evolution of Symbionts with Mixed-Mode Transmission[END_REF] as the mixed-mode transmission of symbionts. This mixed-mode has already been observed in cases where host-symbiont phylogenies were not totally coupled due to occasional horizontal acquisition in holobionts alternatively classified as vertical transmitters [START_REF] Stewart | Lateral Symbiont Acquisition in a Maternally Transmitted Chemosynthetic Clam Endosymbiosis[END_REF]. The horizontal acquisition and the mixed-mode transmission imply some flexibility in the symbiotic association, an important point for the response to environmental variations. Indeed, if the whole holobiont is the unit of selection in a strict vertical transmission, in a mixed-mode transmission, the selective pressure could act on either hosts or symbionts as well as the holobiont level [START_REF] Theis | Getting the Hologenome Concept Right: An Eco-Evolutionary Framework for Hosts and Their Microbiomes[END_REF].

In marine Cnidaria (including hydrozoans, corals, gorgonians, and sea anemones), the animal host can be associated with eukaryotic mutualistic photosynthetic symbionts (see for review [START_REF] Clavijo | Polymorphic adaptations in metazoans to establish and maintain photosymbioses[END_REF]. These symbionts play key roles in the physiology and ecology of their hosts. Among Anthozoa, the main photosynthetic symbionts are dinoflagellates belonging to the Symbiodiniaceae family [START_REF] Lajeunesse | Systematic Revision of Symbiodiniaceae Highlights the Antiquity and Diversity of Coral Endosymbionts[END_REF]. Located intracellularly, these symbionts translocate most of the organic compounds they produce by photosynthesis to their host cells [START_REF] Davy | Cell Biology of Cnidarian-Dinoflagellate Symbiosis[END_REF]. The cnidariandinoflagellate symbiosis is highly diverse with the possible occurrence of multiple dinoflagellate genotypes within a given host. This diversity can be partially explained by diverse ecological factors which can lead to the presence of different dinoflagellate genotypes (see Parkinson & Baums, 2014 for review). In addition, bacterial symbionts are involved in the good working of the cnidarian-dinoflagellate holobiont (e.g. for nitrogen cycle) and contribute to the holobiont diversity [START_REF] Rohwer | Diversity and distribution of coral-associated bacteria[END_REF][START_REF] Suggett | Symbiotic Dinoflagellate Functional Diversity Mediates Coral Survival under Ecological Crisis[END_REF]. The dinoflagellate symbionts can be acquired by the cnidarian hosts via vertical transmission or horizontal acquisition (at larval or adult stages).

Why vertical or horizontal symbiont acquisition modes are associated with a given cnidarian species is still not understood, but in scleractinian corals, only 25% of spawning species (i.e. those that emit their gametes directly into the ocean) have symbionts in their eggs (i.e., vertical transmission) and belong to only a few specific genera (Montipora, Acropora, Pocillopora, and Porites;[START_REF] Baird | Systematic and Biogeographical Patterns in the Reproductive Biology of Scleractinian Corals[END_REF]. In contrast, temperate symbiotic sea a m a ma ly cal "cl acq y m " [START_REF] Muller-Parker | Temperate and tropical algal-sea anemone symbioses[END_REF].

Although horizontal acquisition in adult cnidarians has been poorly studied in the field (except after post-bleaching episodes and mainly for scleractinian corals [START_REF] Boulotte | Exploring the Symbiodinium rare biosphere provides evidence for symbiont switching in reef-building corals[END_REF][START_REF] Silverstein | Change in algal symbiont communities after bleaching, not prior heat exposure, increases heat tolerance of reef corals[END_REF]), [START_REF] Byler | Multiple Symbiont Acquisition Strategies as an Adaptive Mechanism in the Coral Stylophora pistillata[END_REF] demonstrated that some species, such as Stylophora pistillata, display mixed-mode transmission. The flexibility that horizontal transmission may offer to the host is particularly important considering the functional diversity of symbionts at genus or species level [START_REF] Brading | Differential effects of ocean acidification on growth and photosynthesis among phylotypes of Symbiodinium (Dinophyceae)[END_REF][START_REF] Reynolds | Enhanced photoprotection pathways in symbiotic dinoflagellates of shallow-water corals and other cnidarians[END_REF][START_REF] Roberty | Differential antioxidant response between two Symbiodinium species from contrasting environments[END_REF][START_REF] Suggett | Photosynthesis and Production of Hydrogen Peroxide by Symbiodinium (Pyrrhophyta) Phylotypes with Different Thermal Tolerances1[END_REF][START_REF] Tchernov | Membrane lipids of symbiotic algae are diagnostic of sensitivity to thermal bleaching in corals[END_REF], and may also apply to associated bacteria [START_REF] Neave | Differential specificity between closely related corals and abundant Endozoicomonas endosymbionts across global scales[END_REF]. Even within Symbiodiniaceae species, functional genetic differentiation has been measured between two symbiont populations [START_REF] Howells | Corals in the hottest reefs in the world exhibit symbiont fidelity not flexibility[END_REF][START_REF] Howells | Coral thermal tolerance shaped by local adaptation of photosymbionts[END_REF].

Using ITS2 variability, we previously investigated the genetic diversity of the symbionts hosted by Anemonia viridis, a species of symbiotic sea anemone widely distributed in Europe, from the Mediterranean Sea to the North Sea and the west coast of Ireland, in which different cryptic species (independent lineages) have been recently identified [START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF]. In A. viridis, a unique Symbiodiniaceae clade, 'Temperate A (LaJeunesse et al., 2018;

sensu [START_REF] Savage | Molecular diversity of symbiotic algae at the latitudinal margins of their distribution: Dinoflagellates of the genus Symbiodinium in corals and sea anemones[END_REF][START_REF] Visram | Molecular diversity of symbiotic algae of the genus Symbiodinium (Zooxanthellae) in cnidarians of the Mediterranean Sea[END_REF], has been identified, showing a clear genetic diversity within the clade (Casado-Amezúa, [START_REF] Casado-Amezúa | New insights into the genetic diversity of zooxanthellae in Mediterranean anthozoans[END_REF][START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF]. This symbiont diversity was observed with respect to the geographical origin and the host lineage [START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF].

Specifically, in this latter study, we revealed the occurrence of different assemblages of Temperate A Symbiodiniaceae in at least two lineages of A. viridis (English Channel and Mediterranean Sea). However, as the two lineages were found in two far geographically distinct areas, it was unclear whether the symbiotic assemblage was a result of response to the environment vs. a result of genetic drift in a closed system.

A. viridis can reproduce both sexually, emitting gametes in the water column, and asexually, by longitudinal fission of the gastrovascular column (from the foot to the mouth, [START_REF] Shick | A Functional Biology of Sea Anemones[END_REF]. Whatever the reproduction mode, this sea anemone is thought to inherit its symbionts by vertical transmission either by (i) maternal transmission through oocytes [START_REF] Schäfer | Fortpflanzung und Entwicklung von Anemonia sulcata (Anthozoa, Actiniaria). I: Fortpflanzungszyklus und Struktur der Oocyten vor und nach der Besamung[END_REF][START_REF] Utrilla | Reproduction of the anthozoan Anemonia sulcata (Pennant, 1777) in southern Spain: From asexual reproduction to putative maternal care[END_REF] with a mean of 300 symbionts per oocytes (personal observations) or (ii) by partitioning of the symbiont p p la f "m " b w w " a g " a a m . The vertical transmission guarantees the continuity of this important association, as the A. viridis holobiont has the potential to be completely autotrophic [START_REF] Davy | Carbon budgets in temperate anthozoandinoflagellate symbioses[END_REF][START_REF] Muller-Parker | Temperate and tropical algal-sea anemone symbioses[END_REF]. Even if symbiont diversity has been studied in A.viridis [START_REF] Casado-Amezúa | New insights into the genetic diversity of zooxanthellae in Mediterranean anthozoans[END_REF][START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF][START_REF] Savage | Molecular diversity of symbiotic algae at the latitudinal margins of their distribution: Dinoflagellates of the genus Symbiodinium in corals and sea anemones[END_REF][START_REF] Suggett | Sea anemones may thrive in a high CO2 world[END_REF][START_REF] Visram | Molecular diversity of symbiotic algae of the genus Symbiodinium (Zooxanthellae) in cnidarians of the Mediterranean Sea[END_REF], no study of the dynamics of symbiont populations at the adult stage has been conducted so far, although it would fill an important gap in the understanding of the plastic and adaptive capacities of this species. Indeed, many studies have highlighted its strong capacity for resistance to extreme and disturbed environments [START_REF] Richier | Symbiosis-induced adaptation to oxidative stress[END_REF][START_REF] Suggett | Sea anemones may thrive in a high CO2 world[END_REF][START_REF] Ventura | Resilience to ocean acidification: Decreased carbonic anhydrase activity in sea anemones under high pCO2 conditions[END_REF]. In this context, the occasional occurrence of horizontal symbiont acquisition on top of the vertical transmission may provide a mechanism for adaptation, adding to the intergenerational stability of the association the capacity to acquire locally better adapted symbiont genotype [START_REF] Parkinson | The extended phenotypes of marine symbioses: Ecological and evolutionary consequences of intraspecific genetic diversity in coralalgal associations[END_REF].

To determine whether a closed vertical mode of symbiont transmission occurs in Anemonia viridis or if concurrent and complementary horizontal symbiont acquisition within a mixedmode transmission is also possible in this species, we investigated the strength of the link between host genetic diversity and symbiont assemblage. A more extensive sampling of the North Western Mediterranean Sea than in [START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF] was used. We first delimitated reproductively isolated lineages within the A. viridis species complex (using RAD sequencing markers), and then compared their symbiotic contents (based on NGS acquired nuclear ITS2 and microsatellite diversity) to establish the host-symbiont genetic specificity. We also used the pervasive asexual reproduction occurring within these A. viridis lineages to estimate the similarity of the symbiont communities within and among different sea anemone clonal genotypes (as identified with RAD sequencing markers) sampled in distinct geographical locations. We hypothesized that, in strictly closed vertical transmission, ramets (i.e.

representatives of the same clonal genotype) should be more similar in their symbiont content than two sea anemones from different genets (different clonal genotypes), irrespective of their geographical location. In contrast, the occurrence of horizontal symbiont acquisition would be revealed by a strong geographical effect independent of the lineage and/or genotype of the host.

MATERIAL AND METHOD

Mediterranean Sea anemones

We sampled 237 individuals of Anemonia viridis [START_REF] Forskål | Descriptiones animalium, avium, amphibiorum, piscium, insectorum, vermium; quae in itinere orientali observait[END_REF] var smaragdina (characterised by pink tips and a green fluorescence, [START_REF] Mallien | Conspicuous morphological differentiation without speciation in Anemonia viridis (Cnidaria, Actiniaria)[END_REF][START_REF] Wiedenmann | The morphs of Anemonia aff. Sulcata (Cnidaria, Anthozoa) in particular consideration of the ectodermal pigments[END_REF] from five locations in the Mediterranean eastern basin: Banyuls [Ban, 50 sea anemones sampled], Barcarès [Bar,44], Sète [Tha,48], Marseille [Mar,49] and Scandola [Sca,46] (see figure 1 for GPS coordinates). For each sea anemone, a dozen tentacles were cut, fixed in 70% ethanol and preserved at -80°C until DNA extraction. Total DNA was extracted as in [START_REF] Mallien | Conspicuous morphological differentiation without speciation in Anemonia viridis (Cnidaria, Actiniaria)[END_REF] f ll w g a m f " al g " p c l [START_REF] Miller | A simple salting out procedure for extracting DNA from human nucleated cells[END_REF].

Host genotyping

RAD sequencing

The ad hoc libraries were prepared following [START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF]. Sequencing was performed on an Illumina HiSeq 2000 (100-bp single read format) at the Montpellier-GenomiX platform (MGX, Montpellier, France).Raw sequence reads were demultiplexed, filtered and clustered with iPyRAD v0.3041 [START_REF] Eaton | PyRAD: Assembly of de novo RADseq loci for phylogenetic analyses[END_REF] with a clustering threshold set to 90%.

Host datasets

From the iPyRAD outputs, we derived a base dataset by filtering on missing data: the individuals with a low proportion of shared loci were identified using the R package RADami v1.1-2 [START_REF] Hipp | RADami: R package for phylogenetic analysis of RADseq data[END_REF] and removed using VCFTools v0.1.15 [START_REF] Danecek | The variant call format and VCFtools[END_REF]; the loci were then filtered to only keep loci present in at least 95% of individuals (197 individuals and 25603 SNPs) (figure 2).

We then filtered the base dataset to keep only biallelic loci with a minimum allele frequency (maf) of 1% This new dataset was called N197 dataset with 197 individuals and 5299 SNPs.

This dataset was used to detect clonality. By selecting only the 148 individuals for which the symbionts could be genotyped (see below), we defined the N148 dataset. This dataset was used for host/symbiont differentiations comparison.

After clonality detection, we were able to filter the individuals from the base dataset by keeping only one ramet per genet, and we only kept the biallelic loci with a maf of 1%. THis constituted the N65 dataset with 65 individuals and 6498 SNPs that was used for host genetic clustering and lineages identification.

The N65 dataset was further filtered to keep 14 individuals only (one individual per location and genetic cluster), and for the SNPs without missing data, with a maf of 10% and a linkage disequilibrium pruning threshold (LD) of 0.2 (package SNPRelate v1.11.2, Zheng, 2012) which constituted the BFD dataset of 1197 SNPs. This dataset was used for the test of species hypotheses by BFD* (see below).

Detection of clonality

Even if no anemones did share the same multi locus genotypes, we repeatedly detected similar genotypes in the sampling. These genotypes were likely a consequence of asexual reproduction events and have been grouped into multilocus lineages called genets (see [START_REF] Harper | Population biology of plants[END_REF] in this paper (the representatives of each genet being referred to as ramets).

To formally delineate these genets and associate individual samples to them, we computed between individual pairwise distances for all individuals in the N197 dataset (figure 2) with the R package RClone v1.0.2 (Ba ll l, S ck l, & a -Ha , 2016). The observed distributions of pairwise genetic distances thus obtained were compared to simulated distributions of the expected genetic distances under the assumption of purely sexual reproduction (with or without selfing) computed with RClone (using the functions genet_dist and genet_dist_sim) [START_REF] Bailleul | RClone: A package to identify MultiLocus Clonal Lineages and handle clonal data sets in R[END_REF].) This allowed us to identify a gap in distance distribution and define a threshold distance value separating distances within genets from distances among genets. We could thus associate individuals to genets using the function mlg.filter (w "fa _ g b " alg m) from the R package poppr v2.5.0 [START_REF] Kamvar | Poppr: An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction[END_REF]. We then calculated the clonal richness (R) considering the number of genets as the number of different genotypes over the total number of sampled sea anemones (globally, per location and per host genetic cluster).

Structure of host genetic diversity

To identify putative species within A.viridis individuals from the Western Mediterranean basin, we used the N65 dataset (figure 2) and we computed a species tree using SVDquartets [START_REF] Chifman | Quartet Inference from SNP Data Under the Coalescent Model[END_REF] implemented in PAUP* v4.0A165 [START_REF] Swofford | Phylogenetic analysis using parsimony (paup) (Version 4)[END_REF][START_REF] Swofford | PAUP*: Phylogenetic Analysis Using Parsimony (and other methods) 4[END_REF] with an exhaustive quartet sampling and performed a sNMF quantitative hierarchical genetic clustering with the LEA (v2.2.0) R package [START_REF] Frichot | LEA: An R package for landscape and ecological association studies[END_REF], based on individual ancestry coefficients. The optimal number of ancestral clusters, from which the individuals in the dataset could be traced back, was determined using the entropy criterion.

The results were visualized using a STRUCTURE-like representation showing, for each individual, the relative contribution of each ancestral cluster to its multilocus genotype. The ramets not kept in the N65 dataset were attributed by extension to the main cluster of their genet.

To test whether the sNMF clusters and SVD clades were independent species, we compared the likelihood of the species hypothesis they constitute, to the likelihood of a primary species hypothesis based on the geographical locations of the individuals using Bayes factor delimitation with genomic data (BFD*; [START_REF] Leaché | Species Delimitation using Genome-Wide SNP Data[END_REF]. In order to manage computing times, these analyses were performed using the BFD dataset (figure 2).

The BFD* analyses were performed in BEAST2 [START_REF] Bouckaert | BEAST 2: A Software Platform for Bayesian Evolutionary Analysis[END_REF] with 48 path sampler sets of 100000 MCMC repetitions, to sample in 800000 MCMC iterations of SNAPP [START_REF] Bryant | Inferring Species Trees Directly from Biallelic Genetic Markers: Bypassing Gene Trees in a Full Coalescent Analysis[END_REF]. The number of SNAPP iterations was chosen to ensure proper convergence of the species tree inference. To ensure the convergence of the computations while maintaining amenable computing times, the BFD* analyses were repeated on two partitions of the 1197 SNPs of the BFD dataset. The first analysis was run on the 598 first SNPs and the second one the other 599 SNPs of the dataset.

W g F-statistics among and within populations were computed using the R packages StAMPP v1.5.1 [START_REF] Pembleton | StAMPP: An R package for calculation of genetic differentiation and structure of mixed-ploidy level populations[END_REF] and hierfstat v0.04-22 [START_REF] Goudet | Hierfstat: Estimation and tests of hierarchical F-statistics[END_REF].

Symbiont genotyping

Amplicon sequencing

The symbiont populations were genotyped using the nuclear ribosomal DNA internal transcribed spacer 2 (ITS2). We used the primers from [START_REF] Stat | Specificity in communities of Symbiodinium in corals from Johnston Atoll[END_REF] f wa D (5 -GTGAATTGCAGAACTCCGTG -) a 2 2 (5 -CCTCCGCTTACTTATATGCTT -).

The symbiont populations were also genotyped using four microsatellites: three designed from an Expressed Sequence Tag library of A. viridis holobionts [START_REF] Sabourault | Comprehensive EST analysis of the symbiotic sea anemone, Anemonia viridis[END_REF] and one from the study of [START_REF] Catanéo | Development of microsatellite loci from a Mediterranean symbiotic gorgonian, Eunicella singularis (Cnidaria, Anthozoa), and its temperate-A Symbiodinium[END_REF] (see also [START_REF] Abdoullaye | Permanent Genetic Resources added to Molecular Ecology Resources Database 1 August 2009-30[END_REF] about the Temperate A symbionts from the Mediterranean gorgonian Eunicella singularis (with the NCBI accession number FJ917549.1, table 1). These markers are particularly useful to genotype populations instead of species [START_REF] Coffroth | Genetic Diversity of Symbiotic Dinoflagellates in the Genus Symbiodinium[END_REF] as they are more informative of recent evolutionary events [START_REF] Cornuet | Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v.1.0)[END_REF].

PCR amplifications of the symbiont markers were performed by Access Array (FLUIDIGM) PCR and next generation sequencing of the amplicons were obtained by MiSeq Illumina Technology (2x250bp paired-end) at the Brain and Spine Institute (ICM, Paris, France).

Genotyping methods

Because a given sea anemone can host 3 to 20 million symbiont cells per milligram of protein (values from field and lab experiments, [START_REF] Ventura | Resilience to ocean acidification: Decreased carbonic anhydrase activity in sea anemones under high pCO2 conditions[END_REF], the sequences obtained from each anemone were in fact obtained from a pooled sample of symbionts. Their genotyping was therefore not straightforward and required additional analysis steps.

In this context, we performed an analysis of the symbiont diversity within the A. viridis using ITS2 and microsatellite sequence variants. The analysis of ITS2 sequence variant assemblages is an effective method for comparing Symbiodiniaceae populations despite the marker s multicopy nature, which complicates data interpretation (see [START_REF] Hume | SymPortal: A novel analytical framework and platform for coral algal symbiont next-generation sequencing ITS2 profiling[END_REF]. The ITS2 sequences were processed and identified using the SymPortal pipeline [START_REF] Hume | didillysquat/SymPortal_framework: V0[END_REF][START_REF] Hume | SymPortal: A novel analytical framework and platform for coral algal symbiont next-generation sequencing ITS2 profiling[END_REF] based on a minimum entropy decomposition algorithm, as in [START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF].

For the microsatellite markers, the forward and reverse reads were assembled, then trimmed and tagged by microsatellite locus using OBITools scripts [START_REF] Boyer | OBITOOLS: A UNIX -inspired software package for DNA metabarcoding[END_REF]. To filter background noise and to identify reads with actual microsatellite repeats, we used the MicNeSs python script [START_REF] Suez | MicNeSs: Genotyping microsatellite loci from a collection of (NGS) reads[END_REF] with a tolerance of one substitution in the microsatellite motif. Because MicNeSs was designed for diploid organisms, and because we did not want to bias the genetic diversity by forcing a maximum of two microsatellite alleles per within host symbiont population, we chose to consider for each individual, instead of the genotype, the raw observed distribution of the reads containing a given microsatellite locus.

We kept for further analysis only the sea anemones with amplification for at least three out of the 4 symbiont microsatellite markers used. Therefore, as we did not attempt to ascertain any of the microsatellite alleles present in the read collections, this raw observed distribution of the reads is more a fingerprint of the symbiont community composition than a proper collection of individual genotypes. Nevertheless, it can still be used to measure the level of differentiation in symbiont content among different sea anemones.

The symbiont populations were not successfully amplified from all the sampled individuals, and only 148 sea anemones could be genotyped for both host and symbiont microsatellite markers. These individuals composed the N148 dataset (table S1).

Symbiont differentiation

We compared the different compositions of symbiont populations using over the N148 dataset, for the ITS2, the presence/absence of sequence variants among the hosts, and for the microsatellite markers, the observed relative frequency distributions of the reads. For these reasons, the differentiation between hosts in terms of symbiont community composition was expressed both a a Jacca a c (f ITS2 sequence variant assemblage considering the presence and absence of the variants) a a Ma a a distance (for the distributions of microsatellite markers). The distances among hosts/symbiont community compositions were visualized with Neighbour Joining trees with R package ape v5.3 [START_REF] Paradis | APE: Analyses of Phylogenetics and Evolution in R language[END_REF]) with 1000 bootstraps. We measured the correlation between symbiont p p la a c a g c a c by Ma l tests. We also tested by PERMANOVA with the R package vegan v2.5-5 [START_REF] Oksanen | vegan: Community Ecology Package. (Version R package version 2.4-0[END_REF] the differences in symbiont community composition between sea anemones depending on

RESULTS

Clonality in A. viridis var. smaragdina populations

We detected clonality in the dataset with a gap in the distribution of pairwise genetic distances among A. viridis individuals (figure S1), and hence identified 65 genets over the 197 sea anemones sequenced, with 37 genets considered as unique (one ramet) and 28 genets (here called clonal genets) with at least two representatives (i.e. at least two ramets produced from asexual reproduction events, see table S1). These 28 clonal genets could be found in different distant sampling locations (table S1). We hence estimated a global clonal richness of 0.33.

Host genetic structure

To identify independent genetic pools, we performed a sNMF genetic clustering. Based on minimum cross-entropy value (figure S2), we determined K=4 as the optimal number of ancestral genetic clusters among the 65 genets (N65 dataset). Among these four different genetic clusters, we detected admixture, more specifically two of them seemed to be more admixed (K1 and K2, figure 1.A): more than a third of the sea anemones in K2 where not purely descended from this cluster, and had at least 15% of their loci from K1 and, reciprocally, more than half of the sea anemones in K1 admixed with the K2 cluster (figure 1.A and table S1). All these clusters were found in at least two different sampling locations (figure 1.A and 1.B) with K2 being the most common cluster among the Mediterranean sampled populations. The admixture K1/K2 was correlated with the lowest F ST value of 0.119 compared with the other pairwise differentiations (table S2). The K4 cluster appeared as the most differentiated genetic cluster in the Mediterranean Sea with F ST values to other clusters ranging from 0.242 to 0.350

We also used SVDquartets to produce a coalescent taxonomic partitioning from the N65 dataset (figure S3). The SVDquartets analyses showed three main clades (clades 1, 3 and 4) and an unresolved group (clade 2). These clades are partially congruent with the clusters identified by the sNMF analysis: the clades 3 and 4 corresponded to the clusters K3 and K4 respectively, and clade 1 corresponded to a part of the cluster K1 (table S1).

To know which one of the partitions given by sNMF clustering (hypothesis 2 in table 2) or SVDquartet clades (hypothesis 3 in table 2) was the most likely secondary species hypothesis (SSH), we compared them to a null hypothesis (a primary species hypothesis, hypothesis 1) based on geographical origin, using BFD*. In addition to these species hypotheses we considered another one bringing K1 and K2 together (hypothesis 4) because of the admixture between K1 and K2 (figure 1.A). In addition, using two different sets of host SNPs (see Material and Method section), the best SSH, using Bayes Factor comparisons, was the sNMF clusters partition (table 2). We hence identified four different host genetic lineages in A. viridis Mediterranean populations.

Symbiont differentiation

The two types of genetic markers (ITS2 and microsatellites) used to genotype the symbiont populations from 148 A. viridis individuals (N148 dataset) revealed an overall similar picture However, we detected with both types of markers a global effect of the sampling location on the symbiont composition with PERMANOVA analyses (table 3, table S3 and table S4) with similar symbiont content among host clusters at the same location (table S4). Nevertheless, with microsatellite markers, a weak effect of the host genetic cluster was revealed (table 3) mainly due to a differentiation between K2 and K4 occurring at Scandola (table S4).

Symbiont community composition was more similar among sea anemones from the same location (intra-and inter-lineages Ma a a average distances of D = 2.186 and 2.502 respectively, figure 5.A), than within the same host lineage at different locations (intralineage inter-location, D = 2.951) and finally non-related individuals (inter-lineage interlocation, D = 3.063). At a finer genetic scale, we used the clonal host genets found in different locations to evaluate the differentiation of the symbiont community composition within a genet. Six of the identified genets were distributed in at least two locations with at least two ramets within a location (table S5). If vertical transmission was the only mode of transmission of symbionts, we expected that the weakest differentiation would be observed within a genet independently of the geographical origin. However, the differences in the symbiont community composition among individuals from the same location were significantly lower (intra-and inter-genets: D = 1.331, 2.782 resp.) than the -differences among individuals from different locations (intra-and inter-genets: D = 3.223, 3.027 resp., figure 5.B). Hence the genets from the same location had similar symbiont populations.

DISCUSSION

In this study, we questioned whether the symbiont diversity observed in A. viridis adult individuals was the result of closed vertical transmission only or whether secondary horizontal acquisition and/or local selection could also occur and modify the symbiont population in respect to the environment. We showed that sea anemones belonging to the At all locations, lineages were found in sympatry with different degrees of admixture among them. These admixture signals probably explained the incongruence between quantitative genetic clustering (sNMF) and coalescent analysis (SVDquartet), with for example, seven admixed individuals of the cluster K1 belonging to the clade 2 (table S1). As previously found, Banyuls was the most genetically diverse location (with the presence of the four lineages), supporting our supposition of the presence of a strong contact zone between genetic pools in this Mediterranean zone [START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF]. In addition, we now describe an ubiquitous presence of the Med2 (K2) lineage at all sites, suggesting a strong ability for this lineage to occupy different habitats and maybe great dispersal capacities as well. In all lineages and sites, we detected an important clonal richness, which confirmed that A. viridis can deploy genets with a great number of ramets [START_REF] Mallien | Conspicuous morphological differentiation without speciation in Anemonia viridis (Cnidaria, Actiniaria)[END_REF][START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF][START_REF] Sauer | Strategien zeitlicher und räumlicher Einnischung[END_REF][START_REF] Sauer | Alloimmune memory for glycoproteid recognition molecules in sea anemones competing for space[END_REF][START_REF] Wiedenmann | The relationship between different morphs of Anemonia aff. Sulcata evaluated by DNA fingerprinting (Anthozoa, Actinaria)[END_REF], even if, within this sampling, the Med4 (K4) lineage exhibited a much lower clonality. It is interesting to note that the Med2 (K2) lineage showed varying levels of clonality from one site to the other (R varied from 0.05 in Scandola to 0.9 in Sète), in a similar way to what had been observed in another sea anemone, Metridium senile, with different clonality rates in different locations along the N W Atlantic coast (from New York to Nova Scotia) [START_REF] Hoffmann | Variation in contribution of asexual reproduction to the genetic structure of populations of the sea anemone Metridium senile[END_REF]. More surprisingly, the ramets of the Med1 (K1), Med2 (K2) or Med3 (K3) genets were not always located at the same site and some genets were comprised of individuals present in distant sites (tables 4 and S1), suggesting again a potentially high capacity of migration/dispersion in this species, as already alluded to in [START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF]. To disperse, adult sea anemones may detach and move up to several meters away [START_REF] Ayre | Effects of environment and population density on the sea anemone Actinia tenebrosa[END_REF][START_REF] Mackie | Locomotion, Flotation, and Dispersal[END_REF][START_REF] Shick | A Functional Biology of Sea Anemones[END_REF], but long-distance migration could likely occur by rafting attached to f ly f g a al (a w ) a f c al (a p b m ) b a [START_REF] Jackson | Modes of dispersal of clonal benthic invertebrates: Consequences for species' distributions and genetic structure of local populations[END_REF].

The dispersal of asexual larvae as it occurs in some corals could also explain the distribution of clonal genets [START_REF] Gélin | Superclone Expansion, Long-Distance Clonal Dispersal and Local Genetic Structuring in the Coral Pocillopora damicornis Type β in Reunion Island, South Western Indian Ocean[END_REF] but it was not observed yet in A. viridis [START_REF] Utrilla | Reproduction of the anthozoan Anemonia sulcata (Pennant, 1777) in southern Spain: From asexual reproduction to putative maternal care[END_REF]. In fine, even if established sea anemone clones have been considered as highly locally adapted [START_REF] Ayre | Localized Adaptation of Clones of the Sea Anemone Actinia Tenebrosa[END_REF], migration and successful new establishment of ramets are not rare in A. viridis.

Hosts from the same genetic lineage harbour different symbiont communities

With both marker types (ITS2 and microsatellites), no correlation was revealed between the genetic distances among symbionts and the genetic distances among their hosts (figure 3 and 4) suggesting that the closed vertical transmission is not the only mechanism of symbiont acquisition in adult hosts. Our results also showed that the symbiont ITS2 diversity was only, and strongly, linked to the geographical site of sampling. Microsatellite diversity was strongly linked to the sampling location as well (table 3). The microsatellite markers also revealed a weak correlation of symbiont content with the host genetic lineage (table 3), which would be expected considering that symbionts are transferred through the eggs in this species [START_REF] Schäfer | Fortpflanzung und Entwicklung von Anemonia sulcata (Anthozoa, Actiniaria). I: Fortpflanzungszyklus und Struktur der Oocyten vor und nach der Besamung[END_REF][START_REF] Utrilla | Reproduction of the anthozoan Anemonia sulcata (Pennant, 1777) in southern Spain: From asexual reproduction to putative maternal care[END_REF]. This finding was summarized in the figure 6. It should also be noted that only one of the within site comparisons of symbiont microsatellite size distribution between host lineages was significant; at Scandola, where a strong differentiation between the lineages of Med2 (K2) and Med4 (K4) occurred (table S4). This was likely due to a large effect of host clonality on symbiont differentiation (preponderance of the genet R155 at this site, figure 3B and table S5). The presence of this large genet is surprising and could reflect a sampling artefact but could also be the consequence of the dominance of a highly aggressive clone, a phenomenon already identified in sea anemones by [START_REF] Ayre | Aggression, Habituation, and Clonal Coexistence in the Sea Anemone Anthopleura elegantissima[END_REF].

Hosts from the same genet harbour different symbiont communities

If symbionts are vertically transmitted, the individuals belonging to the same genet, present or not at the same site, should be dominated by the same symbiont community. We observed however that individuals from the same genet, but from different locations, were not more similar than individuals from different genets (figure 5.B, intra-genet inter-location vs. inter-genet intra-or inter-location). Considering the high probability that the revealed microsatellite diversity is neutral, this result is strongly in favour of the occurrence of horizontal acquisition of symbionts at the adult stage in A. viridis, and thus against the hypothesis of an entirely closed vertical transmission of symbionts.

A mixed-mode symbiont transmission in Anemonia viridis

This result highlights the potential for a plastic response of A. viridis lineages to the environment by modifying their symbiont community composition. Previous studies stated that no genetic changes in the symbiont community composition occurred in A. viridis in response to environmental change between subtidal and intertidal habitats [START_REF] Bythell | Algal genotype and photoacclimatory responses of the symbiotic alga Symbiodinium in natural populations of the sea anemone Anemonia viridis[END_REF] or different pH conditions [START_REF] Borell | Increasing pCO2 correlates with low concentrations of intracellular dimethylsulfoniopropionate in the sea anemone Anemonia viridis[END_REF][START_REF] Suggett | Sea anemones may thrive in a high CO2 world[END_REF]; or after living for a decade in proximity to cnidarians containing mostly clade C (genus Cladocopium now) symbionts [START_REF] Hartle-Mougiou | Diversity of zooxanthellae from corals and sea anemones after long-term aquarium culture[END_REF].

These studies used ribosomal sequences as a genetic marker, but the flexibility has never been investigated at the symbiont diversity scale using marker analysis by NGS (ITS2 or successive bottlenecks (at the production of the ovule). A mixed-mode of symbiont acquisition with episodes of horizontal acquisition from local environment is the most parsimonious explanation of the distribution of symbiotic diversity we observed. Indeed, a very similar pattern of variability in symbiont population among host ramets had also been observed in Acropora palmata, a coral species that only acquires its symbiont from the environment at each generation [START_REF] Baums | New insights into the dynamics between reef corals and their associated dinoflagellate endosymbionts from population genetic studies[END_REF].

Additional horizontal symbiont acquisition in hosts with maternal vertical symbiont transmission has already been suggested in some corals such as Pocillopora meandrina, Pocillopora damicornis Stylophora pistillata and Seriatopora hystrix [START_REF] Boulotte | Exploring the Symbiodinium rare biosphere provides evidence for symbiont switching in reef-building corals[END_REF][START_REF] Byler | Multiple Symbiont Acquisition Strategies as an Adaptive Mechanism in the Coral Stylophora pistillata[END_REF][START_REF] Magalon | High genetic diversity of the symbiotic dinoflagellates in the coral Pocillopora meandrina from the South Pacific[END_REF][START_REF] Quigley | Relationship between Acropora millepora juvenile fluorescence and composition of newly established Symbiodinium assemblage[END_REF]. In addition, analysing specifically within-colony symbiont diversity in Pocillopora sp., [START_REF] Pettay | Genotypic diversity and spatial-temporal distribution of Symbiodinium clones in an abundant reef coral[END_REF], also strongly suggested that switching may occur between clonal host lineages and that this species can mix vertical and horizontal symbiont acquisition.

Mixed-mode symbiont acquisition is frequent in symbiotic cnidarians and supports the conclusions of [START_REF] Ebert | The Epidemiology and Evolution of Symbionts with Mixed-Mode Transmission[END_REF] that symbionts with mixed-mode transmission may be the most common type of symbionts in animals. . W a k R. f c f c c a acc computing facilities. We also thank P-A.Gagnaire and L.Sromek for the RADseq individual tags; P.Barbry from IPMC who provides access to the CovarisB sonicator and V.Magnone and M-J.Arguel for their help during its utilization. We also thank C.Fauvelot for her attentive reading and her suggestions. The authors thank the anonymous reviewers for their constructive comments.

ACKNOWLEDGEMENTS

TABLES

Table 1: Detail of the symbiont microsatellite markers sequenced by NGS Table 2: Test of A. viridis species hypotheses by BFD*. The four species hypotheses tested were (1) a geographic partitioning (the primary species hypothesis), (2) the sNMF clustering, (3) the SVDquartets partitioning and (4) an intermediate SVDquartets partitioning considering the K1 and K2 admixture. (MLE = maximum likelihood estimator, BF = bayes factor). In bold the secondary species hypothesis retained with the lowest MLE and the degree of significance when compared with the other species hypotheses. Figure 6: Schematic illustration of the host and symbiont differentiation among sites according to our results. The colours of the sea anemones correspond to the host lineages or host lineages, the symbols inside are the symbiont communities.

  (i) the host genetic clusters, or (ii) the sampling location. When these factors had significant effects on the distribution of the genetic diversity of the symbionts, we then computed pairwise PERMANOVAs between each level of the given factor with an FDR correction -location genets. Following the approach of[START_REF] Baums | New insights into the dynamics between reef corals and their associated dinoflagellate endosymbionts from population genetic studies[END_REF] in a coral species, we took advantage of asexual reproduction in A. viridis to test for the tightness of the closed vertical symbiont transmission in A. viridis, by comparing the symbiont diversity among ramets. If symbiont diversity is only the result of mother to offspring transfer, we expected that ramets from the same genet should have a similar symbiont community composition and hence the distances should be shorter compared to genetically distant hosts. The differences were statistically tested with a Kruskal-Wallis test implemented in R (R Core Team, 2018).

  figure 4.A) or with the m c a ll ma k (w a Ma l = 0.06 a a p-value of 0.086,

  microsatellite) or sampling geographically distant A. viridis populations. Such flexibility would be beneficial for an organism living in a diverse and changing environment such as this temperate sea anemone ( f xampl S c & M ll -a k , 2005). However, the mechanism involved in this symbiont flexibility cannot be unambiguously resolved yet as it could still be putatively due to i) a renewal of symbiont content by horizontal acquisition from the environment at the adult stage or; ii) convergent changes in symbiont frequency within large and diverse inherited symbiont populations. In the coral Montipora capitata, comparing ITS2 sequence assemblages in eggs and their parent colonies, Padilla-Gamiño et al. (2012) demonstrated that the eggs harboured a large range of ITS2 sequences and that symbiont assemblages in the eggs were strongly influenced by symbiont composition of the parent colony, which in turn differed and reflected characteristics of their physical environment. In M. capitata, an environmental selection of the inherited symbiont diversity is then more plausible, as was also strongly suggested by Quigley et al. (2019) in another species of the same genus. However, these Montipora species harbour different genera (or species) of Symbiodiniaceae, known to have different ecophysiological profiles. Therefore, the relative proportions of these different symbiont strains within a host could very likely be influenced by local selection. In our case, considering that A. viridis harbours only Temperate A Symbiodiniaceae, to explain the observed pattern of local similarity in symbiont composition among host genets and even lineages by convergent selection without horizontal acquisition, one would have to imagine that each host larva, independently of its mother affiliation to a clone or to a reproductively isolated lineage, should inherit the whole complement of symbiont genetic diversity that would allow such local differentiation. This would mean that adult sea anemones belonging to genetically differentiated lineages still manage to transfer a wide array of microsatellite diversity to their offspring, despite generations of local adaptation and

  g am " a " ( NR-12-ADAP-0016) as a PhD grant to CM and by a doctoral fellowship from the French Ministère de l .Pey and A.Haguenauer for their contribution to the sampling effort. We thank M.Suez and D.Higuet for constructive c

FIGURESFigure 1 :

 1 FIGURESFigure1: sNMF genetic clusters of A. viridis and their geographic distribution in the Mediterranean Sea. A: Bar plot representation of the sNMF analysis of the animal host (N65 dataset). The genetic diversity could originate from 4 ancestral genetic clusters (named K1 to K4, respectively in purple, red, blue and green). The bar plot represents the contribution of each of these clusters to the genetic constitution of each individual in this dataset. The sampling location of each individual is mentioned above the barplot. B: Proportion of each cluster in each sampling location (N148 dataset, tableS1), along with local sample size and clonal richness R. (Banyuls: 2°28 50 N 0 °0 50 E; Ba ca è 42°47'15.828"N 3°2'11.633"E; Sète: °2 2.82 N ° 6 1. E; Ma ll 43°17'44.337"N 5°20'34.348"E; Scandola: 42°21'29.571"N 8°33'19.166"E) *Host genetic lineages correspond to genetic lineages either identified in Porro et al (2019) (Med1 and Med2) and in the present study (Med3 and Med4), see discussion section.

Figure 2 :

 2 Figure 2: Host datasets. Description (number of individuals and SNP markers) of the different datasets used for each analysis. BFD is for Bayes Factor Delimitation and refers to the BFD* analysis to test species hypotheses.

Figure 3 :

 3 Figure 3: Differentiation of A. viridis individuals based on algal symbiont composition (N148 dataset) . NJ ba Jacca a c f ITS2 assemblages among the sea a m ; B. NJ ba Ma a a a c f m c a ll b am g sea anemones. The bootstraps higher than 40 were displayed on the nodes. Colours correspond to the host genetic cluster from the figure 1 (K1: purple, K2: red, K3: blue and K4: green); the symbols correspond to the host genets: the ramets of the same genet shared the same coloured symbol, full circles correspond to the unique genets.

Figure 4 :

 4 Figure 4: Sca pl a Ma l f pa w ost genetic distances against pairwise symbiont distances among A. viridis (N148 dataset): considering, A: the symbiont ITS2 q c a mblag (Jacca a c -w a Ma l = -0.06 and p-value = 0.839); B: the symbiont microsatellite distribu(Ma a a a c -w a Ma l = 0.06 and p-value = 0.086). The dashed line indicates the threshold host distance between intra-genet and inter-genet distances (see figureS1).

Figure 5 :

 5 Figure 5: Intra-and inter-locations differentiation in symbiont microsatellite contents for the different host genetic structures (N148 dataset). A: Comparison of distance distributions within lineages (intra-) or between lineages (inter-) from the same sampling location (intralocations) or between different locations (inter-locations). B: Comparison of distance distributions within genet (intra-) or between genet (inter-) from the same sampling location (intra-locations) or between different locations (inter-locations).The global p-values for the both comparisons (A and B) is < 2.2x10 -16 and tested with Kruskal-Wallis tests. The lowercase letters indicate the statistical differences The boxes are the two first quartiles of distributions, the whiskers are the third and the fourth quartiles. The black horizontal bars (in the boxes) are the medians (with IC95 figured with the notches). The dots correspond to the a ag f Ma a a a c f b a a a g by the

viridis lineages in the Mediterranean Sea As

  we previously identified cryptic species in A. viridis, it was critical to first identify to which host lineage the sampled individuals belonged. After taking clonality into account, we could group the sampled individuals into 4 genetic clusters, recognized by BFD* as four differentiated lineages. As a subset of samples used in this study had already been analysed in[START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF], we could recognize that the K1 and K2 lineages from this study corresponded respectively to the Med1 and Med2 lineages from[START_REF] Porro | The many faced symbiotic snakelocks anemone ( Anemonia viridis , Anthozoa): Host and symbiont genetic differentiation among colour morphs[END_REF].Moreover, having increased the Mediterranean sampling density, we were now able to redefine the distribution of Med1 and Med2, which are no longer single-location lineages, and identify two new lineages, named Med3 and Med4, that corresponded to the clusters K3

same host lineage, but sampled from different locations, harbour different symbiont communities, revealing a mixed-mode symbiont transmission. However, how different are A. viridis host lineages? Multiple A. and K4, respectively (figure 1 and table 2).

Table 3 :

 3 PERMANOVA tests of differentiation of symbiont community compositions based on ITS2 sequences (upper part) or microsatellite distributions (lower part).

	markers	factors	Df	Sums of squares	Mean Sqs	F.Model	R2	Pr(>F)
		Location	4	2.4868	0.6217	2.92662	0.08028	1.00E-05
		Lineage	3	0.6705	0.2235	1.05212	0.02165	0.3437
	ITS2	LocationxLineage	7	1.2646	0.18066	0.85043	0.04083	0.857
		Residuals	125	26.5536	0.21243		0.85725
		Total	139	30.9755			1
		Location	4	196.64	49.16	15.152	0.2829	1.00E-05
		Lineage	3	25.66	8.552	2.636	0.03691	0.00221
	microsatellite	LocationxLineage	7	41.29	5.898	1.818	0.0594	0.00566
		Residuals	133	431.5	3.244		0.62079
		Total	147	695.09			1
	PERMANOVAs were realised over 99999 permutations. The bold values indicate factors with statistically

DATA ACCESSIBILITY

The data will be available on Dryad

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest


Berti Luca

Cemosis

Vincent Chabannes
email: chabannes@math.unistra.fr

Laetitia Giraldi
email: laetitia.giraldi@inria.fr

Christophe Prud'homme
email: prudhomm@math.unistra.fr

MODELING AND FINITE ELEMENT SIMULATION OF MULTI-SPHERE SWIMMERS

Keywords: Stokes in moving domains with articulated rigid bodies, 3-Sphere swimmer, Finite Element

We propose a numerical method for the finite element simulation of micro-swimmers composed of several rigid bodies moving relatively to each other. Three distinct formulations are proposed to impose the relative velocities between the rigid bodies. We validate our model on the three-sphere swimmer, for which analytical results are available.

Introduction

The dynamics of immersed rigid and deformable bodies in low Reynolds number flows has been extensively studied for its applications to suspension phenomena and the motion of biological micro-organisms [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications, Butterworth -Heinemann series in chemical engineering[END_REF][START_REF] Happel | Low Reynolds number hydrodynamics, Mechanics of Fluids and Transport Processes[END_REF].

Different numerical methods are used to study these dynamics, such as the Boundary Element method, which discretizes the boundary integral form of Stokes equations [START_REF] Walker | Boundary behaviours of leishmania mexicana: A hydrodynamic simulation study[END_REF]; Resistive Force Theory, where the hydrodynamical interactions are approximated by suitable coefficients [START_REF] Alouges | Self-propulsion of slender micro-swimmers by curvature control: N-link swimmers[END_REF]; the Finite Element method, where the fluid domain is discretized and fitted [START_REF] Maury | Direct simulations of 2d fluid-particle flows in biperiodic domains[END_REF] or unfitted [START_REF] Glowinski | A distributed lagrange multiplier/fictitious domain method for particulate flows[END_REF] meshes are used for the immersed bodies. In the following, the focus will be on the Finite Element method with body-fitted mesh, where the Arbitrary Lagrangian-Eulerian formalism is used to to solve the motion of the fluid domain. While this choice was dictated by its straightforward extension to non-Newtonian fluids and interactions with immersed elastic bodies, the present article is concerned with the simulation of swimmers composed of multiple rigid bodies that move relatively to each other. The main examples are micro-swimmers composed of spherical bodies, which have been extensively studied via the fundamental solution associated to an immersed rigid sphere [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF][START_REF] Lefebvre-Lepot | A stokesian submarine[END_REF][START_REF] Alouges | Optimally swimming stokesian robots[END_REF]. Theoretical and numerical results make these swimmers ideal benchmarking models.

In this paper we propose a numerical method, based on finite elements and the Arbitrary-Lagrangian Eulerian formulation, that allows the simulation of micro-swimmers composed of rigid parts moving relatively to each other, extending [START_REF] Maury | Direct simulations of 2d fluid-particle flows in biperiodic domains[END_REF] to self-propelled swimmers. We illustrate our method by recovering the motion of the three-sphere swimmer, as reported in [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF].

Mathematical formulation of the problem

In this section we recall the formulation of a generic swimming problem in Stokes flow. Let F 0 ⊆ R d , d = 3, be the initial configuration of the fluid domain and A t : F 0 → F t the smooth function that maps the reference fluid domain onto the domain F t occupied by the fluid at time t. Let S 0 ⊆ R d be the domain occupied by the swimmer at time t = 0 and S t ⊆ R d be the domain it occupies at time t. Let (u, p) be the fluid velocity and pressure defined over F t , µ the viscosity, f the fluid volume forces. We denote by U and ω the translational and angular velocities of the swimmer, by x CM its center of mass and by m and J its volume and geometric inertia tensor. The coupled problem, expressing the kinematic and dynamic coupling of the fluid with the swimmer, reads

               -µ∆u + ∇p = f, in F t , ∇ • u = 0, in F t , u = U + ω × (x -x CM (t)) + u d (t) • A -1 t , on ∂S t , m U = -F f luid , J ω = -M f luid . (1)
where F f luid and M f luid denote the global fluid forces and torques acting on the boundary of the swimmer. The function u d (t) denotes the time derivative of body deformation, and it is the driving mechanism of the whole motion, making the otherwise stationary problem time-dependent.

The map A t (x) is defined as A t (x) = x + φ t (x), where φ t (x) is the solution of the following Laplace problem

∆φ t (x) = 0 in F 0 , φt (x) = U + ω × (x -x CM (t)) + u d (t), on ∂S 0 . ( 2
)
3 The articulated micro-swimmer

The articulated micro-swimmer is composed of n rigid bodies among which a reference body B n is identified. This body is linked to all the other bodies B i , for i ∈ {1 . . . n -1}, by thin and hydrodynamically negligible arms. The length of these links can be changed via "internal motors" that impose a relative speed between the bodies, leading to self-propulsion.

In principle, each B i can have a different shape, but we choose to work with spheres for benchmarking purposes. In fact, the motion of some multi-sphere swimmers can be analytically computed by using the appropriate Green kernel of Stokes equations [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF].

In this section the motion of n independent bodies is first described, by recalling the formulation in [START_REF] Maury | Direct simulations of 2d fluid-particle flows in biperiodic domains[END_REF]. After that, we present the modified formulation that takes into account the relative velocities between the bodies. In the end, we propose few equivalent methods that allow the imposition of the velocity constraints.

Motion of n independent bodies

Using the same notations as before, the problem of n independent bodies moving in a Stokes fluid reads

               -µ∆u + ∇p = f, in F t , ∇ • u = 0, in F t , u = U i + ω i × (x -x CM i (t)), i = 1 . . . n, on ∂B i , m i Ui = -F f luid , i = 1 . . . n, J i ωi = -M f luid , i = 1 . . . n. (3)
In this case u d (t) = 0, as the bodies move independently of each other and no deformation is present. Hence,

A t (x) = x + φ t (x)
, where φ t (x) is the solution of the Laplace problem

∆φ t (x) = 0 in F 0 , φt (x) = U i + ω i × (x -x CM i (t)), on ∂B i . (4)
In this formulation, the body motion is dictated by fluid stresses only. We now address the variational formulation of (3) and look for a solution (u, p,

U i , ω i ) ∈ [H 1 (F t )] d × L 2 (F t ) × [R d ] n × [R d ] n to the weak formulation of the problem. Let (ũ, p, Ũi , ωi ) ∈ [H 1 (F t )] d × L 2 (F t ) × [R d ] n × [R d ] n denote the test functions. The variational formulation reads 2µ Ft D(u) : D(ũ) dx - n i=1 ∂Bi (-pI + 2µD(u)) n • ũ dS - Ft p∇ • ũ dx = Ft f • ũ dx, (5)
Ft p∇ • u dx = 0, (6)
m i Ui • Ũi = - ∂Bi (-pI + 2µD(u)) n • Ũi dS, i = 1 . . . n, (7)
J i ωi • ωi = - ∂Bi (-pI + 2µD(u)) n × (x -x CM i ) • ωi dS, i = 1 . . . n, (8)
where

D(u) = 1 2 (∇u + ∇u T ). Following [8], we choose (ũ, Ũi , ωi ) ∈ [H 1 (F)] d × [R d ] n × [R d ] n satisfying u = U i + ω i × (x -x CM i ) on ∂B i . These test functions form a subspace of [H 1 (F t )] d × [R d ] n × [R d ] n ,
(-pI + 2µD(u)) n • ũ dS = ∂Bi (-pI + 2µD(u)) n • Ũi dS + ∂Bi (-pI + 2µD(u)) n × (x -x CM i ) • ωi dS, (9)
that combines the boundary terms in equations ( 5)-( 7)- [START_REF] Maury | Direct simulations of 2d fluid-particle flows in biperiodic domains[END_REF]. The "compact" reformulation

2µ Ft D(u) : D(ũ) dx - Ft p∇ • ũ dx + n i=1 m i U i • Ũi + J i ω i • ωi = Ft f • ũ dx. ( 10
)
of equations ( 5)-( 7)- [START_REF] Maury | Direct simulations of 2d fluid-particle flows in biperiodic domains[END_REF] shows that boundary terms have been absorbed by the choice of the finite element spaces.

The constraints on relative velocities

Among the n bodies composing the articulated swimmer, we identify B n as the reference body. The velocities U i of all the other bodies B i , i = 1 . . . n -1, are expressed as functions of U n via constraints of the form

U i = U n + W in (t), i = 1 . . . n -1, (11)
where W in (t) represents the relative velocity between B i and B n . The addition of these constraints to (3) completes the formulation of the swimming problem for the articulated swimmers. We notice that the resulting system is a particular instance of the general case presented in [START_REF] Alouges | Self-propulsion of slender micro-swimmers by curvature control: N-link swimmers[END_REF], where u d (t) is given by combining the constraints in [START_REF] Wagner | The extended finite element method for rigid particles in stokes flow[END_REF].

The formulation we just described applies directly to the three-sphere swimmer [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF], an articulated swimmer composed of three aligned spheres. Here the reference body B n = B 3 is the central sphere, which is connected by extensible arms to the other two spheres B 1 and B 2 . The formulation can be applied as well to the planar three-sphere swimmer [START_REF] Lefebvre-Lepot | A stokesian submarine[END_REF] or to the four-sphere swimmer [START_REF] Alouges | Optimally swimming stokesian robots[END_REF], whose spherical bodies are placed on the vertices of an equilateral triangle or tetrahedron, respectively. In those cases, the relative velocity vectors W in should be carefully computed, as each extensible arm connects B i to the barycenter of the swimmer, and not to B n directly as in the case of the three-sphere swimmer.

Discrete and algebraic formulation

Let us consider a geometrical discretization F h of the fluid domain using simplices. The fluid problem is then discretized using an inf-sup stable pair of conforming finite element spaces

X h (F h ) × B h (F h ) for (u h , B h ) and (U j , ω j ) ∈ [R d ] n × [R d ] n .
Since the fluid velocity satisfies u

= U i + ω i × (x -x CM i
) on ∂B i , the degrees of freedom u ∂Bi that lie on ∂B i are treated differently from the remaining ones, that we denote by u I . Indeed, u ∂Bi are expressed as a function of (U i , ω i ).

At first, equations ( 5)-( 8) are discretized by ignoring the constraint u = U i + ω i × (x -x CM i ) on ∂B i . Then, via the operator

P =      I 0 0 0 0 PUi Pωi 0 0 I 0 0 0 0 I 0 0 0 0 I      , (12)
that satifies the equation (u I , u ∂Bi , U i , ω i , p) T = P(u I , U i , ω i , p) T , the change of finite element basis, from the standard formulation to the constrained one, is performed. In [START_REF] Walker | Boundary behaviours of leishmania mexicana: A hydrodynamic simulation study[END_REF], PUi and Pωi are the interpolation operators that allow the expression of u ∂Bi as a function of U i and ω i .

The constraints on the relative velocities between the bodies can be imposed via Lagrange multipliers, through a modification of the operator P or by modifying the matrix that results from the discretization of the fluid problem.

The first method is the least invasive: additional equations and unknowns are added to a pre-existing discretized problem of n independent bodies in a Stokes fluid, as described in [START_REF] Chabannes | High-order fluid-structure interaction in 2d and 3d application to blood flow in arteries[END_REF]. Lagrange multipliers α i ∈ R d , i = 1 . . . n -1 are introduced to impose the constraints on the translational velocities U i onto the differential formulation. The previous constraint will appear in the equations that describe the rigid body motion of the solid bodies: equations ( 7) will be substituted by

m i Ui • Ũi + α i • Ũi = - ∂Bi (-pI + 2µD(u)) n • Ũi dS, i = 1 . . . n -1, (13)
m n Un • Ũn - n-1 i=1 α i • Ũn = - ∂Bn (-pI + 2µD(u)) n • Ũn dS, ( 14
)
α i • (U i -U n ) = α i • W in , i = 1 . . . n -1. ( 15
)
The addition of Lagrange multipliers entails the modification of P by providing an additional identity matrix of size d(n -1) × d(n -1) on the diagonal.

Instead of using Lagrange multipliers to constrain the translational velocities, a modification of the operator P could give the same results. In terms of finite element spaces, this consists in reducing the constrained finite element space to basis functions that satisfy

u = U n + ω n × (x -x CM n (t)) + u d (t), with u d (t) function of W in .
The dn × dn block that corresponds to the translational speeds, that presently is an identity matrix I dn , would be substituted by a block of the form

B 1 { B i { B n { I d . . . 0 d . . . -I d 0 d . . . I d . . . -I d -I d . . . -I d . . . I d . (16)
This modification must be applied directly to the construction of P, because an a posteriori change of the operator would be too costly in terms of operations on the compactly stored matrix. The relative velocities will then be imposed on the right-hand side of the problem, as the known part of the body velocity.

The last modification is even more invasive, and it consists in modifying the system matrix and right-hand side by substituting [START_REF] Lefebvre-Lepot | A stokesian submarine[END_REF] with

(m i U i -m n U n ) • Ũi = W in (t) • Ũi , i = 1 . . . n -1, m n Un • Ũn = - ∂Bn (-pI + 2µD(u)) n • Ũn dS. (17)
4 The three-sphere micro-swimmer

The three-sphere micro-swimmer [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF] is a three-dimensional swimmer composed of three aligned spheres having the same radius R. The two outer spheres are connected to the central one by extensible links, and the propulsion of the swimmer is ensured by changing the lengths of the connecting arms between two fixed values. The arm shrinkage is performed in a non-reversible fashion, in order to break the time-reversal symmetry of the Stokes equations, with a constant relative speed between the central sphere and the approaching one. For example, if the left arm is shrinking and the right arm keeps its length fixed, the translational speed U 3 and U 2 of the central and right sphere coincide, while the speed of the left sphere, moving with relative speed W 13 with respect to the center sphere, will be U 3 + W 13 . The (non-reciprocal) stroke is composed of 4 steps, as shown in Figure 1, where the lengths of the two arms are alternatively modified.

A quantitative example of swimming is presented in [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF]: if D = 10R is the length of each link in its rest position and ε = 4R is the maximal variation for the length of each arm (leaving a link length of 6R when the arms are completely shrunk), the center sphere is displaced by 0.16R in the positive direction at the end of the 4-step stroke. In the first step, the travelled distance is 1.35R in the negative direction; in the second step, it is 1.44R in the positive direction; in the third step, it is 1.44R in the positive direction; in the fourth step, it is 1.35R in the negative direction.

L L B1 B2 B3 L-a L L-a L-a L L-a L L
Using the aforementioned Lagrange multipliers formulation, we are able to recover the displacement at each step of the 4-step stroke reported in [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF]. Figure 2, on the left, translates the steps of the body deformation in terms of relative velocities between the central and lateral spheres. Figure 2, on the right, represents the motion of B 3 during several repetitions of the 4-step stroke. The results were obtained using the library FEEL++ [START_REF] Chabannes | High-order fluid-structure interaction in 2d and 3d application to blood flow in arteries[END_REF] and in particular its toolbox for Navier-Stokes in moving domains including moving rigid bodies. The implementations of Lagrange multiplier and P modification formulations are available in FEEL++ Github repository [START_REF] Vincent Chabannes | Github feel++ repository[END_REF] and can be used to reproduce the results in sequential and parallel.

Conclusion

In this paper, we provide a numerical method to simulate a self-propelled micro-swimmer composed of rigid bodies. Different formulations are proposed, including one based on Lagrange multipliers, to impose the relative motion between its components. The correctness of our formulations is verified on the three-sphere micro-swimmer by comparing the displacements obtained numerically to the ones in [START_REF] Najafi | Simple swimmer at low reynolds number: Three linked spheres[END_REF], even if only the results based on Lagrange multipliers are presented here. Current work includes the treatment of other multi-body swimmers like the planar 3-sphere swimmer [START_REF] Lefebvre-Lepot | A stokesian submarine[END_REF] or the 4-sphere swimmer [START_REF] Alouges | Optimally swimming stokesian robots[END_REF], formulating the ALE framework with mesh adaptation, the extension to other fluid models, e.g. Navier-Stokes and non-Newtonian, and the coupling with elasticity models to handle deformable bodies, i.e. swimmers.

Figure 1 :

 1 Figure 1: Representation of the three-sphere swimmer and its swimming gait. The gait is composed of four strokes in which one of the arms is alternatively shrunk or elongated. The alternation guarantees the non-reciprocity of the motion.

Figure 2 :

 2 Figure 2: The left figure presents the relative speed W 13 , between the central and left sphere, and the relative speed W 23 , between the central and right sphere, as functions of time. The right figure shows, in blue, the position of the central sphere during the 4-step stroke. The red line intersects the trajectory of the central sphere at the red circles, which mark the 0.08R and 0.16R displacements predicted in [9] after 2 and 4 steps composing the swimming stroke.

  and they satisfy the relationship

	∂Bi


Ibtissem Chouba
email: ibtissem.chouba@utt.fr

Lionel Amodeo
email: lionel.amodeo@utt.fr

Farouk Yalaoui
email: farouk.yalaoui@utt.fr

Taha Arbaoui
email: taha.arbaoui@utt.fr

David Laplanche
email: david.laplanche@ch-troyes.fr

-Morocco "New advances and challenges for sustainable and smart industries" A MIXED INTEGER LINEAR PROGRAM FOR HUMAN AND MATERIAL RESOURCES OPTIMIZATION IN EMERGENCY DEPARTMENT

Keywords: Emergency department, resources optimization, mixed integer linear programming

The discrepancy between patient demand and the emergency departments (ED) capacity, that mainly depends on human resources and on beds available for patients, often lead to ED's overcrowding and to the increase in waiting time. In this paper, we focus on the optimization of the human (medical and paramedical staff) and material resources (beds) in the ED of the hospital center of Troyes, France (CHT). We seek to minimize the total number of waiting patients from their arrival to their discharge. We propose a mixed integer linear program solved by a sample average approximation (SAA) approach. The program has been tested on a set of real data gathered from the ED information system. Numerical results show that the optimization of human and material resources leads to a decrease of total number of waiting patients.

INTRODUCTION

Over the last decades, ED overcrowding has known an important and serious growth around the world. This is due to several hospital's internal and external factors such as: high demand (patient overflow) for services, seasonal epidemics, long waiting time and limited resources. Thus, this phenomenon is known to have adverse impacts on the performance and quality of ED services. To overcome these problems, resources' optimization is an excellent approach by which an efficient system can be achieved [START_REF] Cabrera | Simulation optimization for healthcare emergency departments[END_REF].

In this context, human and material resource planning of the healthcare system is widely reported in the literature. Several works have introduced simulation, mathematical, linear programming or heuristics models for the desired optimization [START_REF] Cabrera | Simulation optimization for healthcare emergency departments[END_REF][START_REF] El-Rifai | A stochastic optimization model for shift scheduling in emergency departments[END_REF][START_REF] Daldoul | A stochastic model to minimize patient waiting time in an emergency department[END_REF][START_REF] Chouba | Optimisation du service d'urgences à l'aide d'un modèle de simulation efficace[END_REF][START_REF] Feng | Stochastic resource allocation in emergency departments with a multi-objective simulation optimization algorithm[END_REF][START_REF] Ghanes | A comprehensive simulation modeling of an emergency department: A case study for simulation optimization of staffing levels[END_REF][START_REF] Ahmed | Simulation optimization for an emergency department healthcare unit in kuwait[END_REF][START_REF] Mazier | Scheduling inpatient admission under high demand of emergency patients[END_REF][START_REF] Verweij | The sample average approximation method applied to stochastic routing problems: a computational study[END_REF].

In 2009, Ahmed and Alkhamis [START_REF] Ahmed | Simulation optimization for an emergency department healthcare unit in kuwait[END_REF] proposed a simulationbased optimization approach as a decision support system for the ED in a hospital in Kuwait. The authors aimed to find the optimal staff number that allows to minimize the patients' waiting time. The results provided a decrease of 40% in patients' waiting time. All the levels of the ED process were taken into account in this study. However, optimization model was limited to the staffing of the ED only as other resources were not considered.

In 2012 and 2015, Carbera [START_REF] Cabrera | Simulation optimization for healthcare emergency departments[END_REF] and Ghanes et al. [START_REF] Ghanes | A comprehensive simulation modeling of an emergency department: A case study for simulation optimization of staffing levels[END_REF] were interested to minimize the patient's inpatient stay using simulation models. Ghanes et al. [START_REF] Ghanes | A comprehensive simulation modeling of an emergency department: A case study for simulation optimization of staffing levels[END_REF] proposed a simulation-based optimization model to improve the ED of the French urban hospital Saint Camille taking into account the stochastic aspect of patient's arrival. The new configuration of the ED staff yields a choice of 10% increase in staffing budget resulting 33% reduction of the patient's length of stay (LOS). Also, Carbera [START_REF] Cabrera | Simulation optimization for healthcare emergency departments[END_REF] proposed an agent simulation approach coupled with an exhaustive search optimization to determine the optimal ED staff configuration (emergency physicians, nurses and administrative staff).

In 2017, Feng et al. [START_REF] Feng | Stochastic resource allocation in emergency departments with a multi-objective simulation optimization algorithm[END_REF] presented a new multi-objective simulation optimization algorithm that integrates nondominated sorting genetic algorithm II (NSGA II) and multi-objective computing budget allocation (MOCBA) in order to determine the optimal medical resources allocation solution that minimize the patient's LOS and the total medical wasted cost.

El-Rifai et al. [START_REF] El-Rifai | A stochastic optimization model for shift scheduling in emergency departments[END_REF] proposed a stochastic mixed-integer linear program (MILP) solved by a sample average approximation (SAA) approach to allocate the physicians and nurses of the ED of the hospital center of Lille in France while taking into account the stochastic nature of patient's arrivals and service times. The authors aim to minimize the total expected patients' waiting time. Then, to evaluate the resulting personnel schedules, a discreteevent simulation model is used. In this work, only three queues are considered: two assessment queues (patients before and after having auxiliary examinations) operated MOSIM'20 -November 12-14, 2020 -Agadir -Morocco by physicians and a third queue operated by nurses for treatment of patients.

In 2018, based on the method proposed by El-Rifai et al. [START_REF] El-Rifai | A stochastic optimization model for shift scheduling in emergency departments[END_REF], Daldoul et al. [START_REF] Daldoul | A stochastic model to minimize patient waiting time in an emergency department[END_REF] used a stochastic MILP solved by SAA to optimize human (physicians and nurses) and also materials (beds) resources organization in the ED at the La-Rabta Tunisian public hospital. They aim to minimize the patient's waiting time, under the uncertainty of the number of patients' arrivals and their service time. In contrast to El-Rifai et al. [START_REF] El-Rifai | A stochastic optimization model for shift scheduling in emergency departments[END_REF], the authors consider the entire process of an ED from the triage of patient to bed assignment. In the two cases study [START_REF] El-Rifai | A stochastic optimization model for shift scheduling in emergency departments[END_REF][START_REF] Daldoul | A stochastic model to minimize patient waiting time in an emergency department[END_REF], numerical results proved that the patient's waiting time can be decreased by optimizing the number of the human and material resources.

Through the literature review, we note that not all works consider all levels of the process or optimize all the medical and paramedical resources. Also, few studies have utilized the SAA approach to solve the stochastic problems in ED [START_REF] El-Rifai | A stochastic optimization model for shift scheduling in emergency departments[END_REF][START_REF] Daldoul | A stochastic model to minimize patient waiting time in an emergency department[END_REF][START_REF] Mazier | Scheduling inpatient admission under high demand of emergency patients[END_REF].

In this paper, we focus on the optimization of human (medical and paramedical) and material (beds) resources of the ED of the CHT, France. We aim to increase patients' satisfaction by reducing the number of waiting patients from the triage process to the discharge or the assignment to a bed in the ED. Our model captures all the complexity of the ED process that varies according the patient's state (more details in Section 2). In addition, the uncertainty related to patient's arrival and to the service time are taken into consideration. To achieve this, we develop a stochastic MILP solved by SAA which is one technique of Monte Carlo sampling methods [START_REF] Verweij | The sample average approximation method applied to stochastic routing problems: a computational study[END_REF]. The proposed model is an extension of the problem proposed by Daldoul et al. [START_REF] Daldoul | A stochastic model to minimize patient waiting time in an emergency department[END_REF].

The main contributions of our approach compared to existing work, particularly that of Daldoul et al. [START_REF] Daldoul | A stochastic model to minimize patient waiting time in an emergency department[END_REF], are twofold. First, we consider all medical and paramedical resource in the optimization model. Second, to our best knowledge, our work is the first to consider the patient's arrivals forecasts that are introduced as input data in the optimization model. The expected values of patients' arrivals are obtained through forecasting models developed in a previous work [20]. These models proved a very good performance (up to 91,24 % for the annual total flow forecast) and a robustness during epidemic periods [START_REF] Afilal | Forecasting the emergency department patients flow[END_REF]. Using these forecasting models, we adapt or relocate necessary resources to meet demand on the long or short term runs.

MODELLING

Each process in the ED can be modeled as a queue where jobs to be served are the patients and the hospital staff are the servers that handle these jobs [START_REF] El-Rifai | A stochastic optimization model for shift scheduling in emergency departments[END_REF]. Presented as a queuing system, Figure 1 illustrates the patient's process in the ED of CHT. Two main circuits exist: the short circuit and the long circuit. Within the short circuit, they are mostly autonomous patients, who do not need major paramedical care. Thus, we consider five queues: the first queue (q=1) is for the triage process operated by a triage nurse; the second (q=2) and the third (q=3) queues are for medical consultation in short and long circuit respectively. These queues are operated by the medical (physician, internal) and paramedical (nurse, caregiver) staff. The fourth queue (q=4) is dedicated to auxiliary exams and is operated by nurses. The fifth queue (q=5) is devoted to patients waiting to be assigned to a bed in ED.

In the triage process (St,1), the nurse identifies the patient's state and the pathway they go through: either the short (α1, St,1) or the long (α2, St,1) circuit. Once an appropriate box becomes free, the physician makes a first assessment (St,2 or St,3 for the short and the long circuit respectively) and decides if the patient needs auxiliary exams to confirm or to refine their diagnosis (β1St,2 and γ1St,3), to be assigned to a bed (β2St,2and γ2St,3) or to be discharged to go home (β3St,2 and γ3St,3). In the case where tests are requested, a second assessment is required in the short or the long circuit (λ1St,4 and λ2St,4 respectively).

MATHEMATICAL FORMULATION

A mathematical formulation based on the details discussed previously is presented in this section. The objective is to increase the patient's satisfaction through the optimization of the number of human and material resources and the minimization of the number of waiting patients. i,n,t,q L (): Random variable representing the service capacity of the staff, i.e. the number of patients that n resources of type i can serve during period t for each queue q, under scenario .

Decision variables

, , , 1, n resources of type i woking at period t. 0 else. i n t q if there are exactly x  =   ,, i t q p : The number of resources of type i at period t, q

∈ [1..Q], Q≠5
,, i t q s : The number of resources of type i starting at period

t, i ∈ [1..I], q ∈ [1..Q], Q≠5
,, i t q e : The number of resources of type i ending at period t, i ∈ [1..I], q ∈ [1..Q], Q≠5 t sb : The number of occupied beds at the start of period t t eb : The number of liberated beds at the start of period t t b : The number of beds at period t , tq w (): The number of waiting patients in queue q at the start of period t, under scenario .

, tq S (): The number of patients served in a queue q during period t, under scenario .

Resource optimization stochastic model:

, , 1 1 1,1 1, ( ) : ( ( ( )) (1)
( ) ( )

T Q t q t q t q q P Min E W S W A      = =   -   =   1,q 1,q (2) ( ) 0 q [2..Q] (3) ( ) 0 W S   =   =  t,1 t 1,1 t t 1,1 t,2 t 1,2 1 t 1,1 1 t 1,1 t 1,2 t,3 t 1,3 q [1..Q] (4) ( ) ( ) ( ) ( ) t [2..T] (5) ( ) ( ) ( ) ( ) ( ) t [2..T] (6) ( ) ( W W A S W W S S S W W             - - - - - - -  = + -   = + + -   = 2 t
) ( ) ( ) ( ) t [2..T] (7) ( ) ( ) ( ) ( ) ( ) t [2..T] (8) ( ) ( ) ( ) ( ) ( ) S S S W W S S S W W S S S                     - - - - - - - - - - - + + -   = + + -   = + + - i,1 t,q t,q t,1 i, , ,1 i, , ,1 1 t [2..T] (9) ( ) ( ) t [1..T] q [1..Q] (10) ( ) ( )* t [1..T] i [1..I] N n t n t n S W S L x     =              ,2 ,3 t,2 , , ,2 , , ,2 1 t,3 , , ,3 , , ,3 1 (11)
( ) ( )* t [1..T] i [1..I] (12) ( ) ( )* t [1..T] i [1..I] i i N i n t i n t n N i n t i n t n S L x S L x     = =             3,4 t,4 3, , , 4 3, , ,4 1 t,5 (13)
( ) ( )* t [1..T] (14)
( ) t [1..T] N n t n t n t S L x S b    =        (15) ,, 1 , , , , 1 i
[1..I] q [1..Q] q 5 (16) ( ) 0 i [1..I] q [1..Q] q 5 ( 17
) T i t q i t T i t q i t q t p TT s p = =       - =        , , , ', , ', ,1, ' 2 ,1, , , ,1, ,1, , ', ' ( ) t [2.
.T] i [1..I] q [1..Q] q 5 (18) i [1..I] q [1..Q] q 5 (19) t i t q i t q i t q i q t i q i t q i q i q i t q t t p s e p p p s e e

= = = - +        = + -        mod (t+LBD, T+1) , , ,1, 1 mod
(t+UBD+1, T+1) , ', , , ' 1 1 t [1..T] i [1..I] q [1..Q] q 5 (20) t [1..T] i [1..I] q [1..Q] q 5 (21) i t q i q i t q i t q t t p s e p b Nb + =+ -                =   1 ( 22
) eb NB = 1 1 (23) 0 ( 24
) t t t t sb b b eb sb - = = + - ,q ,q , , , , , 1 , , , 1 t [2..T] ( 25
) * t [1..T] i [1..I] q [1..Q] q 5 (26) 1 i i N i n t q i t q n N i n t q n n x p x = =   =        =   ,, , , ,t, t [1..T] i [1..I] q [1..Q] q 5 (27) 1 t [1..T] i [1..I] ( 28
)
1 i t q i t q i q p p s              + t [1..T] i [1..I] (29)    
The objective function (1) aim to minimize the expected value (E[.]) of the total number of waiting patients while respecting to the distribution of the random value of patient arrival and service times.  designates one scenario of patient arrival and service time in the ED. Constraint (2), (3), and (4) initialize the number of patients in the queues in each scenario. Constraint ( 5) to ( 9) update the number of patients in each queue at every period. Constraint [START_REF] Afilal | Forecasting the emergency department patients flow[END_REF] specifies a logical bound for the number of patients served during a period, that cannot be larger than the number of patients in the queue. Constraints (11) through (15) limit the number of patients MOSIM'20 -November 12-14, 2020 -Agadir -Morocco served in each queue at every period. Constraint (16) limits the total number of staffing hours available for each type of resource. Constraints (17), (18), and ( 19) are logical constraints that define the range of values that the variables p, s, and e can take. Constraint (20) ensures that at least si,t,q resources of type i will not end their shift from t to (t + LBD), for each queue q. Moreover, constraint (21) ensures that the last pi,t,q resources of type i end their shift from t + 1 to (t + UBD + 1). Constraints ( 22), (23), and (24) initialize the number of existing beds, liberated and occupied. Constraint (25) defines the number of available beds at every period. Equations ( 26) and ( 27) assign a value to the variable x, which is used to limit the number obtained in constraints (11) through (15). Finally, constraints (28) and ( 29) are operational constraints to ensure that there are different moments of the day where resources of type i start working and that there is at least one resource of each type in every period.

NUMERICAL RESULTS

In this section, we present the numerical results of the MILP solved by the SAA method. The main idea of the SAA method is to approximate the expected objective function of the stochastic problem by a sample average estimate derived from a random sample. The resulting sample average approximating problem is then solved by deterministic optimization techniques. The process is repeated with different samples to obtain candidate solutions along with statistical estimates of their optimality gaps [START_REF] Verweij | The sample average approximation method applied to stochastic routing problems: a computational study[END_REF].

To solve the model, we use the commercial solver IBM ILOG CPLEX 12.5. All the parameters, historical performances and input data are based on real data gathered from the ED's information system and covering the whole year 2019. Based on these data, we use the input analyzer of the ARENA software to generate the distributions of patient's arrivals and services' times in each queue for each period of the scheduling horizon. We assume that the scheduling time horizon is equal 6 hours and the time period (time step) is equal to 1 hour.

To evaluate the performance of the proposed approach, we varied K, the number of scenarios, to obtain K independent random samples of the patient's arrival (At 1 …At k ) and the services' times (Lt 1 …Lt k ). The different values of K are [START_REF] Feng | Stochastic resource allocation in emergency departments with a multi-objective simulation optimization algorithm[END_REF][START_REF] Afilal | Forecasting the emergency department patients flow[END_REF]20,30,40,50]. For each of these values, we run the model with seven different replications to produce seven Monte Carlo optimum solutions. The expected total patient's waiting number stabilizes for approximately K = 30 scenarios. We notice that the bigger K is, the better the solution tends towards convergence.

Optimization results are presented in Table 1. We consider 30 scenarios and we present the results of the proposed configurations. First, Table 1 provides the performance (total number of waiting patients) of the current configuration of the ED system which corresponds to the average of the different real scenarios. Second, the performance of the reproduction model (Configuration 1), in which the same number of resources available in the ED are considered, is displayed to evaluate the robustness of our model. Third, the performance of the optimization model (Configuration 2) is shown. Using this model, we measure the impact of the optimizing the number of medical and paramedical resources on each queue q. In our experiments, we consider that the maximum number of both medical and paramedical staff per period is equal to 3 (Physician P1=3, Internal P2=3, Nurse P3=3, Caregiver P4=3), and that the number of beds is 17.

Table 1 shows that the total number of waiting patients associated with Configuration 1 (104 patients) is close the real total number of waiting patients (116 patients). For Configuration 2, our model improves the performance by 24% by decreasing the total waiting patients from 104 patients to 79 patients. Optimization results demonstrate that reducing the total number of waiting patients requires adding resources in certain periods and queues.

Incorporating forecasts of patients' arrival

In our work, instead of real arrivals, we introduce patients' arrivals forecasts (Afilal et al. [START_REF] Afilal | Forecasting the emergency department patients flow[END_REF]) of the ED of CHT as input data for each period in the optimization model. Following the analysis of patients' arrival data, we noticed that the arrival rate follows a typical behavior whatever the day of the week and the hour of the day. As a result, we can calculate the estimated number of patient arrivals per hour (estimated arrival value* percentage of arrivals per hour). For the resolution, the expected arrival value will be set K times for each period and we generate K independent random samples of the services time.

We test the model, with data of 29-November 2019 (Figure 2) where the mean absolute error between arrival forecasts and observation equal to zero and the number of patients' arrivals was 10% higher than the average daily arrival. This allow as to test the sensitivity and the robustness of the model during an epidemic period for example.

Obtained results are summarized in Table 2. We compute the total number of waiting patients number and the number of required physicians.

Results demonstrate that the total number of waiting patients increases with the increase in the number of patients. Thus, by comparing the optimal configurations, we observe that the number of physicians remains stable in long circuit, except in the first two. However, for the short circuit, we have to add two physicians in average.

This paper is about the optimization of human and material resources in the ED of CHT. We proposed a generic approach to appropriately adjust the number of human and material resources from the patient arrival into the ED to the discharge or the assignment to emergency beds, taking into account uncertainties related to the patients' arrival and the services times. We proposed a stochastic MILP solved using a SAA method to determine the optimal allocation of medical and paramedical staff and beds, thus minimizing the total number of waiting patients in all the queues. The program was tested on a set of real data gathered from the ED's information system. The experimental results demonstrated that using the proposed approach could reduce the number of waiting patients. In contrast to most existing works, this paper considers optimizing all the humans and materials resources used in the overall patient process. Moreover, to the best of our knowledge, this is the first work that considers the use of forecasts of patients' arrivals as input data to help the decision maker to schedule daily the human and material resources.

In future work, we will consider testing this approach on other emergency structures in order to evaluate their performances and propose new resource organizations and to take into consideration other uncertainties such as the availability of medical staff, the availability of equipment, etc.

Figure 1 :

 1 Figure 1 : Patient flow in ED service

3. 1

 1 Parameters T: Planning horizon t: Different period of the day [1..T] I: The number of different types of human resources i: A type of human resource available [1..I] Q: The total number of queues q: The queue number [1..Q] Ni,q: The total number of each human resource i, i ∈ [1..I], q ∈ [1..Q], Q≠5 Nb: The total number of beds n: Number of resources of each type [1..Ni,q] αj: Coefficient associated with the patient flow served in the triage box, j ∈ [1..4] βj: Coefficient associated with the patient flow served in the box of consultation in short circuit, j ∈ [1..4] γj: Coefficient associated with the patient flow served in the box of consultation in long circuit, j ∈ [1..2] λj: Coefficient associated with the patient flow served in the box of auxiliary exams, j ∈ [1..2] TTi: The number of total working hours available to schedule resources of type i LBD: Lower bound of the shift length UBD: Upper bound of the shift length At (): Random variable representing the number of patients arriving at each period t, under scenario .

17 17 17 17 17 17 17 17 17 17 17 17 17 14 14 13 Performance 116 patients 104 patients

MOSIM'20 -November 12-14, 2020 -Agadir -Morocco Table 1 : Optimization results for human and material resources