Fabrice Darinot
email: fabrice.darinot@espaces-naturels.fr

Quentin Le Petitcorps
email: quentin.lepetitcorps@gmail.com

Véronique Arnal
email: veronique.arnal@cefe.cnrs.fr

Aurélie Coulon
email: aurelie.coulon@mnhn.fr

Le Quentin

Petitcorps

Arnal Ve´ronique

Coulon Claudine Aure´lie

Montgelard

C Montgelard
email: claudine.montgelard@cefe.cnrs.fr

Effects of landscape features and flooding on the genetic structure of a small wetland rodent, the harvest mouse (Micromys minutus)

Keywords:

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

Introduction

Wetlands are one of the ecosystems with the highest biodiversity in the world, but they are fragile and their species and habitats are among the most threatened. Indeed, the Millennium Ecosystem Assessment showed that environmental degradations are the most prominent within wetland ecosystems [START_REF] Finlayson | Ecosystems and human well-being: wetlands and water: synthesis[END_REF], with about half of the global wetland area lost as a result of human activities during the twentieth century [START_REF] Zedler | Wetland resources: status, trends, ecosystem services, and restorability[END_REF]. Water is the most important component of wetland ecosystems, controlling soil characteristics and associated plant and animal communities. Wetlands are inundated or saturated near the surface long enough to influence vegetation and fauna, resulting in species that are adapted to tolerate wet environments and that are not found in other habitats. Among them, many species of mammals depend, to diverse degrees, on wetland habitats for survival.

Dispersal plays an important role on global population survival in wetlands. Indeed, wetlands often occur as discrete patches within a matrix of upland habitat, where many wetland-dependent organisms like amphibians or invertebrates are organized as metapopulations [START_REF] Gibbs | Importance of small wetlands for the persistence of local populations of wetland-associated animals[END_REF][START_REF] Schooley | Metapopulation dynamics of wetland species[END_REF]. Their local populations are often small and thus vulnerable to extinctions caused by demographic events, diseases or natural events such as changing water levels, succession of wetland vegetation and habitat destruction by floods. In this context, dispersal among local populations is necessary to maintain the global metapopulation. The temporal fluctuations in water levels are an important factor which particularly affects small mammal assemblages [START_REF] Amori | The influ-_ ence of spring floods on small mammal communities in the Nemunas River Delta, Lithuania[END_REF]. Floods, as an extreme hydrological constraint, can even jeopardize population sustainability, depending on the level of dispersal that is possible between local populations.

However, dispersal has been poorly studied in wetland species, except for seed dispersal in hydrochorous plants [START_REF] Cunnings | Fluvial seed dispersal of riparian trees: transport and depositional processes[END_REF][START_REF] Fraaije | Spatial patterns of water-dispersed seed deposition along stream riparian gradients[END_REF]. For example, the Round-tailed muskrat (Neofiber alleni) in Florida, USA [START_REF] Schooley | Enhancing the area-isolation paradigm: habitat heterogeneity and metapopulation dynamics of a rare wetland mammal[END_REF] was shown to function as a metapopulation. But overall, little is known about the genetic structure of wetland mammal populations in these naturally patchy habitats. Especially, the influence of waterways on gene flow in small mammals is still poorly documented [START_REF] Mergey | Identifying environmental drivers of spatial genetic structure of the European pine marten (Martes martes)[END_REF], while wetlands often depend on rivers that maintain hygrophilous habitats and strongly structure the landscape. Moreover, only a few studies have addressed the effect of flooding on genetic structure of small mammal populations [START_REF] Matocq | Population genetic structure of two ecologically distinct Amazonian spiny rats: separating history and current ecology[END_REF][START_REF] Kittlein | Interactions between environmental factors can hide isolation by distance patterns: a case study of Ctenomys rionegrensis in Uruguay[END_REF][START_REF] Rocha | Seasonal flooding regime and ecological traits influence genetic structure of two small rodents[END_REF]. It is assumed that floods may have similar consequences on the patterns of genetic diversity as wildfires in heterogeneous landscapes, especially in promoting gene flow [START_REF] Pereoglou | Landscape genetics of an early successional specialist in a disturbance-prone environment[END_REF]. Indeed, because fires increase the extent of suitable habitat for some small mammals tied to open areas, isolated prefire populations expand and individuals disperse to colonize newly suitable habitat. Similarly, floods prevent the establishment of trees in wetlands and thus promote the maintenance of large areas of suitable open habitats, with increased connectivity between populations. Conversely, floods could generate massive deaths resulting in demographic bottlenecks and, as a result, reduced genetic diversity and gene flow levels.

In this study, we investigated population structure and gene flow in the harvest mouse (Micromys minutus), a threatened small mammal inhabiting wetlands. The harvest mouse is a tiny and elusive rodent, with adults weighting 7-8 g, widespread in the Palaearctic and Indomalayan regions, from Spain and Great Britain to Japanese 1 islands (Wilson and Reeder 2005). Whereas the IUCN has assigned the ''least concern'' status to this species [START_REF] Aplin | Micromys minutus[END_REF], the harvest mouse has the growing attention of scientists and naturalists in some countries like Switzerland and Great Britain where populations have declined in the past decade [START_REF] Perrow | What future for the Harvest mouse?[END_REF][START_REF] Blant | Nouvelles donne´es sur la re´partition de la souris des moissons (Micromys minutus Pallas, 1771) en Suisse occidentale et implications pour la gestion de son habitat[END_REF]. Recently, several programs of biological conservation survey and species management were planned in these countries in order to restore the populations of harvest mouse. Although its vernacular name refers to harvests, the original habitat of this species is linked to wetlands and the harvest mouse was present in Europe even before the cereal crops appeared in the Neolithic (Franck 1957). Hygrophilous and dense vegetations like reed beds, sedge and macrophyte habitats are now recognized as their main habitat [START_REF] Amori | The influ-_ ence of spring floods on small mammal communities in the Nemunas River Delta, Lithuania[END_REF]. A recent review using the diet of birds of prey confirmed that in Europe, the harvest mouse is most abundant in large wetlands [START_REF] Darinot | The Harvest mouse (Micromys minutus Pallas, 1771) as prey: a literature review[END_REF]. Furthermore, this rodent uses all three dimensions of its habitat to feed, reproduce and rest, often moving up the herbaceous compartment. As a result, this species can be expected to be less sensitive to flooding than strict terrestrial species. Finally, the harvest mouse is known to be a fast colonizer of new habitats [START_REF] Churchfield | Community structure and habitat use of small mammals in grasslands of different successional age[END_REF]), but its dispersal abilities are still unknown.

The aim of our study was to investigate the spatial scale and the landscape determinants of genetic structure in this wetland species. More specifically, the objectives of our study were to assess (i) the genetic structure of harvest mouse populations, at a regional spatial scale over several countries of Western Europe, and at a local scale in a flooded wetland in France; (ii) the effects of landscape features and floods on gene flow, at a local scale, in an area (Marais de Lavours, France) crossed by a large and embanked river (the Rhoˆne). A high level of genetic differentiation was expected between Western European populations of the harvest mouse, as it was shown in other small mammal species widely distributed in Europe, such as the common vole (Microtus arvalis) [START_REF] Heckel | Genetic structure and colonization processes in European populations of the Common vole, Microtus arvalis[END_REF] or the field vole (Microtus agrestis) [START_REF] Jaarola | Phylogeography of Field voles (Microtus agrestis) in Eurasia inferred from mitochondrial DNA sequences: phylogeography of Field voles[END_REF]. At the local scale, the study area is wide (45 km 2 ) and shows a low level of landscape fragmentation, for which a weak level of genetic structure was expected. Because the harvest mouse is known to be a good disperser, we hypothesized that it may be able to escape rising waters and then rapidly recolonize flooded areas. We hence hypothesized that floods have a weak effect on the population genetic structure of this species, similarly as in post-fire specialist mammals. The combination of both genetic structure and landscape genetic analyses will provide clues to inform conservation strategies of the species: the former will inform on the current genetic status of the different study populations (isolated or part of a larger group); the latter will enhance our understanding of which environmental factors most affect gene flow, and hence predict the potential effect of future land management actions on gene flow in the harvest mouse.

Material and methods

Study area, habitat description and sampling

The study was performed at two spatial scales: (1) a regional scale (Fig. 1a), including four sites distributed along 110 km of the Rhoˆne River in France (Prodon marsh near Geneva lake, Lavours marsh, Chautagne marsh, Saugey meander) and four sites 70-890 km away from the former (Potsdam pond in Northern Germany, Neuchaˆtel lake in Switzerland, Birieux pond, Meyranne marsh in Southern France); and (2) a local scale (Fig. 1b), focusing on Lavours marsh and Chautagne marsh, the latter being separated from the former by the Rhoˆne River. At each site, sampling was conducted in marshy habitats characterized by various plant communities, varying according to water table depth and flooding intensity (list of stations provided in Table 1, Online Resource 1).

Individual sampling was conducted between 2014 and 2016, except in Neuchaˆtel Lake (2012), using livetrapping. In each site, several stations were sampled with an average of 300 trap-nights (sampling details are provided in Table 2, Online Resource 1). Ear biopsies were taken in the field, sampling a piece of tissue of approximatively 3 mm 2 , without any injury to the animals (no bleeding in particular). Tissues were stored in 80% ethanol for DNA analysis. Tissues from Potsdam pond came from ten tail tips kindly given by the Tu¨bingen University (Deutschland), whereas tissues from the Neuchaˆtel Lake came from 13 livers (specimens numbers IZEA 7613,7617,7619,7620,7623,7624,7629,7748,7749,7571,7582,7584,7586) courtesy of the Cantonal Museum of Zoology of Lausanne (Switzerland).

Production of polymorphic microsatellite markers

We used six of the seven polymorphic microsatellite loci previously described for the harvest mouse [START_REF] Jiang | Polymorphic Microsatellite Markers for the Harvest Mouse (Micromys minutus) in Taiwan[END_REF] see Table 1 in Online Resource 2). To increase this number, a microsatellite library was specifically developed for M. minutus by the biotechnology company Genoscreen (Lille, France www. genoscreen.fr) using the method of highthroughput pyrosequencing (454 GS-FLX Titanium, Roche Diagnostics) as described in [START_REF] Malausa | High-throughput microsatellite isolation through 454 GS-FLX titanium pyrosequencing of enriched DNA libraries: pyrosequencing of SSRenriched DNA libraries[END_REF]. From Genoscreen we received potential primer pairs for 101 microsatellite loci (83-di, 14-tri, 3-tetra and 1penta-nucleotides), among which 40 markers (22-di, 14-tri, 3-tetra, 1-penta-nucleotides) were selected and tested for amplification and polymorphism. For each primer pair, melting molecular temperature (MT) was determined with OligoCalc [START_REF] Kibbe | OligoCalc: an online oligonucleotide properties calculator[END_REF], then tested for amplification with standard polymerase chain reaction (PCR) as follows: in a final volume of 10 ll containing 5 ll of Master Mixsolution from Qiagen TM (2X), 1 ll of each forward and reverse primers at 2 lM; 2 lL water and 1 lL of genomic DNA. All microsatellite loci were amplified with the following cycling: an initial denaturation step of 15 min at 95 C, followed by 30 cycles of denaturation (30 s at 94 C), optimal annealing (90 s at the optimal temperature; see Table 1, Online Resource 2) and elongation (1 min at 72 C), and ending with a final elongation step of 30 min at 60 C. DNA amplification was evaluated on agarose (2%) gel electrophoresis.

Genotyping

Genomic DNA was extracted from ear tissues using the Red-Extract-N-Amp kit Sigma-Aldrich. Multiplex PCRs were conducted in a 10 ll reaction volume containing 3 ll of H2O, 5 ll of the Master Mixsolution from Qiagen, 1 ll of each multiplexed primer at 2 lM and 1 ll of DNA. All microsatellite loci were amplified using the same cycling as above. Forward primers were end-labeled with fluorescent dyes FAM TM , VIC, NED TM , PET (Applied Biosystems). Amplified DNA was quantified by agarose (2%) gel electrophoresis. Then 3 ll of diluted PCR product were mixed with 0.2 ll size standard GeneScan TM 500 LIZ (Applied Biosystems) and 15 ll formamide (Hi-Di [START_REF] Peakall | genalex 6: genetic analysis in Excel. Population genetic software for teaching and research[END_REF]. Allelic richness was corrected for sample size with the rarefaction approach implemented with HP-RARE 1.0 [START_REF] Kalinowski | hp-rare 1.0: a computer program for performing rarefaction on measures of allelic richness[END_REF]. Hardy-Weinberg equilibrium was tested by estimating the differences between observed and expected heterozygosities using the program GenePop 4.6 [START_REF] Rousset | genepop'007: a complete re-implementation of the genepop software for Windows and Linux[END_REF], and p-values were estimated by the Markov chain method with 100,000 dememorizations, 1000 batches and 50,000 iterations per batch.

The level of genetic divergence among sites was estimated by calculating pairwise FST according to [START_REF] Weir | Estimating F -statistics for the analysis of population structure[END_REF], and their upper and lower confidence intervals by bootstrapping among loci (10,000 permutations), using the functions genet.dist and boot.ppfst in the R package hierfstat [START_REF] Goudet | HIERFSTAT, a package for R to compute and test hierarchical F-statistics[END_REF]). We tested for isolation-by-distance among the eight sites, using the maximumlikelihood population-effects (MLPE) approach [START_REF] Clarke | Confidence limits for regression relationships between distance matrices: estimating gene flow with distance[END_REF]: we ran a mixed model where the response variable was Fst/(1-Fst), the predictor variable was the logarithm of the geographic distance among sites and site identity was included as a random factor. Mixed models were run using the mlpe_rga function (package ResistanceGA, Peterman 2018) in R 3.6.2 (R Core Team 2019).

We also looked for the presence of genetic clusters with two approaches. First, we used STRUCTURE 2.3.4 [START_REF] Falush | Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies[END_REF]. Twenty independent runs were performed for each value of K (the number of genetic clusters) between 1 and 12 (to test for the possibility of existence of several genetic clusters within the study areas) using the admixture model and correlated allele frequencies, with prior information corresponding to the location of samples (Locprior option). This option is recommended when using a reduced number of markers or small sample sizes [START_REF] Porras-Hurtado | An overview of STRUCTURE: applications, parameter settings, and supporting software[END_REF]. For each run, Markov Chain Monte Carlo (MCMC) were carried out with 106 repetitions and a 5.10 5 burn-in period. The best number of genetic clusters was inferred through both the analysis of data likelihood for each value of K (LnL(K)) and Evanno's delta(K) method [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: a simulation study[END_REF], using STRUCTURE HARVESTER v0.6.8 [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF]. The results of genotype clustering from STRUCTURE were analyzed with CLUMPAK [START_REF] Kopelman | CLUMPAK: a program for identifying clustering modes and packaging population structure inferences across K[END_REF], which separates groups of runs representing distinct solutions and identifies an optimal cluster label alignment across different values of K, simplifying the comparison of clustering results across K. It was shown that when several hierarchical levels of population structure exist, the delta(K) method identifies the uppermost level [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: a simulation study[END_REF]. As a result, we conducted a hierarchical STRUCTURE analysis [START_REF] Coulon | Congruent population structure inferred from dispersal behavior and intensive genetic surveys of the threatened Florida Scrub-Jay (Aphelocoma coerulescens)[END_REF][START_REF] Janes | The K=2 conundrum[END_REF] to explore population subdivision. To do so, we repeated the same procedure described above (this time, with ten independent runs for each value of K between 1 and 5) on each of the genetic clusters obtained at the previous step, until no further structure was present. Because the results inferred by STRUCTURE can be altered in the presence of isolation-by-distance [START_REF] Meirmans | The trouble with isolation by distance[END_REF][START_REF] Perez | Assessing population structure in the face of isolation by distance: are we neglectingthe problem?[END_REF], we also used a discriminant analysis of principal components (DAPC), a multivariate method which allows to identify and describe clusters of genetically similar individuals [START_REF] Jombart | Discriminant analysis of principal components: a new method for the analysis of genetically structured populations[END_REF]. The DAPC was performed with the adegenet 2.0.0 R package [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF]). The number of groups K was first determined with the find.clusters function (keeping 60 principal components, based on the plot of cumulated variance). We then used a cross validation procedure (with the function xvalDapc) to assess the most appropriate number of principal components to use for the DAPC in itself. Finally, the DAPC was performed on the groups identified at the first step, based on 30 principal components and 5 discriminant functions.

Landscape genetic analyses at the local spatial scale

We first described the fine scale genetic structure within Lavours marsh and Chautagne marsh by testing for the presence of isolation by distance at this spatial scale (with the same MLPE approach as for the IBD test performed at the global scale). We also performed a spatial autocorrelation analysis, with Genalex [START_REF] Peakall | genalex 6: genetic analysis in Excel. Population genetic software for teaching and research[END_REF], based on distance classes of 500 m.

Then, we tested whether the Rhoˆne River acted as a barrier in the harvest mouse, by comparing the genetic distance of stations separated or not by the Rhoˆne, controlling for the geographic distance among them. The MLPE approach was performed on three matrices: a genetic distance matrix (Fst/(1-Fst)); a binary ''barrier'' matrix in which pairs of stations on the same side of the Rhoˆne took a value of 0, and pairs of stations on the opposite sides took a value of 1; and a geographic distance matrix (log(geographic distance among stations)). Because geographic distance and the « barrier» state are correlated, we first ran an MLPE model of genetic distance on geographic distance. We then performed a MLPE of the residuals of this first model on the barrier variable-this two-step approach allows us to assess the effect of the Rhoˆne River once the effects of the geographic distance are taken into account. This analysis was performed on a subset of 1 stations (listed in Table 1 in Online Resource 1), so that the distribution of geographic distances between the stations separated or not by the Rhoˆne was balanced (i.e. we only kept the pairs of stations that were separated by a distance ranging between 1.8 and 2.3 km).

Finally, in order to assess more finely the effects of land cover and flooding on gene flow, we took advantage of the more intense sampling in the Lavours marsh. The resistance of the different types of land cover on one hand and of the different degrees of flood intensity on the other hand, was assessed using the optimization approach implemented in ResistanceGA [START_REF] Peterman | ResistanceGA: an R package for the optimization of resistance surfaces using genetic algorithms[END_REF]. ResistanceGA is an R package which relies on a genetic algorithm to optimize resistance surfaces (i.e. GIS layers in which each cell is assigned a value representing the degree to which the corresponding land cover impedes movements-gene flow, here), looking for the set of resistance values that maximizes the correlation between genetic and resistance distances. This approach allows estimating without a priori the effect on gene flow of each type of landscape element represented in the resistance layer. Genetic distances were estimated with FST and resistance distances were calculated according to circuit theory [START_REF] Mcrae | Circuit theory predicts gene flow in plant and animal populations[END_REF][START_REF] Mcrae | Using circuit theory to model connectivity in ecology, evolution, and conservation[END_REF]), using the function Commute Distance. A linear mixed effect model with a maximum likelihood population effect parameterization was fitted to the data. In order to assess which environmental characteristics most affect gene flow in this species, we which included a large number of phytosociological units that were synthesized into nine habitat classes, taking into account the ecology of the harvest mouse: open water, large maize field, small stream with vegetation, mixed farming, wet meadows, reeds, poplar grove, forest, artificial soil (Fig. 2). We used vegetation maps at two time points to account for the potential delay between changes in land cover and their translation into the genetic structure of a population. From 1966 to 2012, the main land changes consisted in landscape fragmentation, with a reduction of wet meadows, an increase of woodlands and the emergence of largescale corn crops. We quantified those changes in habitat area and configuration between 1966 and 2012 with the Landscapemetrics R package [START_REF] Hesselbarth | Landscapemetrics: an open-source R tool to calculate landscape metrics[END_REF].

We also optimized resistance based on (iii) a map of winter flood intensity (Fig. 3), built from maps of the flood hazards of the Rhoˆne upstream on the administrative districts of Lavours marsh (French administration risk prevention units) and which consisted in four classes: unflooded artificial soil, water level below 1 m, between 1 and 1.5 m, and above 1.5 m. Finally, we used (iv) a map of emergent vegetation height, manually built (by expert knowledge) in combining vegetation and winter flood intensity maps. We identified five classes of vegetation: unflooded artificial soil, 1 completely flooded vegetation, emergent vegetation height below 1 m, between 1 and 1.5 m, above 1.5 m. The range of resistance values allowed in the iterations was 1-3500. We performed five runs for each model of resistance and we used AIC (Akaike's information criterion) to compare the performance of the different resistance surface transformations to explain observed pairwise genetic differentiation [START_REF] Peterman | Ecological resistance surfaces predict fine-scale genetic differentiation in a terrestrial woodland salamander[END_REF].

Results

New polymorphic microsatellite markers Among the 40 microsatellites tested, only 9 polymorphic loci were retrieved because 18 loci could not be amplified, 1 locus was unreadable with GeneMapper and 12 loci were monomorphic. In total, 15 polymorphic microsatellite markers were hence used in this survey (Table 1 in Online Resource 2): 9 new ones and 6 previously described by [START_REF] Jiang | Polymorphic Microsatellite Markers for the Harvest Mouse (Micromys minutus) in Taiwan[END_REF]. None of the 15 loci exhibited null alleles.

Genetic diversity

For the eight studied populations, genotypic data revealed a rather low level of diversity for all 15 loci: 76 alleles were detected in total with an average of 7.2 alleles per locus. Furthermore, the six loci from [START_REF] Jiang | Polymorphic Microsatellite Markers for the Harvest Mouse (Micromys minutus) in Taiwan[END_REF] presented twice the mean number of alleles per locus as compared to ours. Observed heterozygosity (Ho) at each locus ranged from 0.014 to 0.613, with an average value of 0.36 for all loci. Expected heterozygosity (He) at each locus ranged from 0.018 to 0.612, with an average value of 0.35 for all loci. However, microsatellites from [START_REF] Jiang | Polymorphic Microsatellite Markers for the Harvest Mouse (Micromys minutus) in Taiwan[END_REF] showed higher values than ours, with a mean Ho of 0.57 vs 0.22 and a mean He of 0.55 vs 0.22, respectively. A significant deficit in heterozygotes with respect to the Hardy-Weinberg equilibrium (HWE) was detected for the loci Mm10 and MME05 in the whole sampling (Table 1 in Online Resource 2). When applying the Bonferroni correction for multiple tests, no significant linkage disequilibrium was observed for each pair of loci. Mean allelic richness per population (AR) varied between 1.49 and 3.96, whereas seven populations exhibited private alleles (Table 1). Observed (Ho) and expected (He) heterozygosities per population ranged from 0.127 to 0.556 and from 0.159 to 0.560, respectively. Three populations significantly deviated from HWE. Genetic structure at the regional spatial scale Pairwise FST values were very low (FST\0.03) between Lavours marsh and Chautagne marsh. For all other sites, pairwise FST ranged from 0.14 to 0.71 (Table 2 in Online Resource 2) and all FST values were significant. There was a significant isolation by distance pattern at this regional spatial scale (slope = 0.1 ± 0.007).

The first analysis performed by STRUCTURE, on the whole dataset (309 individuals), showed that LnL(K) reached a plateau at K = 10 groups whereas the delta(K) approach indicated K = 2 as the most likely number of clusters (Fig. 4a). For K = 2, the 20 repetitions gave the same clustering with Lavours marsh and Chautagne marsh belonging to one group and all the other populations to the other one, with Saugey meander in an intermediate position. When looking for potential substructure within each of these two clusters, the analysis conducted on Birieux, Prodon, Neuchatel, Meyranne and Postdam (80 individuals) revealed that each of them represented a distinct cluster (Fig. 1a in Online Resource 2). The substructure analysis performed on Lavours, Chautagne and Saugey populations only (229 individuals) confirmed that Saugey meander clustered apart from the others (Fig. 1b in Online Resource 2). Finally, the STRUCTURE analysis performed on Lavours marsh and Chautagne marsh (211 individuals) revealed a very subtle level of genetic structuring, as shown by the large variability in LnL(K) values for K of 2 or more (Fig. 4b), with the station 12 of Lavours marsh that appeared genetically different from other stations of Lavours and Chautagne marshes. Basically, the 123 same results were obtained from K = 10 which identified in majority (14 out of 20 repetitions) Postdam, Meyranne, Neuchatel, Prodon, Birieux, Saugey and station 12 of Lavours marsh as separate clusters, whereas all stations of Lavours (12 excepted) and Chautagne marshes were not differentiated.

Results from the DAPC mainly confirmed the STRUCTURE analyses. On the whole dataset (309 individuals), the best number of clusters was K = 11 (Figs. 2 and3

Landscape genetic analyses at the local spatial scale

There was a weak but significant isolation-by-distance pattern at the local spatial scale (slope = 0.02 ± 0.016).

There was significant spatial autocorrelation in interindividual genetic distances for the first distance class only (\500 m) (Fig. 4 in Online Resource 2). When controlling for geographic distance, there was no effect of the Rhoˆne River on gene flow (« barrier» estimate: 0.002 ± 0.013). In the Lavours marsh, there was a significant effect of landscape elements on dispersal movements of harvest mice as indicated by the four resistance surfaces ''habitat 1966'', ''habitat 2012'', ''flood intensity'' and ''emergent vegetation height'' that explained more variation in the pairwise FST than Euclidean distance alone (Fig. 5, and Table 1, Fig. 1, Fig. 2 in Online Resource 3). The best model was the ''emergent vegetation height'', even if the difference was small relative to ''habitat 2012 and ''habitat 1966''. The model ''flood intensity'' alone was less supported than the other three models. In the ''emergent vegetation height'' model, resistance was the highest for the emergent vegetation class between 1 m and 1.5 m (mean resistance value = 377 ± 60), which corresponds to reed beds (Fig. 6d). Caution should be taken when interpreting STRUCTURE outputs performed in the presence of isolationby-distance [START_REF] Meirmans | The trouble with isolation by distance[END_REF][START_REF] Perez | Assessing population structure in the face of isolation by distance: are we neglectingthe problem?[END_REF]. However, results of the DAPC and STRUCTURE analyses were consistent, which suggests these genetic structure results are robust. An interesting finding was that allelic richness and heterozygosity indices, which reflect the genetic diversity of populations, were the highest in populations located in very large wetlands that are associated with numerous suitable habitats for the harvest mouse: Birieux pond (AR = 3.96, Ho = 0.556), Saugey meander (AR = 3.02, Ho = 0.433) and Lavours marsh (AR = 2.66, Ho = 0.353). Conversely, allelic richness and heterozygosity indices were the lowest in populations located in small or less favorable areas, like Potsdam pond (AR = 2.5, Ho = 0.400), Prodon marsh (AR = 1.74, Ho = 0.164) and Meyranne marsh (AR = 1.49, Ho = 0.127). At a finer scale, some genetic proximity seems to exist between both Lavours and Chautagne populations and Saugey meander, which are only 37 km away from each other, following the Rhoˆne stream. As the surrounding landscape is hilly and most probably unfavorable to the harvest mouse, it is likely that the Rhoˆne functions, or has functioned, as a corridor for this species (see paragraph below), by promoting gene flow through populations located in suitable vegetation on both banks like reed beds, small wetlands associated with the river, and small channels of vegetation.

At a local scale, our major finding is that Lavours (stations 1 to 11) and Chautagne marshes (stations 14 to 17) as well as station 13 of Lavours marsh constitute one single population in the studied landscape, an area over 45 km 2 . We found no signal of isolation by distance between extreme stations (i.e. between stations 1-16, 3-15, 11-16) in this area, and no barrier effect of the embanked Rhoˆne between them. Finally, only the station 12 from Lavours marsh appeared to be differentiated from other localities. All stations from Lavours but station 12 belong to the same watershed. Station 12 is isolated from the Lavours watershed by large forest areas, which may explain why the harvest mice in station 12 belong to a distinct population, a result that could also explain the HWE deviation observed for the Lavours marsh (including stations 1-13; see Table 1 in Online Resource 1). Several studies revealed that large or embanked rivers can represent partial barriers to dispersal for small or medium mammals, such as the Rhine River in Germany for the wild cat (Felis sylvestris) [START_REF] Hartmann | Potential barriers to gene flow in the endangered European wildcat (Felis silvestris)[END_REF], large rivers of Mexico for the painted spiny pocket mouse (Liomys pictus) (Garrido-Gardun˜o et al. 2016), the North Oconee River for the white-footed mouse (Peromyscus leucopus) in the United States [START_REF] Klee | Riverine peninsulas: an experimental approach to homing in White-footed Mice (Peromyscus leucopus)[END_REF] or large streams of South Africa for the natal multimammate mouse (Mastomys natalensis) [START_REF] Russo | Landscape determinants of fine-scale genetic structure of a small rodent in a heterogeneous landscape (Hluhluwe-iMfolozi Park, South Africa)[END_REF]. But in 1 some other species, it was shown that even wide rivers do not reduce gene flow, such as the Jurua River, a wide affluent of the Amazon River, for the arboreal spiny rat (Mesomys spp.) [START_REF] Patton | Gene genealogy and differentiation among arboreal spiny rats (Rodentia: Echimyidae) of the Amazon basin: a test of the riverine barrier hypothesis[END_REF], mountain rivers in the Alps for the common shrew (Sorex araneus) [START_REF] Lugon-Moulin | Do riverine barriers, history or introgression shape the genetic structuring of a Common shrew (Sorex araneus) population?[END_REF], or the 60 m wide Kolba¨ck River in central Sweden [START_REF] Wyttenbach | Genetic structuring and gene flow in a hybrid zone between two chromosome races of the Common shrew (Sorex araneus, Insectivora) revealed by microsatellites[END_REF] for the common shrew. As our analyses did not indicate any barrier effect at the scale of our study area, we can conclude that harvest mice are probably able to cross the Rhoˆne and this was most likely the case in the past. Before it was first embanked, 240 years ago, banks had gentle slopes with dense herbaceous cover, the stream was wider than today, with a lot of small islands covered by willow thickets and herbaceous vegetation, which might have constituted suitable halts for the mice. Moreover, meanders and braided channels were numerous and probably easy to cross for such good swimmers as harvest mice. Nowadays, the Rhoˆne is embanked and the width of the stream is between 100 and 350 m, with some rare islands, which lead us to hypothesize that this river now acts as a barrier to gene flow although this remains to be demonstrated. Indeed, a time lag effect does not seem likely, given the large number of generations since the embankment happened which can be estimated at about 700 generations [START_REF] Darinot | Dispersion et structure ge´ne´tique d'une population de Rat des moissons (Micromys minutus PALLAS, 1771) soumise a` des inondations re´gulie[END_REF]. Moreover, we can hypothesize that the effective population size on both sides of the Rhoˆne could be large enough to prevent the genetic signature of a barrier effect to develop. A recent study on an endangered butterfly, the False ringlet (Coenonympha oedippus), that was conducted in exactly the same area as our survey, also showed that in the past, and perhaps also today, the Rhoˆne was not a barrier to gene flow and that the population of this butterfly probably spread over the Lavours and Chautagne marshes as a whole [START_REF] Despre´s | Inferring the biogeography and demographic history of an endangered butterfly in Europe from multilocus markers[END_REF]. Although, such a genetic homogeneity is not often observed in small mammals, another example is provided by the common vole (Microtus arvalis), which is not a wetland species and which was documented to form one population over 500 km 2 [START_REF] Gauffre | Spatial genetic structure of a small rodent in a heterogeneous landscape[END_REF]. Although no such study exists for wetland small mammals, we can however mention a study performed on the water vole (Arvicola sapidus), which conversely showed strong spatial genetic structuring, a result consistent with spatially restricted dispersal [START_REF] Berthier | Genetic structure of the cyclic fossorial Water vole (Arvicola terrestris): landscape and demographic influences[END_REF]. In other species of Muridae with good dispersal capabilities, such as the wood mouse (Apodemus sylvestris), a substructure possibly exists even at a microgeographic scale (\3 km) but in this case it was attributed to restricted gene flow due to habitat fragmentation [START_REF] Berckmoes | Effects of environmental pollution on microsatellite DNA diversity in Wood mouse (Apodemus sylvaticus) populations[END_REF][START_REF] Booth | Spatial genetic structuring in a vagile species, the European wood mouse[END_REF]). The most suitable habitats for the harvest mouse are wetland habitats, but it is not a strict wetland species as it also lives in a wide variety of interstitial and wetland connected habitats [START_REF] Harris | Breeding season, litter size and nestling mortality of the Harvest mouse, Micromys minutus (Rodentia: Muridae), in Britain[END_REF]). Thus, the lack of population genetic structure at a local scale is most likely linked to this ecological valence, allowing this rodent to use a range of habitats which make up the wetland habitat matrix.

Effects of landscape features on gene flow

Surprisingly, results from ResistanceGA revealed little differences in resistance surfaces between 1966 and 2012. In 2012, the landscape was more fragmented than in 1966 (Table 1 in Online Resource 4) with a number of landscape units growing from 259 to 458, and an average area of these units decreasing from 12.7 to 7 ha. This fragmentation was probably not strong enough to subdivide the population, probably because of high dispersal abilities (demonstrated by their low level of local genetic structuring) and rather high ecological valence of the harvest mouse. Some habitats, such as large maize fields and poplar groves, which did not exist in 1966, are now recognized as one of the causes of the destruction of wetlands in Europe and their effect on harvest mouse dispersal had to be assessed. Maize fields were assigned very low resistance to gene flow despite the fact that they constitute a poor habitat for the harvest mice [START_REF] Darinot | The Harvest mouse (Micromys minutus Pallas, 1771) as prey: a literature review[END_REF]. This result might be explained by the fact that these crops are crossed by numerous small ditches covered with very suitable vegetation (tall grasses and reeds) for the species. Large maize fields are therefore a composite habitat, representing a strong support for dispersal and thus permeable to the species movements. Poplar groves form a special type of artificial afforestation, with a very developed and dense herbaceous stratum formed by the remaining meadows that existed before tree planting. This herbaceous stratum is still somehow attractive to harvest mice, even if this habitat is degraded due to the abundance of nitrophilous plants. This could explain why ResistanceGA assigned a lower resistance to gene flow to poplar groves than to natural forests that do not have any suitable vegetation for the harvest mouse. The ''streams with vegetation'' class that corresponds to the Se´ran River (Fig. 1b) and several large ditches with tall vegetation and reed on the banks, showed resistance values in 1966 and 2012 that did not significantly differ from the resistance of reeds (Fig. 1A and 1B in Online Resource 3). These two habitat types (ditches and reed) hence similarly negatively affect harvest mouse gene flow: despite the fact that they represent suitable habitats, they were assigned high resistance values, most likely due to their role of refuge during winter floods, as explained in the next section. Indeed, even if they are in large parts different, the resistance models based on habitat on one hand and on flooding on the other hand are partly redundant, the height of emergent vegetation being dependent on the type of vegetation (i.e. land cover).

The ''forest'' class resistance was rather low in 1966, which can be explained by the fact that the woods then formed linear habitats along the streams, i.e. riparian forests, which could be easily crossed by the mice. The ''forest'' class resistance was twice as high in 2012 as in 1966 while at the same time the forest area has increased from 794 to 1154 ha, making large wooded areas in the center of the marsh: unlike riparian forests, these wooden areas are very little permeable to mice movements and therefore to gene flow.

Effects of flooding on gene flow

Interesting results on the effect of winter flooding on gene flow at the landscape scale were provided by the outcomes of ResistanceGA. These floods occur between November and March, several times each winter, with varying degrees of intensity. Among the four models tested, the best one was the ''emergent vegetation height'' model, which indicated that winter floods have an effect on dispersal, with stronger resistance in areas where the height of vegetation that emerges from water is the highest. These areas correspond to the highest and tallest reed beds, which, as discussed above, also constitute the most resistant habitat surfaces in 1966 and 2012 (Fig. 1 in Online Resource 3). Moreover, the study of the distribution of private alleles within the Lavours marsh population revealed that they are only present in the most heavily flooded stations, which is consistent with low gene flow levels (locus Mm3 in stations 1, 3, 6 and locus Mm10 in station 4, see Fig. 1b). These observations (moving upwards and presence of private alleles in the most heavily flooded stations) strengthen the hypothesis that mice do not disperse during periods of floods. When floods arise, small mammal behavior is known to depend on their mobility, attachment to their habitat and ability to find refuges [START_REF] Zhang | Recovery of a rodent community in an agro-ecosystem after flooding[END_REF]). Species that are present in the floodplain can have two kinds of response: some move toward higher ground in order to escape floods, whereas others stay in the flooded area with an amphibious way of life and seek refuges in above water vegetation [START_REF] Sheppe | The annual cycle of small mammal populations along the Chobe River, Botswana[END_REF]. We demonstrated in complementary experiments based on live-trapping and camera-traps that some harvest mice respond in the latter manner, remaining in the reed during the flood, with no apparent escape movements nor effort to disperse (unpublished results). The same behavior was observed in the Salt marsh harvest mice (Reithrodontomys raviventris halicoetes), in the San Francisco Bay Estuary (U.S.A.), which remain in vegetation above the water during high tides [START_REF] Smith | Effects of natural and anthropogenic change on habitat use and movement of endangered Salt Marsh Harvest Mice[END_REF].

Another outcome of ResistanceGA is that submerged habitats or those with less than 1 m of emergent vegetation are the least resistant to gene flow, meaning that more dispersal occurred. These habitats mainly correspond to hygrophilous meadows that are suitable for the harvest mouse during the breeding season, but during floods, they can probably no longer persist because of the lack of emergent vegetation and animals are forced to escape, thus promoting dispersal.

Contrary to our initial hypothesis, our results suggest that the effect of flooding on gene flow in harvest mouse is quite different from that of fire. Indeed, floods are able to prevent wet meadows and reed beds from being closed by the colonization of trees, and thus they maintain large areas of suitable 1 habitat with potentially high connectivity. Unlike a fire which will completely remove the vegetation in place, a flood affects very little its characteristics, and in particular its cover, density and height. Thus, the attributes of the vegetation, which correspond to different habitats for the harvest mouse, condition the mouse response to flood, which is absolutely not true for a fire which kills all mammals that remain in the vegetation, regardless of its characteristics. This is shown by the higher resistance assigned to reed beds than to wet meadows. Although floods can result in some mortality of mice due to drowning, hypothermia or increased predation by birds of prey [START_REF] Wijnhoven | Flooding ecology of voles, mice and shrews: the importance of geomorphological and vegetational heterogeneity in river floodplains[END_REF], it is not as systematic as fires.

It would be interesting to identify refuge areas at the periphery of the wet meadows, which are used by migrating mice during floods. Dispersal appears to be high in the wet meadows, implying that individuals fleeing the flood do not return breeding in their native habitat. This hypothesis should be verified with individual movement monitoring. Finally, additional studies in the reed bed (unpublished results) have shown that some mice moved over long distances (several hundred meters) in the summer, outside the flood period. These movements were not highlighted in the wet meadows. This raises the question of dispersal during the reproduction period of the harvest mouse: does this dispersal exist and how does land cover affect it?

Conclusions

In general, the hydrological functioning of wetlands and in particular the flooding events have largely unknown effects on the genetics of the populations of the animal species that develop there. Our work is one of the few devoted to assess the effects of landscape and floods on the genetics of a small wetland mammal. Our results showed that the genetic structure of a wetland small mammal population can be quite homogeneous, despite major disturbances such as floods. Concerning the apparent lack of effect of river embankment on the dispersal of a small mammal such as the harvest mouse, our results must be taken with great care and they deserve to be deepened. It is by no means certain that the artificialization of rivers has no impact on gene flow in small mammals. Our study indicates that artificial vegetation (maize fields and poplar groves) can be permeable to the harvest mouse movements when suitable habitat (small ditches or herbaceous stratum) still persist. This knowledge is not only interesting from scientific considerations, but should also be useful for the development of plans to safeguard wetland ecosystems.

Fig. 1

 1 Fig. 1 Geographic distribution of the eight sites sampled for harvest mice (a: regional scale; b: local scale; N: number of sampled

Fig. 2

 2 Fig. 2 Maps of the 1966 and 2012 vegetation in the Lavours marsh

Fig. 3

 3 Fig. 3 Map of winter flood intensity in the Lavours marsh (averaged on the period 1944-2014)

  in Online Resource 2) with the individuals from Meyranne, Neuchatel, Prodon and Birieux clearly representing separated genetic clusters (8, 7, 10 and 6, respectively in Fig. 2 in Online Resource 2). The individuals sampled in Lavours marsh, Chautagne marsh and Saugey meander exhibited shared membership probabilities, and the corresponding clusters (1, 3, 4, 5, 9, 11) overlapped on the Fig. 4 Results of the STRUCTURE analysis performed on a the whole sampling and b Lavours-Chautagne stations. In both cases, the top panel represents the natural logarithm of the likelihood (LnL(K)) and of delta(K) computed following the Evanno's method (2005) for K = 1 to 12 for (a) and 1 to 5 for (b). The bottom panel is the bar plot for K = 2 from CLUMPAK, representing the proportion of the genome of each individual (one bar = one individual) assigned to each cluster (one color = one cluster). The bar plot was based on twenty independent runs for (a), 10 for (b)

Fig. 5

 5 Fig. 5 Boxplots of mean delta AIC (over 5 repetitions) of the mixed models explaining genetic distances by resistance distances based on different representations of the landscape:

Fig. 6

 6 Fig. 6 Boxplots of resistance values per environmental variable for each model. a habitats in 1966, b habitats in 2012, c flood intensity, d emergent vegetation height. Note the different scales in the four y axes

Table 1

 1 Genetic diversity of Micromys minutus in eight study sites

		Population genetic				HWE test
			parameters
	Population	N	AR	AP	Ho	He	Fis	Df	Chi-sq	P
	Lavours marsh	180	2.66	4	0.348	0.380	0.112	28	65.60	0.000***
	Chautagne marsh 31	2.51	0	0.374	0.372	-0.027	22	23.45	0.377
	Saugey meander	18	3.02	0	0.433	0.436	0.020	24	15.36	0.910
	Birieux pond	37	3.96	8	0.556	0.560	0.102	28	43.80	0.029*
	Prodon marsh	9	1.74	0	0.164	0.159	-0.075	12	4.72	0.967
	Neuchaˆtel lake	13	2.33	1	0.294	0.347	0.155	26	28.00	0.358
	Meyranne marsh 11	1.49	0	0.127	0.161	0.187	12	19.98	0.067
	Potsdam pond	10	2.50	5	0.400	0.393	-0.031	22	218.83	0.012*

N sample size, AR allelic richness, AP private allele, H0 observed heterozygosity, He expected heterozygosity, Fis fixation index, HWE Hardy-Weinberg equilibrium with probability (P) noted * when significant (P\0.05), ** when very significant (P\0.01) and *** when highly significant (P\0.001) 1 DAPC. The individuals from Postdam (cluster 2) partly overlapped with the group Lavours-ChautagneSaugey and the cluster 6 of Birieux. Finally, the main difference between STRUCTURE and DAPC analyses concerned the station 12 of Lavours that appeared as differentiated from other Lavours and Chautagne stations only with STRUCTURE.

Acknowledgements The authors would like to thank for people who help them in obtaining tissues of M. minutus for DNA isolation, for the description of new polymorphic microsatellite markers and genotyping of populations: Gre´goire Massez and Sylvain Ceyte (Re´serve Naturelle Nationale des Marais du Vigueirat), Ge´rard Hommay (Institut National de la Recherche Agronomique), Franc¸ois Le´ger (Office National de la Chasse et de la Faune Sauvage), Andrea Schu¨ster and Katarina Foerster (University of Tu¨bingen), Jacques Gillie´ron and Hubert du Plessix (Prodon marsh), Maurice Benmergui (Office National de la Chasse et de la Faune Sauvage, Birieux pond), Olivier Glaizot (Cantonal Museum of Zoology of Lausanne, Neuchaˆtel lake), Emilie Wichroff and Re´mi Bogey (Syndicat du Haut-Rhoˆne, Saugey meander), Manuel Bouron (Conservatoire des Espaces Naturels de la Savoie, Chautagne marsh). This study was supported by the French Ministry of Environment and the Departmental Council of Ain.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Arnaud Liefooghe
email: arnaud.liefooghe@univ-lille.fr

Sébastien Verel
email: verel@univ-littoral.fr

Tinkle Chugh
email: t.chugh@exeter.ac.uk

Jonathan Fieldsend
email: j.e.fieldsend@exeter.ac.uk

Richard Allmendinger
email: richard.allmendinger@manchester.ac.uk

Kaisa Miettinen
email: kaisa.miettinen@jyu.fi

Feature-based benchmarking of distance-based multi/many-objective optimisation problems: A machine learning perspective

Keywords: Multi/many-objective distance problems, Feature-based performance prediction, Automated algorithm selection

come   L'archive ouverte pluridisciplinaire

Introduction

Given a collection of problems and algorithms, the algorithm selection problem [START_REF] Rice | The algorithm selection problem[END_REF] is concerned with identifying an algorithm that is most suitable, in terms of some performance criteria, for a given problem at hand. An additional component is the availability of features characterising a problem. One can first extract problem features to generate a feature space, and then act on it as opposed to on the problem space. Significant research has been carried out on algorithm selection during the past two decades. Amongst others, tackling more efficiently a variety of continuous and mixed discrete/continuous [START_REF] Kerschke | Automated algorithm selection: Survey and perspectives[END_REF][START_REF] Smith-Miles | Cross-disciplinary perspectives on meta-learning for algorithm selection[END_REF], combinatorial [START_REF] Smith-Miles | Measuring instance difficulty for combinatorial optimization problems[END_REF] and multi-objective (continuous) optimisation problems [START_REF] Liefooghe | Landscape features and automated algorithm selection for multi-objective interpolated continuous optimisation problems[END_REF], as well as supervised learning [START_REF] Muñoz | Instance spaces for machine learning classification[END_REF]. For algorithm selection, extracting problem features that drive algorithm performance is critical; doing this efficiently is the focus of fitness and exploratory landscape analysis [START_REF] Malan | A survey of advances in landscape analysis for optimisation[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF]. Furthermore, understanding which and how problem features drive algorithm performance is valuable information when designing artificial problems of different complexity for benchmarking and tuning of algorithms. This information can then be used, e.g., to develop problem generators that can create test problems that meet user-defined problem characteristics. Such generators exist, for example, for many-objective distance-based optimisation [START_REF] Fieldsend | A visualizable test problem generator for many-objective optimization[END_REF] and cluster analysis [START_REF] Shand | HAWKS: Evolving challenging benchmark sets for cluster analysis[END_REF].

Our focus is to advance the area of feature design and algorithm selection for multi-and many-objective optimisation problems. This is motivated by the prominence of problems in practice, combined with our limited understanding on suitable multi-objective problem features [START_REF] Kerschke | Automated algorithm selection: Survey and perspectives[END_REF][START_REF] Malan | A survey of advances in landscape analysis for optimisation[END_REF]. The most relevant work is given in [START_REF] Liefooghe | Landscape features and automated algorithm selection for multi-objective interpolated continuous optimisation problems[END_REF], where features from landscape analysis (adapted from [START_REF] Liefooghe | Landscapeaware performance prediction for evolutionary multiobjective optimization[END_REF]) are applied on a benchmark set of 1 200 randomly-generated bi-objective interpolated continuous optimisation problems [START_REF] Zȃvoianu | Comparative run-time performance of evolutionary algorithms on multi-objective interpolated continuous optimisation problems[END_REF]. The study concluded that combining a classification model with a range of landscape features used as predictors can deliver a similar accuracy to predicting algorithm performance based on parameters used to generate the problems. It also investigated the relative importance of features for performance prediction and algorithm selection.

In this paper, we follow a similar approach to investigate the predictive power of parameters (problem features) used to generate distance-based multi/manyobjective point problems (DBMOPPs) proposed in [START_REF] Fieldsend | A visualizable test problem generator for many-objective optimization[END_REF]. More precisely, we (i) generate 945 problems with different characteristics (as defined by 7 problem features), then (ii) test the correlation between the problem features and the performance of three popular multi-objective evolutionary algorithms and one baseline approach (random search), and finally assess the problem features as predictors for (iii) algorithm performance prediction and (iv) algorithm selection on machine learning regression and classification, respectively. This is the first study of DBMOPPs and the generator/problem features proposed in [START_REF] Fieldsend | A visualizable test problem generator for many-objective optimization[END_REF].

The paper is organised as follows. Section 2 introduces DBMOPPs, together with the generator and its parameters (problem features) used to control the generation of such problems. Section 3 gives the experimental setup and discusses algorithm performance. Section 4 presents an experimental analyses of the problem features for automated performance prediction and algorithm selection. Finally, Section 5 concludes the paper and discusses further research.

Distance-based Multi/Many-objective Problems

In multi-objective optimisation (with 2 or 3 competing objectives) and manyobjective optimisation (with 4+ objectives), a set of solutions is typically sought that approximates the optimal trade-off combinations between the objectives, given any constraints on decision vectors. Formally, the tuple of a decision vector x, and its evaluation under the objective vector f () define a solution s = (x, f (x)). A solution s is said to dominate another s if s performs better than s on at least one objective, and no worse on all others. The maximal set of non-dominated decision vectors is known as the Pareto set, and its image in the objective space is known as the Pareto front.

A wide range of scalable test problem frameworks have been developed for multi-and many-objective optimisation, which are used (together with set quality indicators) to assess the performance of optimisers. These encapsulate a range of known problem characteristics (e.g. [START_REF] Brockhoff | Using well-understood singleobjective functions in multiobjective black-box optimization test suites[END_REF][START_REF] Deb | Scalable test problems for evolutionary multiobjective optimization[END_REF][START_REF] Huband | A review of multiobjective test problems and a scalable test problem toolkit[END_REF]). In addition, means for generating instances of problems have been created to prevent "tuning" of optimisers to particular suites of tests, to the detriment of performance on practical problems.

Frameworks for DBMOPPs have been developed over the last decade. They incorporate the range of features exhibited in other test suites (constraints, neutrality, multi-modality, dominance resistance regions, local fronts, etc.) and enable direct visualising of the search space in a plane. Initial work in this area includes [START_REF] Köppen | Fuzzy-Pareto-dominance and its application in evolutionary multi-objective optimization[END_REF][START_REF] Köppen | Substitute distance assignments in NSGA-II for handling many-objective optimization problems[END_REF], and [START_REF] Fieldsend | A visualizable test problem generator for many-objective optimization[END_REF] includes a summary of the features incorporated into this test problem design approach over the last 15 years. Arbitrarily many objectives can be defined. If |x| > 2 the decision vector is projected into two dimensions via a pair of orthogonal vectors prior to function evaluation. We now describe the construction of DBMOPPs, and the generator we recently developed. We will then extensively investigate its characteristics.

Properties and Features of DBMOPPs. Point-based distance problems are parameterised by sets of attractor vectors, where the minimum distance to a member of the ith set, V i , defines the ith objective value:

f i (x) = min v∈Vi dist(v, x). (1)
Further complexity can be added by imposing regions of constraint violation which locally adjust the distance function, and thereby can introduce discontinuities and neutrality to particular objectives, amongst other modifications.

In [START_REF] Fieldsend | A visualizable test problem generator for many-objective optimization[END_REF], we introduced a DBMOPP 6 instance generator, where problems can incorporate a range of properties, a subset of which are listed in Table 1 which we consider here. An example 3-objective problem with local fronts and dominance resistance regions is illustrated in Figure 1, along with its local dominance landscape [START_REF] Fieldsend | A feature rich distance-based many-objective visualisable test problem generator[END_REF], PLOS-net [START_REF] Liefooghe | On Pareto local optimal solutions networks[END_REF] and PLON [START_REF] Fieldsend | Visualising the landscape of multi-objective problems using local optima networks[END_REF] network visualisations. In this work, we focus on box-constrained problems.

Experimental Setup

This section describes our experimental setup covering the approach adopted to generate problems, algorithms and their settings, and performance metrics.

Dataset

We generated 945 problems with different parameter settings. A random latin hypercube sample [START_REF] Carnell | lhs: Latin hypercube samples[END_REF] of size 1 000 was generated, among which 55 problems were discarded due to restrictions in the benchmark generator. For each parameter setting, 30 instances (folds) were independently created using the generator and thus the total number of instances was 945 × 30 = 28 350. Features used to create the problems are provided in Table 1. All 30 folds for a given problem had the same complexity in terms of features. However, folds could be different from each other because of the randomness in the generator. For each fold, we ran different algorithms (one run per instance): NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF], IBEA [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF] with the ε indicator, MOEA/D [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF] with the Chebyshev scalarising function and random search for up to 50 000 evaluations. For a fair comparison, we kept the same initial population for a fold and used the same population size. It was selected based on the number of objectives and is given in Table 2. We employed an out of the shelf implementation with simulated binary crossover and polynomial mutation with probability of 0.8 and 1 n var and distribution index of 20 and 20, respectively; κ = 0.5 for IBEA. For each algorithm, we calculated the normalised hypervolume (hypervolume of final solutions/hypervolume of the Pareto front 7 ) after 5 000, 10 000, 30 000 and 50 000 evaluations for all 30 folds of each problem. Normalised hypervolume values were then averaged for each of the 945 problems and 4 algorithms. The code is available at: https://github.com/ tichugh/Feature_Analysis_DBMOPP_EMO_2023, and the corresponding dataset at: https://doi.org/10.5281/zenodo.7155803.

Algorithm Performance

We analysed results with R [23] using the caret [START_REF] Kuhn | Building predictive models in R using the caret package[END_REF], rpart [START_REF] Therneau | rpart: Recursive partitioning and regression trees[END_REF] and ggplot2 [START_REF] Wickham | ggplot2: Elegant graphics for data analysis[END_REF] packages. In Figure 2 (left), we show the average normalised hypervolume (and the 95% confidence interval) for each algorithm with respect to the search budget over all considered 945 problems. Figure 2 (right) gives the proportion of problems where each algorithm obtained the best average performance, over the 30 folds, for the considered budgets. We observe that IBEA was consistently the best-performing algorithm, whatever the budget, and outperformed others in at least 50% of problems (for a budget of 5 000 evaluations), and at most 80% (10 000 evaluations). It was followed by NSGA-II (almost as good as IBEA for 7 1 000 members drawn from the Pareto front plus all non-dominated points found by the union of the algorithms' approximation sets for each instance. The reference point for hypervolume was 1.1 × maximum of objective values on the Pareto front and estimated via Monte Carlo [START_REF] Everson | Full elite sets for multi-objective optimisation[END_REF] with 50 000 samples for 4+ objectives. 50 000 evaluations). However, one should note that NSGA-II was not significantly better than random search for 5 000 evaluations. MOEA/D was efficient for the budget 5 000, but the increase of budget did not improve its performance as much as for the others. Indeed, the average hypervolume obtained by MOEA/D went from 0.77 for a budget of 5 000 to 0.79 for 50 000. A similar observation was reported in [START_REF] Fieldsend | A visualizable test problem generator for many-objective optimization[END_REF]. We conjecture that MOEA/D is significantly impacted by an increase in local Pareto fronts and multi-modality. Random search was dominated by other algorithms for the two lowest budgets. However, surprisingly, it surpassed MOEA/D for the highest budgets. In fact, random search was even the best for 1 problem for a budget of 30 000, and 2 problems for 50 000. Interestingly, whatever the budget, there is no algorithm that outperforms the others for all problems. For the smallest budgets, IBEA and MOEA/D share the success almost equally on more than 95% of the 945 problems. For the largest budgets, NSGA-II and IBEA share the success on more than 99% of problems.

Experimental Study

This section uses a machine learning perspective to investigate the relationship of problem features and algorithm performance, the predictive power of features for performance prediction, and classification for feature-based algorithm selection.

Problem Features vs. Algorithm Performance

We first investigate how problem features impact search performance. Figure 3 shows how the normalised hypervolume of algorithms is individually impacted by each of the 7 problem features. Due to space restrictions, we report only two budgets. In addition, Figure 4 gives the Spearman's rank correlation coefficient between each problem feature and algorithm performance for all budgets. A larger hypervolume indicates a better performance and thus a positive correlation means that the problem feature has a favourable effect on algorithm performance.

For a given problem feature, the trend is similar for all algorithms and budgets (there is no feature with a positive effect for one algorithm and a negative effect for another algorithm). The same goes for the budgets. However, the strength of correlation is at times different. For instance, although a larger number of objectives (n obj) means a worse performance for all algorithms and budgets, it is more impactful for NSGA-II than for other algorithms for large budgets. The more variables (n var), the worst the performance. However we see in Figure 3 (top-left) that, for a budget of 5 000, NSGA-II is good for a small number of variables, and becomes worse than random search as the number of variables increases. Correspondingly, while the number of objectives and hav- ing non-identical Pareto sets (nonident ps) negatively impact performance, the number of dominance resistance regions (n resist regions) has a positive effect. For a budget of 30 000, NSGA-II is as good as IBEA with few objectives and few dominance resistance regions, but IBEA gets better as both numbers increase. Besides, the fewer local fronts (n local fronts) and disconnected Pareto sets (n discon ps), the better the performance of all algorithms with all budgets. The varying density (var density) has a minor impact on performance.

Overall, problem features in the 945 problems often imply a significant difference between algorithms independently of other features.

Performance Prediction by Regression

The previous results concerned the individual effect of problem features on performance. We now investigate their combined effect by constructing a regression model for predicting the hypervolume reached by the algorithms under different budgets. We thus end up with 4 (algorithms) × 4 (budgets) = 16 models aiming at predicting performance separately for each algorithm and budget, using the problem features as predictors. We consider random forest [START_REF] Breiman | Random forests[END_REF][START_REF] Liaw | Classification and regression by randomforest[END_REF] with default parameters, a well-established state-of the-art ensemble learning method that constructs multiple decision trees for regression. We start by evaluating the prediction accuracy of the trained models using 30 independent replicates of 10-fold cross-validation. We report the repeated cross-validated coefficient of determination (R 2 ) for each algorithm and budget in Figure 5. We observe that the smallest R 2 obtained over all folds and repetitions is above 0.7, and the median R 2 is always above 0.85, whatever the algorithm and budget. This suggests that more than 85% of the variance in hypervolume values across all problems is explained by the model, and thus by problem features. The slight drop in the prediction accuracy for IBEA and NSGA-II as the budget increases is not significant and the R 2 values remain quite satisfactory. The prediction accuracy for random search is particularly high, regardless of the budget.

For each regression model, we also compute the importance of predictors, commonly measured as the mean decrease of prediction accuracy with random forest [START_REF] Breiman | Random forests[END_REF][START_REF] Liaw | Classification and regression by randomforest[END_REF]. The higher the value, the more important the predictor. Figure 6 shows the relative importance of problem features, scaled between 0 and 100. Overall, the most important problem features remain quite consistent with their correlation with algorithm performance reported in Section 4.1. The numbers of variables, objectives, dominance resistance regions and the presence of nonidentical Pareto sets appear on top of the list, whereas the varying density only has a marginal contribution to the prediction accuracy. Interestingly, the presence of non-identical Pareto sets and the number of disconnected Pareto sets get more important as the budget increases. However, noticeable differences appear for MOEA/D, for which the numbers of local fronts and dominance resistance regions are highly important, regardless of the budget. They even surpass the number of variables as the most important features for larger budgets. Thus, it is interesting to note that even though the problem features coming from the problem generator have a similar effect on performance overall, their strength may be quite different depending on the considered algorithm and budget.

Algorithm Selection by Classification

Random Forest Classification. We now focus on feature-based algorithm selection to answer the following question: Given a problem, what is the recommended algorithm for solving it under a particular budget? The interest is no longer in which problem features have more impact on the performance of algorithms, but which features best distinguish algorithms from each other. We still apply random forest, this time for classification. We construct 4 classification models (one per budget) to predict the best algorithm using problem features as predictors. Their classification accuracy, based on 30 repetitions of 10-fold crossvalidation, is reported in Figure 7 (left). The lowest accuracy obtained over all folds and repetitions is 0.64, and the median accuracy is always above 0.75 for all budgets. This means that the classifier is able to predict the best algorithm in at least 75% of the cases, which is significantly better than a random classifier (with an accuracy of 25%), or a dummy classifier that would always select IBEA, the most frequent best algorithm for any budget, which outperformed other algorithms in 53%, 80%, 63%, and 58% of problems respectively, for budgets of 5 000, 10 000, 30 000, and 50 000 evaluations, as reported in Section 3.2. The accuracy of random forest is slightly higher for a budget of 10 000. We attribute this to the fact that IBEA outperforms other approaches more often under this budget, which makes the classification problem easier. We also report in Figure 7 (right) the relative hypervolume deviation of the feature-based random forest classifier (Auto) from the virtual best algorithm, that is the ideal method, an oracle, that always selects the best algorithm. This measure is termed regret, as it indicates how far the obtained hypervolume is from an ideal classifier. It is compared against always selecting each one of the considered algorithms. Notice the log scale in the plot. The regret obtained by the feature-based classifier ranges from 0.0016 to 0.0044 and deviates from an ideal classifier by less than 0.5%. This is less than IBEA, the most frequent best algorithm for all budgets, by an order of magnitude. Compared to an ideal classifier, the relative performance of NSGA-II increases with the budget, while MOEA/D moves away from it. In Figure 8, we report the relative importance of problem features for the random forest classifiers. As suggested earlier, the importance of the number of variables decreases with the budget. By contrast, the importance of the numbers of objectives and of dominance resistance regions increases with the budget, and even exceeds the number of variables for larger budgets. The numbers of disconnected Pareto sets and local fronts are less important for algorithm selection than for performance regression. In fact, they have a low importance similar to the presence of non-identical Pareto sets and to the varying density for some budgets. Thus, although they have a significant effect on algorithm performance, the impact on all algorithms is the same. The feature analysis for classification shows that it is possible to recommend an algorithm based on problem features with a fairly high accuracy. Moreover, the features have a different impact on the choice of the algorithm depending on the budget. Decision trees. We conclude with a basic classifier for algorithm selection based on a decision tree. Its construction follows the well-established CART algorithm [START_REF] Breiman | Classification and regression trees[END_REF][START_REF] Therneau | rpart: Recursive partitioning and regression trees[END_REF], whose segmentation criterion is the Gini diversity index and which generates binary decision trees (i.e. a node has two children at most). In Figure 9, we show decision trees for a budget of 5 000 (left) and 30 000 evaluations (right). Numbers below each node are the number of times NSGA-II, IBEA, MOEA/D and random search are each the best algorithm, respectively, followed by the proportion of problems covered by the node. There are only three values on the first rows for a budget of 5 000 since random search is never the best. Although the accuracy is slightly lower than that of a random forest (0.75 and 0.79, respectively, for a budget of 5 000 and 30 000), we argue that this can provide a useful recommendation system for algorithm selection, with only three levels of decision in these examples.

The tree levels are consistent with the importance of features depicted by random forest, but their joint effect appears more explicitly. Under the smaller budget, the first decision is based on the number of variables. For example, with less than 13 variables, more than 2 local fronts, and less than 10 objectives, IBEA is the best. For a larger budget, and as expected from our previous comments, the feature that appears on top of the tree is the number of objectives. In this case, NSGA-II is the best with few objectives and few variables or few local fronts. Conversely, looking at the rightmost branch of the tree, IBEA is clearly the best with 4+ objectives and 3+ dominance resistance regions. This covers 48% of the problems. Such decision trees justify why an algorithm is recommended based on feature values. In addition, it points out the problem characteristics for which search mechanisms should be refined to improve performance.

Conclusions

We adopted a machine learning perspective to carry out the most extensive feature-based benchmark study of distance-based multi/many-objective optimisation problems to date. We generated 28 350 instances based on 945 problem configurations by varying the complexity controlled by 7 features. Random forests and decision trees were then used to understand correlations between the problem features and algorithm performance, predict algorithm performance, and automate the task of algorithm selection for a given problem and budget at hand. We find that, although the considered problem features affect the performance of algorithms in distinctive ways, when used as predictors in a random forest classifier we can predict the best algorithm with an accuracy of 75% or more. Thus, problem features can control the complexity of a problem, and lend themselves to selecting an algorithm when faced with a previously unseen problem. This is the first automated algorithm selection study for continuous problems with more than 2 objectives. We observed that the number of objectives is (i) negatively correlated with algorithm performance, (ii) one of the most important problem features for predicting algorithm performance, especially for larger budgets, and (iii) a key feature to making an accurate algorithm selection when faced with an unseen problem. Future work could investigate if considering additional problem and landscape features can help increase prediction accuracy further. Some expected algorithm behaviours with respect to specific problem features could also be corroborated based on a fine-grained analysis of the data produced in this work. At last, it would be worth studying the sensitivity of algorithm parameters, and their joint impact with problem features on performance.

Fig. 1 .

 1 Fig. 1. Top-left: an example 3-objective problem with regions creating local fronts (green triangles), dominance resistance regions (points and lines) and Pareto set (red triangle). Top-right: its corresponding local dominance landscape -black regions are locally dominance neutral, shaded grey regions denote basins (for all basin members all neighbouring dominating moves lead to the same dominance neutral region), and white regions denote locations where immediate neighbours lead to different basins (saddles). Bottom-left: PLOS-net and Bottom-right: PLON network visualisations.

Fig. 2 .

 2 Fig.2. Hypervolume (left) and proportion of problems where each algorithm obtains the best average hypervolume (right) with respect to the search budget.

Fig. 3 .

 3 Fig. 3. Hypervolume vs. each problem feature and 5 000 and 30 000 evaluations.

Fig. 4 .

 4 Fig. 4. Correlation between problem features and algorithm performance.

Fig. 5 .

 5 Fig.5. Coefficient of determination (R 2 ) of regression models trained by algorithm and budget, calculated using 10-fold cross-validation with 30 repetitions.

Fig. 6 .

 6 Fig. 6. Relative importance of problem features for regression models.

Fig. 7 .

 7 Fig. 7. Accuracy (left) and regret (right) of classification models trained by budget, calculated using 10-fold cross-validation with 30 repetitions.

Fig. 8 .

 8 Fig. 8. Relative importance of problem features for classification models.

Fig. 9 .

 9 Fig.9. Algorithm selection for budgets of 5 000 (left) and 30 000 evaluations (right).

Table 1 .

 1 Problem features considered from the DBMOPP generator[START_REF] Fieldsend | A visualizable test problem generator for many-objective optimization[END_REF].

	6 Available	in	Matlab	(https://github.com/fieldsend/DBMOPP_generator),
	and	in	Python	(https://github.com/industrial-optimization-group/
	desdeo-problem/tree/master/desdeo_problem/testproblems/DBMOPP).

Table 2 .

 2 Setting of the population size according to the number of objectives.

	number of objectives	2	3	4	5	6	7	8	9 10
	population size	100 105 120 126 132 112 156 90 275

Acknowledgements

This research is part of the thematic research area DEMO, Decision Analytics utilising Causal Models and Multiobjective Optimisation, jyu.fi/demo, at the University of Jyvaskyla.