Bertrand Gauffre
email: bertrand.gauffre@inrae.fr

Alexandre Boissinot

Vincent Quiquempois

Raphael Leblois

Pierre Grillet

PG Sophie Morin

Damien Picard

Cécile Ribout

Olivier Lourdais

Agricultural intensification alters marbled newt genetic diversity and gene flow through density and dispersal reduction Running title: Newts landscape genetic in farmlands

Keywords: Landscape genetics, demographic inferences, dispersal, genetic diversity, microsatellites, Triturus marmoratus

Recent agricultural intensification threatens global biodiversity with amphibians being one of the most impacted groups. Because of their biphasic life cycle, amphibians are particularly vulnerable to habitat loss and fragmentation that often result in small, isolated populations and loss of genetic diversity. Here, we studied how landscape heterogeneity affects genetic diversity, gene flow and demographic parameters in the marbled newt, Triturus marmoratus, over a hedgerow network landscape in Western France. While the northern part of the study area consists of preserved hedged farmland, the southern part was more profoundly converted for intensive arable crops production after WWII. Based on 67 sampled ponds and ten microsatellite loci, we characterized regional population genetic structure and evaluated the correlation between landscape variables and i) local genetic diversity using mixed models and ii) genetic distance using multiple regression methods and commonality analysis. We identified a single genetic population characterized by a spatially heterogeneous isolation-bydistance pattern. Pond density in the surrounding landscape positively affected local genetic diversity while arable crop land cover negatively affected gene flow and connectivity. We used demographic inferences to quantitatively assess differences in effective population density and dispersal between the contrasted landscapes characterizing the northern and southern parts of the study area. Altogether, results suggest recent land conversion affected T. marmoratus through reduction in both effective population density and dispersal due to habitat loss and reduced connectivity.

Introduction

Habitat loss and fragmentation induced by human activities is a major threat to global biodiversity [START_REF] Fisher | Landscape modification and habitat fragmentation: a synthesis[END_REF][START_REF] Foley | Global consequences of land use[END_REF]. This often results in small, isolated populations that are more vulnerable to loss of genetic diversity and fitness decrease through inbreeding and fixation of deleterious alleles [START_REF] Frankham | Genetics and extinction[END_REF]. Patterns of neutral genetic diversity mainly depend on the balance between genetic drift and gene flow [START_REF] Hutchison | Correlation of pairwise genetic and geographic distance measures : Inferring the relative influences of gene flow and drift on the distribution of genetic variability[END_REF]. While decrease in effective population size results in lower genetic diversity due to increased genetic drift, gene flow is a source of genetic variability.

Gene flow homogenizes allele frequencies and maintains population connectivity, which buffers against the negative effects of isolation and inbreeding, and allows for future responses to environmental change [START_REF] Frankham | Genetic Management of Fragmented Animal and Plant Populations[END_REF]. Understanding how genetic diversity and genetic differentiation respond to current landscape structure and past land conversion is critically important for conservation decisions and could allow robust predictions about species responses to global change [START_REF] Nimon | Identification of management units using population genetic data[END_REF][START_REF] Scoble | A case for incorporating phylogeography and landscape genetics into species distribution modelling approaches to improve climate adaptation and conservation planning[END_REF][START_REF] Sgrò | Building evolutionary resilience for conserving biodiversity under climate change[END_REF].

In European agricultural hedgerow landscapes, land conversion has caused loss of permanent habitats and their connection in space [START_REF] Benton | Farmland biodiversity: Is habitat heterogeneity the key?[END_REF]. While land conversion threatens all taxa, amphibians constitute one of the most impacted ones [START_REF] Ceballos | Vertebrates on the brink as indicators of biological annihilation and the sixth mass extinction[END_REF][START_REF] Crawford | Altered functional connectivity and genetic diversity of a threatened salamander in an agroecosystem[END_REF][START_REF] Cushman | Effects of habitat loss and fragmentation on amphibians: A review and prospectus[END_REF][START_REF] Hof | Additive threats from pathogens, climate and land-use change for global amphibian diversity[END_REF][START_REF] Mccartney-Melstad | Amphibian molecular ecology and how it has informed conservation[END_REF]. This sensitivity is primarily due to i) their specific habitat requirements and complex life cycle involving spatially distinct breeding and foraging habitats [START_REF] Karraker | Amphibian production in forested landscapes in relation to wetland hydroperiod: A case study of vernal pools and beaver ponds[END_REF][START_REF] Sztatecsny | The abundance of premetamorphic newts (Triturus cristatus, T. marmoratus) as a function of habitat determinants: An a priori model selection approach[END_REF], and ii) their low dispersal capacity [START_REF] Hillman | Physiological vagility and its relationship to dispersal and neutral genetic heterogeneity in vertebrates[END_REF][START_REF] Smith | Dispersal and the metapopulation in amphibian and paradigm ecology are all amphibian conservation: populations metapopulations?[END_REF]. Among amphibians, urodeles (newts and salamanders) have lower mobility [START_REF] Smith | Dispersal and the metapopulation in amphibian and paradigm ecology are all amphibian conservation: populations metapopulations?[END_REF], and are particularly vulnerable to temperature and water constraints [START_REF] Riddell | Thermal cues drive plasticity of desiccation resistance in montane salamanders with implications for climate change[END_REF][START_REF] Riddell | Terrestrial Salamanders maintain habitat suitability under climate change despite trade-offs between water loss and gas exchange[END_REF]. Most European newt species live in water during their larval stage and for reproduction but they are generally terrestrial for the rest of the year. Newts mainly reach breeding ponds following corridors connecting ponds with woodlands, and avoid arable crops [START_REF] Marty | Directionality of pre-and postbreeding migrations of a marbled newt population (Triturus marmoratus): Implications for buffer zone management[END_REF]. Therefore, land conversion for intensive agriculture, including the loss of hedgerows, shelters and corridors, could affect populations by decreasing effective population size and connectivity, both resulting in loss of genetic diversity and gene flow.

Landscape genetics combines tools from population genetics, spatial statistics and landscape ecology to relate landscape features directly to population structure, genetic diversity, and gene flow [START_REF] Manel | Ten years of landscape genetics[END_REF][START_REF] Manel | Landscape genetics: combining landscape ecology and population genetics[END_REF][START_REF] Storfer | Putting the 'landscape' in landscape genetics[END_REF]. Detection of barriers and genetic clusters is a first step in landscape genetics but, because many populations are continuously distributed, barriers may only exist at large biogeographic scales. However, at smaller spatial scales, gene flow and functional connectivity might be gradually modulated by landscape and environmental heterogeneity (isolation by resistance, IBR, McRae, 2006). In addition, because species dispersal abilities are generally limited, geographic distance is expected to play a significant role in the explanation of genetic differentiation (isolation by distance, IBD, [START_REF] Rousset | Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance[END_REF][START_REF] Wright | Isolation by distance[END_REF]. IBD is generally considered as the null model in landscape genetics and methods such as causal modelling are employed to compare different models of landscape resistance and IBD [START_REF] Cushman | Gene flow in complex landscapes: testing multiple hypotheses with causal modeling[END_REF]. More recently, multiple regression methods are becoming more popular as they allow to identify landscape variables influencing genetic differentiation in a more robust statistical framework [START_REF] Shirk | A comparison of regression methods for model selection in individual-based landscape genetic analysis[END_REF][START_REF] Wang | Examining the full effects of landscape heterogeneity on spatial genetic variation: A multiple matrix regression approach for quantifying geographic and ecological isolation[END_REF], and rank them according to their importance in shaping genetic structure while accounting for correlations among them [START_REF] Prunier | Multicollinearity in spatial genetics: separating the wheat from the chaff using commonality analyses[END_REF]. While matrix-based analyses (i.e. relating genetic differentiation to landscape structure) are the most popular in landscape genetics, node-based analyses (i.e. relating local genetic diversity to surrounding landscape structure) might provide complementary information on landscape genetics relationships [START_REF] Flavenot | Impact of quarrying on genetic diversity: an approach across landscapes and over time[END_REF]. However, both approaches are mostly correlative and only allow identifying landscape genetics relationships without assessing the demographic parameters (e.g. effective population size and dispersal) by which landscape influences genetic patterns. This is particularly important regarding two related issues. First, landscape genetics is mostly used in conservation studies where populations are rarely at a demographic equilibrium [START_REF] Segelbacher | Applications of landscape genetics in conservation biology: Concepts and challenges[END_REF]. Furthermore, genetic patterns at a large spatial scale might be integrative of many generations that are likely to have experienced numerous changes in population size and dispersal. It is thus challenging to assess whether landscape genetic patterns result from variations in effective population size or changes in dispersal [START_REF] Richardson | Navigating the pitfalls and promise of landscape genetics[END_REF]. Landscape genetics studies often interpret increases in measures of genetic differentiation in terms of reduced gene flow and loss functional connectivity, neglecting the effect of local genetic drift on genetic structure. Thus, demographic inferences of population parameters, i.e. effective population size and dispersal rate and distance, in a spatially explicit framework are required to improve the landscape genetics toolbox. Indeed, methods based on models specifying population size and dispersal function might allow inferring whether genetic diversity (or genetic differentiation) is affected by variations in population size, variations in dispersal or both. The second issue is methodological. There is actually a large body of methods in population genetics producing demographic inferences on Ne and dispersal (migration) rate (m) but they are based on the simple, not spatially explicit, island model, while most natural populations (and landscape genetics approaches) follow the IBD model. There is a real lack of use of demographic inference methods based on the IBD model, including methods within the ABC framework or likelihood based methods [START_REF] Bertorelle | ABC as a flexible framework to estimate demography over space and time: Some cons, many pros[END_REF][START_REF] Rousset | Likelihood-based inferences under isolation by distance: Two-dimensional habitats and confidence intervals[END_REF]. To investigate landscape genetics effects, comparative approaches involving multiple landscapes or gradients of landscape composition are recommended [START_REF] Goldberg | Comparative landscape genetics of two pond-breeding amphibian species in a highly modified agricultural landscape[END_REF]. These approaches should be particularly appropriate to infer and compare demographic parameter estimates in contrasted landscapes, which can provide further information for conservation efforts [START_REF] Beebee | The amphibian decline crisis: A watershed for conservation biology?[END_REF].

In this study, we investigated the effect of landscape composition and land conversion on genetic diversity, genetic structure and demographic parameters in the marbled newt, Triturus marmoratus (Latreille 1800). This endangered species (listed on Annex III of the Bern Convention and Annex IV of the European Habitats Directive) is distributed over the north Iberian Peninsula and western France. It prefers areas with bushes, hedgerows and trees, and avoids pastures and open areas [START_REF] Jehle | The terrestrial summer habitat of radio-traked great crested newts (Triturus cristatus) and marbled newts(T.marmoratus)[END_REF][START_REF] Jehle | Post-breeding migrations of newts (Triturus cristatus and T. marmoratus) with contrasting ecological requirements[END_REF][START_REF] Trochet | Postbreeding movements in Marbled newts ( Caudata , Salamandridae ): A comparative radiotracking study in two habitat types[END_REF]. In contrast to those sampled in agricultural lands, newts in woodlands moves greater distances [START_REF] Trochet | Postbreeding movements in Marbled newts ( Caudata , Salamandridae ): A comparative radiotracking study in two habitat types[END_REF]. Therefore, T. marmoratus might be particularly impacted by agricultural intensification. More specifically, we hypothesize a negative influence of habitat loss with the conversion of permanent grasslands and meadows to arable crops on T. marmoratus population size and connectivity as recently suggested [START_REF] Costanzi | Agricultural landscapes and the Loire River influence the genetic structure of the marbled newt in Western France[END_REF].

Determining which landscape features influence T. marmoratus population genetic patterns would help to understand population connectivity and to identify key habitats for preserving genetic variation.

We studied T. marmoratus in a 6000 km² agricultural landscape in Western France. The study area is heterogeneous with preserved hedgerow farmland landscape with permanent grasslands and meadows (hereafter grasslands) in the northern part, while recently opened landscape dedicated to arable crops (mostly cereals) dominate the southern part (see Figure 1). We sampled 67 ponds distributed in these two contexts. Previous studies showed that pond disappearance was associated with a decrease in grasslands and an increase in arable crop land covers [START_REF] Curado | Amphibian pond loss as a function of landscape change -A case study over three decades in an agricultural area of northern France[END_REF]. In the study area, ponds were initially created to water the livestock and are thus associated to grasslands. When grasslands are converted to arable crops ponds are usually drained to increase cultivated land [START_REF] Arntzen | Amphibian decline, pond loss and reduced population connectivity under agricultural intensification over a 38 year period[END_REF]. We characterized landscape changes (hedgerow loss) over the last 60 years in the study area using historical landscape data and analyzed genetic data collected in order to address the following goals:

(i) Identify the effect of landscape composition at various spatial scales on within pond genetic diversity. We expected genetic diversity to increase with pond density and forested land cover (woodlands and hedgerows) in the local landscape and to decrease with arable crop land cover;

(ii) Characterize population genetic structure and test whether it was homogeneous across the study area. In particular, we expected the northern and southern parts of the study area to show contrasted levels of genetic structuring with stronger IBD in the southern part;

(iii) Identify landscape features that affect genetic differentiation. We expected roads, arable crops and urbanized land cover to negatively affect gene flow while the opposite was expected for pond density and forested land cover.

(iv) Infer population demographic parameters, including effective population density (equivalent to effective population size in the continuous population model), dispersal rate and the shape of dispersal distribution between contrasted landscape contexts. Enhanced dispersal and larger population effective densities were expected in more connected areas.

Materials and methods

Study area and sample collection

The study area covers the entire Deux-Sèvres department (6000 km²), in Western France. It is a farmland area with the northern part dominated by preserved hedged farmland with permanent grasslands and meadows dedicated to livestock farming while the southern part is more open and mostly dedicated to intensive cereal crops production (Figure 1). Sampling occurred during the newt breeding season from February to May, in 2014 and 2015. A total of 734 individuals from 67 ponds were sampled (15 ponds in 2014 and 52 in 2015). At night, individuals were localized using a flashlight and captured using a dip net. Newts were then put into a bucket half filled with water. DNA samples were obtained by rubbing the bottom of the mouth, cheeks and tongue with cotton swab to collect epithelial cells [START_REF] Pidancier | Buccal swabs as a non-destructive tissue sampling method for DNA analysis in amphibians[END_REF]. Cotton swabs were stored in an alcohol solution (solution EDTA buffer: Ethylene Diamine Tetra-Acetic). All individuals were released once the pond sampling was completed.

Genetic data collection

DNA extraction from cotton swabs was performed using Macherey-Nagel Genomic DNA extraction kit. We used a set of 10 microsatellite markers (see [START_REF] Costanzi | Agricultural landscapes and the Loire River influence the genetic structure of the marbled newt in Western France[END_REF], including nine developed specifically for the marbled newt [START_REF] Costanzi | Characterization of nine new microsatellite loci for the marbled newt, Triturus marmoratus[END_REF] and one (Tcri27) originally developed for the crested newt [START_REF] Dopaminergic Krupa | Microsatellite loci in the crested newt (Triturus cristatus) and their utility in other ttax[END_REF]. Details of microsatellites amplification are presented in Appendix S1 in Supplementary Information. Alleles were scored using GENEMAPPER TM v 4.0 (Applied Biosystems) and checked manually.

There was a relatively high amount of missing data (13.41 % over the 734 samples) related to low yield of DNA extraction in some pond samples, likely due to low quality DNA conservation. A total of 612 individuals from 64 ponds were successfully genotyped at 8 or more loci. From this global dataset, ponds less than 150 m apart were pooled and ponds with less than 7 individuals were removed. The final dataset included 550 individuals from 39 ponds (Table S1.2, Supplementary Information). For sake of consistency, node-based and matrix-based landscape genetics analyses as well as demographic inferences were all performed on this population-based dataset.

Genetic diversity

We tested linkage disequilibrium between all pairs of loci on the global 612 individuals dataset using GENEPOP4.2 [START_REF] Rousset | GENEPOP'007: A complete re-implementation of the GENEPOP software for Windows and Linux[END_REF]. Tests for deviations from Hardy-Weinberg equilibrium for each locus and linkage equilibria for each locus pair were performed using GENEPOP globally, and in each of the 39 ponds. Significance of tests was assessed following a false discovery rate correction for multiple tests [START_REF] Benjamini | Controlling the false discovery rate : A practical and powerful approach to multiple testing[END_REF] with a nominal significance level of 5%. Observed and expected heterozygosities (H o and H e ) and Weir and Cockerham's estimate of F IS [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] were calculated using GENEPOP globally and for each pond. Using the rarefaction procedure implemented in FSTAT 2.9.3.2 [START_REF] El Mousadik | High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L.) Skeels] endemic to Morocco[END_REF][START_REF] Goudet | FSTAT (Version 1.2): A computer program to calculate F-statistics[END_REF], we calculated allelic richness corrected for sample size (A r ) for each locus in all 39 pond.

Landscape features and changes overtime

We characterized the landscape over the study area using CORINE Land Cover 2012 (resolution 1/100 000) and BDTOPO® 2015 (® IGN) in QGIS3.0.1 (QGISDevelopment Team, 2018). Eight landscape features were selected based on their functional role for T.

marmoratus:

(i) The area of five classes of land cover were extracted, including woodlands, which represents the terrestrial habitat of T. marmoratus, hedgerows and grasslands, which might constitute alternative terrestrial habitats and favorable habitats for T. marmoratus movement, and arable crops and urbanized land covers, which are expected to constitute unfavorable features for T. marmoratus.

(ii) Rivers and roads (linear features) were extracted, as they can constitute barriers to movements of some amphibians.

(iii) Ponds, which represents the breeding habitat of T. marmoratus were exhaustively inventoried in the Deux-Sèvres department (N = 17400) and manually georeferenced for the purpose of this study.

In order to evaluate habitat change since 1950 we also characterized past and present landscape structure based on hedgerows linear digitalization from aerial photos taken in 1950 and 2015, hereafter Hedgerow 1950 and Hedgerow 2015 . This historical analysis was focused on the landscape surrounding the 39 ponds. We manually digitalized linear distance of hedgerows within 2 km buffers around the ponds. Our approach was restricted to hedgerows because some parameters (such as pond density) were not available on historical images. Yet, hedgerows are integrative of landscape quality and correlated with pond density and meadows. Indeed, Pearson correlation coefficient calculated over the entire Deux-Sèvres department with a 5 km grid indicated strong correlation between hedgerow and pond density and between hedgerows and meadows (Pearson's r = 0.85, P-value < 0.001 and Pearson's r = 0.91, P-value < 0.001, respectively).

Effects of local landscape composition on genetic diversity

Local genetic diversity depends on both population size and landscape connectivity. We investigated which landscape features affected the genetic diversity of T. marmoratus. Within circular buffers centered on the 39 ponds, we calculated the area of four different classes of land cover (woodlands, grasslands, arable crops and urbanized), the length of roads and rivers, and the number of ponds. Considering hedgerows, we used the linear data digitalized in 2 km buffers around the ponds (Hedgerow 2015 ) as it was more precise compared to the land cover database and allowed us to include historical data to derive an index of landscape intensification (Intensification index ) as 1-(Hedgerow 2015 / Hedgerow 1950 ). Intensification index was fixed to 0 for ponds located in large forest patches, where woodland represented > 75% of the total buffers area, because there are no hedgerows in those patches. There was no intensification index < 0 (i.e. no increase in hedgerow from 1950 to 2015 in any buffers). To explore the scale at which landscape features influence genetic diversity, we calculated landscape composition in buffer of different radius (500, 1000 and 2000 m). The maximal distance was set to 2000 m to avoid artificially large correlation due to overlapping of buffers and because buffers tended to be out of the mapped area at larger scales. To address possible confounding effects due to correlations between landscape features, we computed pairwise Pearson correlations between landscape features for each buffer radius. Pairs of landscape features with Pearson correlation above the 0.70 threshold were considered highly correlated [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF]. We thus excluded Hedgerow 2015 from the set of landscape descriptors for all models, grasslands from the set of landscape descriptors for 500 model, and Intensification index from the set of landscape descriptors for 1000 m model (see details on correlations in Appendix S2, Supplementary Information). A r and H o , calculated for each locus and each pond 39 ponds, were alternatively used as response variable. A r is expected to be more sensitive to recent reductions in population size since it is more impacted by loss of rare allele than H o [START_REF] Schwartz | Genetic monitoring as a promising tool for conservation and management[END_REF].

We used linear mixed effect models to quantify the contribution of landscape features on genetic diversity using the 'lmer' function implemented in the LME4 package [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF] in R 4.0.2 (R Development Core Team, 2014). We treated selected landscape features as fixed factors and they were z-transformed (i.e. standardized to zero mean and a standard deviation of one). Locus identity was included as a random intercept to account for differences in allelic diversity among loci. First, we compared R² of the full models (interactions were not considered) at 500, 1000 and 2000 m to assess the best spatial scale. Second, for the best spatial scale, we ran models with all possible combinations of explanatory variables using dredge function implemented in the MUMIN package [START_REF] Barton | MuMIn: Multi-model inference[END_REF]. Models were fitted with maximum likelihood to compute reliable Akaike information criterion (AIC) scores, and ranked according to AICc values (the small-sample-size corrected version of AIC). To take into account uncertainty in model selection we then used model averaging on models with Delta AICc values < 4 to get estimates of final parameters. Models were validated a posteriori by checking plots of residuals. We also checked for autocorrelation in the model residuals.

Regional genetic structure

To assess the extent to which ponds might be isolated from each other and form separate populations, the population genetic structure was investigated by two complementary approaches.

(i) We used the Bayesian clustering method implemented in STRUCTURE Version 2.2 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF], with the admixture model and correlated allele frequencies [START_REF] Falush | Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies[END_REF] to determine the number of genetic clusters (K). Ten replicate runs of 2 × 10 6 Markov chain Monte Carlo (MCMC) iterations, after an initial burn-in period of 5 × 10 5 iterations, were performed for values of K ranging from one to ten. Results were summarized using the standard pipeline on the CLUMPAK web server [START_REF] Kopelman | CLUMPAK : a program for identifying clustering modes and packaging population structure inferences across K[END_REF]. The most likely number of clusters (K) was explored using the estimated logarithm of likelihood (LnP(D)) and the [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: a simulation study[END_REF] ΔK method that finds the point of greatest change in the distribution of LnP(D) with STRUCTURE HARVESTER Version 0.6.92 [START_REF] Earl | STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF]. We used all 612 individuals for this individual-based analysis.

(ii) The level of pairwise genetic differentiation between the 39 ponds was quantified by pairwise F ST [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF], and tested using the exact probability test for population differentiation implemented in GENEPOP. IBD was then analyzed by regressing pairwise estimates of F ST /(1-F ST ) against the logarithm of the geographical distances [START_REF] Rousset | Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance[END_REF], and tested using a Mantel test (10,000 permutations). To visualize IBD and whether it was homogeneous across our study area, we used LOCALDIFF 1.5 [START_REF] Duforet-Frebourg | Nonstationary patterns of isolation-bydistance: Inferring measures of local genetic differentiation with bayesian kriging[END_REF], a method that infers local genetic differentiation based on Bayesian kriging. The pairwise F ST /(1-F ST ) values between the 39 ponds was used as input measure of pairwise dissimilarity. The method requires fictive neighboring populations to be introduced at the vicinity of sampling sites (i.e. short distance compared to the dimension of the region under study), as a means to provide measures of local genetic differentiation that are comparable between sampling sites. Hence, four fictive neighboring populations were introduced at 2 km from the 39 ponds. The output of LOCALDIFF generates a minimum convex polygon that encompasses all sampling sites, where warmer colors represent areas with stronger local genetic differentiation. We used Pearson correlations to test the relationship between estimates of genetic diversity (A r , H e and H o ) within ponds and estimates of local genetic differentiation computed with LOCALDIFF.

Effects of landscape composition on genetic structure

To investigate the effect of landscape features on spatial genetic structure, we performed multiple regression of distance matrices (MRDM) [START_REF] Legendre | Modeling brain evolution from behavior: A permutational regression approach[END_REF] and multiple matrix regression with randomization MMRR [START_REF] Wang | Examining the full effects of landscape heterogeneity on spatial genetic variation: A multiple matrix regression approach for quantifying geographic and ecological isolation[END_REF]. Pairwise genetic distances were computed between the 39 ponds using F ST /(1-F ST ) as genetic distance. To compute pairwise resistance distances a specific layer was created for each landscape feature. We overlaid a 250 m grid on these raster data, and calculated the percentage per grid cell of each categorical feature (Woodlands, Hedgerows, Arable crops, Urbanized, and Grasslands), attributed a value of one for each grid cell intersected by linear features (roads and rivers) and calculated the number of ponds in each cell of the grid. To allow comparison among landscape features, these layers were finally rescaled to range from 1 to 100 and used in CIRCUITSCAPE 4.0.5 [START_REF] Mcrae | Circuitscape 4 user guide[END_REF] to compute pairwise resistance distances between ponds.

MRDMs are similar to classical multiple ordinary least-square regressions, except that the significance of model fit and beta weights (β) is assessed through permutations [START_REF] Legendre | Modeling brain evolution from behavior: A permutational regression approach[END_REF]. Beta weights could be heavily impacted by multi-collinearity among landscape variables (Ray-Mukherjee et al., 2014), so we used a commonality analysis (CA) (J. G. [START_REF] Prunier | Multicollinearity in spatial genetics: separating the wheat from the chaff using commonality analyses[END_REF][START_REF] Ray-Mukherjee | Using commonality analysis in multiple regressions: A tool to decompose regression effects in the face of multicollinearity[END_REF] to dissect the complexity of landscape features relative contribution to the model fit. CA is a variance-partitioning procedure that disentangles the individual vs. shared contribution of each predictor to R². This procedure is particularly useful when the predictors are themselves correlated [START_REF] Prunier | Multicollinearity in spatial genetics: separating the wheat from the chaff using commonality analyses[END_REF].

Pearson's correlations among pairwise resistance distances computed for each landscape feature ranged from 0.02 to 0.69, while variance inflation factor (VIF) ranged from 1.62 to 2.96 (Table S4 in Supplementary Information). All pairwise resistance distances were kept in the model as we considered Pearson's correlations and VIF to be high above thresholds of 0.7 and 10, respectively [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF][START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF]. To determine the contribution of pairwise resistance distances relatively to IBD, pairwise geographic distance was also included in the model. CA-MRDM was run using the R packages ECODIST [START_REF] Goslee | The ecodist package for dissimilarity-based analysis of ecological data[END_REF] and YHAT [START_REF] Nimon | Identification of management units using population genetic data[END_REF] following the CAonDM script provided by [START_REF] Prunier | Multicollinearity in spatial genetics: separating the wheat from the chaff using commonality analyses[END_REF]. Significant levels were assessed with 10,000 permutations after sequential Bonferroni correction [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF]. Commonalities and 95% confidence intervals were computed using a bootstrap procedure with 1000 replicates based on a random selection of 10% of samples (out of 39 populations) without replacement [START_REF] Prunier | Multicollinearity in spatial genetics: separating the wheat from the chaff using commonality analyses[END_REF]. Alternatively, we ran MMRR according to [START_REF] Wang | Examining the full effects of landscape heterogeneity on spatial genetic variation: A multiple matrix regression approach for quantifying geographic and ecological isolation[END_REF]. MMRR was implemented in R using the package PopGenReport [START_REF] Adamack | PopGenReport: Simplifying basic population genetic analyses in R[END_REF]. We used pairwise geographic and the eight pairwise resistance distances as the explanatory variables and 999 permutations to assess the additive effect of both independent factors (geographical distance and landscape features). All matrices were standardized using the 'scale' function implemented in R before running the MMRR. Finally, we tested each resistance distances and geographic distance in eight separate models.

Inferences on demographic parameters

Spatial variation in IBD was explored by the maximum likelihood method implemented in MIGRAINE that infers model parameters using importance sampling algorithms (De Iorio, Griffiths, Leblois, & [START_REF] De Iorio | Stepwise mutation likelihood computation by sequential importance sampling in subdivided population models[END_REF], extended to consider IBD as a model for population structure [START_REF] Rousset | Likelihood-based inferences under isolation by distance: Two-dimensional habitats and confidence intervals[END_REF][START_REF] Rousset | Likelihood and approximate likelihood analyses of genetic structure in a linear habitat: Performance and robustness to model mis-specification[END_REF]. Based on LOCALDIFF results, ponds from the northern and southern parts of the study area were analyzed separately to infer area specific demographic parameters. Our specific aim here, was to assess whether differences in IBD between these two parts of the study area result from differences in effective population size and / or differences in dispersal. We excluded the three northermost ponds to get parameter estimates from a more homogeneous landscape and with a continuous distribution of sampled ponds. Based on a geometric distribution for dispersal and a K allele model for mutation, MIGRAINE provides point estimates, 95% coverage confidence intervals (CIs) and two-dimensional parameter likelihood profiles for several parameters: 1) The scaled local population size θ=2×N genes ×μ, where N genes is the local population size expressed in number of genes and μ the mutation rate per locus per generation. 2) The scaled emigration rate per generation: γ = 2×N genes ×m, where m is the total emigration rate per generation for a local population. 3) The parameter of the geometric distribution of dispersal (g). 4) The neighborhood size, Nb = 4π×N×σ 2 , where N is the effective population size and σ 2 the mean squared parent-offspring dispersal distance. All MIGRAINE runs were performed under a twodimensional model of IBD. The spatial binning of samples was performed using MIGRAINE-GUI and resulted in 16 x 10 and 20 x 11 bins for the northern and southern parts of the study area, respectively. Bin size, 3275 x 3275 m, was set constant over the two analyses to facilitate results interpretation. We used the following computing parameters: 2000 trees, 100 points and 3 iterations. We translated the parameters inferred from MIGRAINE into effective population size using a mutation rate range commonly used for microsatellites: 1×10 -4 [5×10 -3 -5×10 -5 ] [START_REF] Sun | A direct characterization of human mutation based on microsatellites[END_REF].

Results

1) Landscape characteristics and changes through time

Pond density was almost five times higher in the northern part with 3.79 ponds / km² compared to 0.8 ponds / km² in the southern part of the study area (Figure 1). Consistently, arable crops represented 41.16 % and 63.88 % and meadows 48.59 % and 21.49 % of land cover in the northern and southern parts of the study area, respectively. Based on the digitalized hedgerow in 2 km buffers around the 39 ponds, we found a significant difference in the proportion of hedgerow loss between the northern and southern parts of the study area (two samples t-test, P-value < 0.001). On average loss was 34 % and 53 % in the northern and southern parts of the study area, respectively. However, no significant effect of landscape composition on H o was found at any spatial scale (Table S2.5, Supplementary Information).

4) Regional pattern of genetic structure

The most likely value of K from the STRUCTURE analysis based on the method of [START_REF] Evanno | Detecting the number of clusters of individuals using the software structure: a simulation study[END_REF] was two (Table S3, Supplementary Information). Consistently, the estimated logarithm of likelihood for data was highest for K = 2. However, for most individuals, the estimated membership coefficients in each cluster Q was low (i.e. > 60 % individuals had Q < 0.7) and inspection of STRUCTURE barplot was more consistent with an IBD pattern than to the presence of two or more genetic clusters (Figure S1). Indeed, a pattern of IBD was supported by the positive relationship between F ST /(1-F ST ) and the logarithm of geographic distance (Figure 3, Slope = 0.027, 95% CI: (0.018 -0.037), Mantel test P-value < 0.001) over the 741 pairs of ponds. The average F ST value between ponds was 0.066 and for 621 pairs (84%) genetic differentiation was significant. The analysis performed with LOCALDIFF indicated a heterogeneous pattern of IBD across the study area. The southern part was characterized by strongest values of local genetic differentiation, whereas smaller local genetic differentiation was found in the northern part of the study area (Figure 4). We measured a significant negative correlation between allelic richness (Ar) and estimates of local genetic differentiation computed with LOCALDIFF (cor = -0.59, P-value < 0.01), while the correlation was nonsignificant when H o (cor = -0.19, P-value = 0.24) was used as a measure of genetic diversity (Figure 5).

5) Effect of landscape composition on genetic structure

The MRDM was significant and explained 35.7% of the variance in F ST /(1-F ST ) (Table 2).

After sequential Bonferroni corrections, only arable crop land cover had a significant positive effect on genetic differentiation (β crops = 0.33, P-value = 0.016). F ST /(1-F ST ) increased by 0.33 standard deviations with a one standard deviation change in arable crops resistance distance, all other predictors being held constant. Commonality coefficients showed that arable crop land cover uniquely contributed 8.2% of the total variance in F ST /(1-F ST ) and to 22.8% of the 35.7% of the variance explained by the regression model. Unique contributions of other landscape features were negligible or counterbalanced by their common contribution with other predictors (indicative of classical suppressor effect according to Prunier et al. 2017) (Table 2).

The MMRR model including geographic distance and all resistance distances explained nearly 30 % of the variability in genetic differentiation (R² = 0.291, P < 0.001). It revealed that arable crop land cover was the only landscape feature significantly influencing genetic differentiation (β crops = 0.009, P = 0.033, Table 2 and Figure S2, Supplementary Information).

6) Demographic inferences

Outputs from the MIGRAINE software confirmed a clear contrast in demographic parameters between the northern and southern parts of the study area (Table 3). Nb was found to be significantly higher in the northern part (557.5, 95% CI (206.8 -22134)) compared to the southern part (68.77,(95)(96)(97)(98)(99)(100)(101)(102)(103)(104)(105)(106)(107)2)) of the study area as was g, the parameter of the geometric distribution describing dispersal (Table 3). Using 10 -4 as mutation rate for microsatellites it was possible to derive a theoretical area specific effective population density (De) per bin and migration rate among bins (m) (see Table 3 for figures) and draw a theoretical dispersal kernel (Figure 6). The latter indicated higher dispersal rate and distance in the northern part of the study area.

Discussion

In this study, we identified landscape effects on genetic structure, genetic diversity and gene flow in T. marmoratus an amphibian species with supposed limited mobility. Our results support the hypothesis that land conversion (i.e. loss of habitats and their connectivity) that occurred in the last 60 years impacted genetic diversity and functional connectivity. In addition, we were able to conduct demographic inferences to characterize effective population density and dispersal according to landscape context. Both demographic parameters varied with landscape composition.

Population genetic structure

In Western France, at a local spatial scale (up to 6.5 km) previous genetic studies on T. marmoratus reported distinct genetic clusters with STRUCTURE, an average pairwise F ST of 0.11 (range: 0.007 -0.303) and no significant IBD (Jehle, Burke, & Arntzen, 2005;[START_REF] Jehle | Contemporary gene flow and the spatio-temporal genetic structure of subdivided newt populations (Triturus cristatus, T. marmoratus)[END_REF]. In our study, at distance < 6.5 km, the average pairwise F ST was only 0.025 (range: -0.03 -0.07) with 9 out of 25 pairs of ponds significantly differentiated. A single population exhibiting IBD was found, and genetic structure was not affected by rivers or motorways crossing the study area, similar to other newt studies [START_REF] Costanzi | Agricultural landscapes and the Loire River influence the genetic structure of the marbled newt in Western France[END_REF][START_REF] Luqman | No distinct barrier effects of highways and a wide river on the genetic structure of the Alpine newt (Ichthyosaura alpestris) in densely settled landscapes[END_REF][START_REF] Prunier | A 40-yearold divided highway does not prevent gene flow in the alpine newt Ichthyosaura alpestris[END_REF]. Based on 11 T. marmoratus populations sampled across Western France, [START_REF] Costanzi | Agricultural landscapes and the Loire River influence the genetic structure of the marbled newt in Western France[END_REF] identified strong genetic structure at a large scale (> 100 km), including distinct genetic clusters in the area of our study. Their result is likely due to IBD which might create spurious genetic clusters when geographical sampling is clumped [START_REF] Blair | A simulation-based evaluation of methods for inferring linear barriers to gene flow[END_REF][START_REF] Frantz | Using spatial Bayesian methods to determine the genetic structure of a continuously distributed population: clusters or isolation by distance[END_REF]. This supports that continuous distribution of samples with individual-based sampling is the best strategy to uncover unbiased spatial genetic structure. Our results also side with the conclusion of [START_REF] Smith | Dispersal and the metapopulation in amphibian and paradigm ecology are all amphibian conservation: populations metapopulations?[END_REF] that although amphibians are predominantly philopatric with poor dispersal capacities, they could move distances much greater than anticipated. Another relevant finding of our study is the clear heterogeneous IBD pattern across the study area as evidences by the LOCALDIFF analysis. Results from our landscape genetics analyses strongly suggest this contrast might be related to differences in landscape composition between the northern and southern parts of the study area as previously suggested [START_REF] Costanzi | Agricultural landscapes and the Loire River influence the genetic structure of the marbled newt in Western France[END_REF].

Landscape influence on genetic diversity and genetic structure

Genetic diversity is affected by effective population size and connectivity [START_REF] Flavenot | Impact of quarrying on genetic diversity: an approach across landscapes and over time[END_REF]. It is an important parameter in conservation genetics since reduced genetic diversity could translate to lower fitness and subsequently vortex of extinction [START_REF] Fagan | Quantifying the extinction vortex[END_REF]. A review analyzed 19 studies that directly quantified genetic diversity -fitness relationships in amphibians, among which 15 provided evidence that levels of genetic diversity affected important traits such as growth or survival [START_REF] Allentoft | Global amphibian declines, loss of genetic diversity and fitness: A review[END_REF].

Pond density was the only landscape feature affecting local genetic diversity in our study with a positive influence on allelic richness. Interestingly, we did not find any effect of woodlands land cover, the main terrestrial habitat of the species. This result emphasizes the importance of the breeding habitat in newt life cycle and suggests larger effective population density with increased pond density. Higher pond density may also facilitate dispersal and resulting gene flow, as it is easier for a dispersing newt to encounter a new pond. This result might be related to reduction of pond density associated with intensification that affected both parts of the study area, though to a higher degree in the south. Both MRDM and MMRR identified arable crop land cover as the only landscape variable affecting significantly gene flow, likely through reduced functional connectivity. It is also supported by the LOCALDIFF analysis that showed stronger IBD in the southern part of the study area. It is interesting to note that nodeand matrix-based landscape genetics analyses provided complementary results. In the nodebased analysis, we only considered the local landscape (i.e. up to 2 km² buffers around sampled ponds). The landscape feature selected in the model, pond density, was more representative of local landscape quality in terms of habitat availability. For the matrix-based analyses, we considered the full landscape matrix and the landscape feature selected in MRDM and MMRR models, arable crop land cover, was more indicative of reduced landscape permeability between habitat patches. However, since gene flow depends on both effective population density and dispersal, population decline in the southern side of the study area could also explain the negative effect of arable crop land cover on genetic differentiation found with MRDM and MMRR. The negative relationship between genetic diversity and estimates of local genetic differentiation from LOCALDIFF supports the hypothesis that genetic drift is a predominant micro-evolutionary process driving genetic differentiation [START_REF] Coleman | Balancing genetic uniqueness and genetic variation in determining conservation and translocation strategies: A comprehensive case study of threatened dwarf galaxias, Galaxiella pusilla (Mack) (Pisces: Galaxiidae)[END_REF]. Altogether, our study exemplifies the difficulty to disentangle the influence of effective population density and dispersal on gene flow from patterns of genetic diversity and genetic structure. More specifically, it is challenging to assess whether landscape structure and land conversion affect natural populations through reduction in effective population size or dispersal.

Demographic inferences on dispersal and population effective density

Several lines of evidence support a recent demographic decline in the southern part of the study area. Landscape influence on local genetic diversity was only found for allelic richness, while it was not significant for heterozygosity that is less affected by recent demographic changes [START_REF] Schwartz | Genetic monitoring as a promising tool for conservation and management[END_REF]. Consistently, correlation between estimates of local genetic differentiation from LOCALDIFF and genetic diversity were significant for allelic richness and not for heterozygosity, suggesting T. marmoratus populations are declining to a greater degree in the southern side of the study area. This supports the hypothesis that landscape genetic effects found in the present study reflect recent demographic changes, predominantly in the southern part of the study area, where land conversion for intensive arable crops production induced substantial T. marmoratus aquatic and terrestrial habitat loss (we estimated a 53 % decrease in hedgerows since 1950). Demographic inferences conducted with MIGRAINE were consistent with this hypothesis and indicated contrasted values between the two parts of the study area, in line with previous results. First, Nb, the neighborhood size, was significantly lower in the south, with Nb = 68.77 and 557.5 in the southern and northern parts of the study area, respectively. While Nb depends on both effective population size and dispersal [START_REF] Neel | Estimation of effective population size in continuously distributed populations: there goes the neighborhood[END_REF], MIGRAINE analysis also allowed us to infer area-specific effective population density and dispersal separately. Although not significantly different, MIGRAINE estimate of area-specific effective population density was almost 35% lower in the south, with 2.56 and 3.86 individuals per km² in the southern and northern parts of the study area, respectively. This contrast is consistent with the positive effect of local pond density on genetic diversity, since pond density was five times lower in the southern part of the study area. Overall, our effective population density estimates might seem low. However, in amphibians, the effective population size (Ne) is usually much lower than the census population size (N) (see Table 4 in [START_REF] Schmeller | Demographic and genetic estimates of effective population and breeding size in the Amphibian Rana temporaria[END_REF]. Indeed, T. marmoratus Ne/N ratio was estimated to range from 0.05 to 0.65 in five ponds of Western France [START_REF] Jehle | Contemporary gene flow and the spatio-temporal genetic structure of subdivided newt populations (Triturus cristatus, T. marmoratus)[END_REF]. In the same study, the effective population size was estimated in the order of 100-200 individuals in a 26.25 km² study site (~ 3.8 to 7.6 individuals per km²), thus consistent with our effective population density estimate in the northern part of the study area.

Interestingly, landscape structure in the northern part of the study area and in the Jehle et al (2005) study area were similar. Finally, area-specific MIGRAINE estimates of dispersal rate m and of the parameter of the dispersal function g also showed contrasted situation between the two parts of the study area. Dispersal rate was 46 % lower in the south, and the significant difference of area-specific g estimates indicated shorter dispersal range in this less connected part of the study area, consistently with a previous radiotracking study [START_REF] Trochet | Postbreeding movements in Marbled newts ( Caudata , Salamandridae ): A comparative radiotracking study in two habitat types[END_REF].

Our landscape scale estimates of dispersal rates and distances might appear much higher than was previously thought (Jehle et al. 2005). These authors concluded that dispersers could not travel distances > 1km but their estimates of pond specific immigration rate were only based on three closely located ponds as potential sources. Our analysis based on the IBD model and the sampling of more than 39 ponds over a large area appear more robust to infer landscapescale dispersal distances and rates. In their review paper on amphibians dispersal, [START_REF] Smith | Dispersal and the metapopulation in amphibian and paradigm ecology are all amphibian conservation: populations metapopulations?[END_REF] concluded that for salamanders ponds may receive migrating individuals from distances up to 8-9 km. Although T. marmoratus was not included in this review, our results suggest some individuals might occasionally travel similar distances.

Conclusion and recommendation for conservation

Our study contributes to a growing body of literature suggesting that agricultural intensification is harmful for pond breeding amphibians [START_REF] Crawford | Altered functional connectivity and genetic diversity of a threatened salamander in an agroecosystem[END_REF][START_REF] Curado | Amphibian pond loss as a function of landscape change -A case study over three decades in an agricultural area of northern France[END_REF][START_REF] Joly | Habitat matrix effects on pond occupancy in newts[END_REF][START_REF] Marty | Directionality of pre-and postbreeding migrations of a marbled newt population (Triturus marmoratus): Implications for buffer zone management[END_REF][START_REF] Boissinot | Amphibian diversity in farmlands: Combined influences of breeding-site and landscape attributes in western France[END_REF].

Combining demographic, movement, and genetic data is needed to fully understand spatial population dynamics for conservation [START_REF] Cayuela | Demographic and genetic approaches to study dispersal in wild animal populations : A methodological review[END_REF][START_REF] Wood | Combining genetic and demographic monitoring better informs conservation of an endangered urban snake[END_REF]. We demonstrated the necessity to move from site-specific to landscape level analyses to understand the population dynamics of T. marmoratus. Our results underline the need to base conservation planning at the landscape level [START_REF] Cushman | Effects of habitat loss and fragmentation on amphibians: A review and prospectus[END_REF]. In particular, increasing connectivity among populations appears to be a major issue for T. marmoratus and likely other amphibians in agricultural landscapes. The ultimate goal of conservation for amphibians should be long-term regional persistence by addressing issues at both local (notably quality of breeding-site) and landscape scale [START_REF] Boissinot | Amphibian diversity in farmlands: Combined influences of breeding-site and landscape attributes in western France[END_REF][START_REF] Semlitsch | Differentiating migration and dispersal processes for pond-breeding amphibians[END_REF]. Measures of local genetic differentiation correspond to F ST /(1-F ST ) calculated between the sampled ponds and 4 fictive neighboring populations, not shown, located at 2 km. Warmer colors indicates higher local genetic differentiation, i.e. higher genetic distances between sampled ponds and interpolated fictive neighbors.

Tables and figures

Figures 1 .

 1 Figures 1. Map of sampling area in the Deux-Sèvres department with (a) rivers network and ponds and (b) main land use categories according to CORINE Land Cover 2012 and BDTOPO 2015. Dots indicate sampling locations with red and blue colors referring to the population dataset and the northern and southern parts of the study area, respectively.

Figure 2 .

 2 Figure 2. Relation between Allelic richness (A r ) and pond density in a 1000 m buffer around the 39 ponds. Black dots and gray dots represent ponds from the northern and southern sides of the study area, respectively.

Figure 3 .

 3 Figure 3. Correlation between pairwise genetic differentiation (F ST /(1 -F ST )) and logarithm of the geographic distance between pairs of ponds. Filled dots indicate pairs of ponds that were significantly different (GENEPOP exact probability test for population differentiation).

Figure 4 .

 4 Figure 4. Map of local genetic differentiation inferred by the Bayesian kriging method implemented in LOCALDIFF . Shape represent the minimum convex polygon that encompasses the 39 ponds.

900 901 Table 1 :

 9011 Results of model averaging on models of allelic richness with delta AIC < 4 for the spatial scale with the highest R² (i. e. 1000 m buffer size). Estimates and P-values are presented for each landscape descriptor tested. Significant relations are in bold.

	Landscape variable	Estimate	Std. Error	z-value	P-value
	Arable crops	-0.065	0.037	1.750	0.080
	Ponds	0.107	0.035	3.021	0.003
	Woodlands	-0.026	0.060	0.441	0.659
	Rivers	0.021	0.034	0.604	0.546
	Buildings	0.020	0.036	0.554	0.579
	Grasslands	0.027	0.052	0.511	0.610
	Roads	-0.008	0.039	0.197	0.844

Table 3 :

 3 Inferences on demographic parameters by the software MIGRAINE. Point estimate values with data range outputted in brackets are shown. Scaled local population size (2Nµ), scaled emigration rate (2Nm), parameter of the geometric distribution of dispersal (g) and neighborhood size (Nb) are indicated. De, the effective density per 3250 x 3250 m bin and m, the migration rate among bins were estimated using 10 -4 [5 x 10 -4 -5 x 10 -5 ] as value and range for the mutation rate.

	Parameters	South	North
	2Nµ	0.011 [0.008 -0.014]	0.016 [0.012 -0.022]
	2Nm	25.4 [18.68 -38.78]	70.85 [47.12 -141.5]
	g	0.15 [0.003 -0.31]	0.41 [0.031 -NA]
	Nb	68.77 [48.34 -107.2]	557.5 [206.8 -22134]
	De	27 [5.4 -54]	40.8 [8.15 -81.5]
	m	0.235 [0.118 -1]	0.435 [0.217 -1]

Acknowledgements

This work was financially supported by the Conseil Régional du Poitou-Charentes (Projet "Bocage et Biodiversité"), the Fundation Lisea (Projet "Bocage") and the CNRS. This publication has been written with the support of the AgreenSkills+ fellowship programme, which has received funding from the EU's Seventh Framework Programme under grant agreement N° FP7-609398 (AgreenSkills+ contract). We address special thanks to people who accompanied us on the field. More specifically, we wish to thank Leïla Migault, Jérémie Souchet, Fabrice Connort, Jérôme Lallemand, Florian Doré, Alexandre Baduel, Nathan Braconnier, Jean-Pierre Lassaire, Stéphane Bonnet, Antoine Pasquier, Sabrina Lefebvre, Frédéric Migault, Morgan Pourrageau, Yannick Besson, Alain Texier, Clément Gouraud, Sébastien Billaud, Sabrina Maiano, Cédric Auburtin, Justine Vidal, Virginie Lassaire, Loïc Nau, Anthony Leguen and Eddy Leguen. We are also grateful to the local farmers in the study area for allowing us to study the ponds on their land. The study was carried out in compliance with all applicable animal care guidelines, and all required French permits were obtained (from the Préfecture des Deux-Sèvres).

Data Accessibility

The data sets with T. marmoratus sampling locations and microsatellite genotypes can be found on Dryad -https://doi.org/10.5061/dryad.mkkwh710s.

2) Genetic markers and diversity

The number of alleles per locus ranged from four to 16 in the global dataset (Table S1.1, Supplementary Information), with an average value of 8.5. Expected heterozygosity for each locus ranged from 0.047 to 0.895, with an average value of 0.576. Observed heterozygosity ranged from 0.046 to 0.813, with an average value of 0.486. There was a marginally significant overall heterozygosity deficit in the global dataset (P-value = 0.056, F IS = 0.127) that could be related to the presence of spatial structure ("Wahlund effect"). Indeed, only one of the 39 ponds with at least 7 individuals showed significant departure from HWE (Table S1.2, Supplementary Information). After FDR correction, there was no pair of loci showing significant LD in the global dataset.

3) Effect of landscape composition on ponds genetic diversity

Within pond allelic richness corrected for sample size (A r ) ranged from 2.684 to 3.773 (Table S1.2, Supplementary Information). R² comparison among linear mixed models based on landscape descriptors calculated in 500, 1000 and 2000 m buffers indicated 1000 m was the best spatial scale to explain A r (marginal R² = 0.005, 0.007 and 0.005 for 500, 1000 and 2000 m buffers, respectively). Several plausible best models (i.e. with delta AICc < 4) were identified at 1000 m (Table S2.4, Supplementary Information). Model averaging indicated that pond density in the surrounding landscape had a significant and positive effect on within pond A r while the negative effect of arable crop land cover was only marginal (Table 1 andFigure 2). Consistently, models at 500 m and 2000 m also reported a significant and positive effect of the amount of ponds in the surrounding landscape.

Regarding H o , R² comparison among linear mixed mode indicated 500 m was the best spatial scale (marginal R² = 0.006, 0.001 and 0.004 for 500, 1000 and 2000 m buffers, respectively).

Author Contributions


Ngoc-Son Nguyen

Hélène Magoariec
email: helene.magoariec@ec-lyon.fr

Eric Vincens

Bernard Cambou

On definition of a relevant meso-scale for upscaling the mechanical behavior of 3D granular materials

Keywords: 3D granular materials, DEM, Mechanical behavior, Upscaling, Meso-scale, Merging criteria

In this paper, the difficulty to adapt the definition of a meso-scale for 2D granular materials to derive the macroscopic mechanical behavior of 3D granular materials by upscaling is discussed. A starting point for a definition of a 3D meso-scale may be found in previous works where the transport of small particles through a set of defined pores and constrictions has been investigated. In this framework, different methods to derive a set of pores on which can be defined 3D meso-domains have been compared. Some clues are given about their relevancy to address the macroscopic mechanical behavior of 3D granular materials using an uspcaling from a meso-scale.

Introduction

The macroscopic mechanical behavior of 2D granular materials has been derived from a meso-scale using an uspcaling approach [START_REF] Satake | A discrete-mechanical approach to granular materials[END_REF][START_REF] Kuhn | Deformation measures for granular materials[END_REF][START_REF] Nguyen | Analysis of structure and strain at the meso-scale in 2D granular materials[END_REF][START_REF] Nguyen | Local stress analysis in granular materials at a mesoscale[END_REF][START_REF] Cambou | Granular Materials at Meso-scale: Towards a Change of Scale Approach[END_REF][START_REF] Nguyen | Towards a new approach for modeling the behavior of granular materials: A mesoscopic-macroscopic change of scale[END_REF]. Such a 2D meso-scale can be defined at the level of closed loops composed of pairs of contacting particles. Nguyen et al. [START_REF] Nguyen | Analysis of structure and strain at the meso-scale in 2D granular materials[END_REF][START_REF] Nguyen | Local stress analysis in granular materials at a mesoscale[END_REF] presented a method to construct these closed loops for 2D granular assemblies of disks based on Delaunay triangulation. First, this method consists in paving the studied domain into triangles, each of which joining the centers of three neighboring particles. In this Delaunay network, a branch is considered as closed if it joins two contacting particles or as open otherwise as illustrated in Figure 1. Every couple of adjacent triangles that share an open branch are then merged together. In doing so, each set of merged triangles forms a closed loop denoted mesodomain.

The local fabric of each meso-domain is described by a loop tensor proposed by Satake et al. [START_REF] Satake | A discrete-mechanical approach to granular materials[END_REF], from which the elongation direction and degree of each mesodomain are derived [START_REF] Nguyen | Analysis of structure and strain at the meso-scale in 2D granular materials[END_REF]. The strain and stress tensors associated to each meso-domain are defined from the particle movements and contact forces, respectively [START_REF] Kuhn | Deformation measures for granular materials[END_REF][START_REF] Nguyen | Analysis of structure and strain at the meso-scale in 2D granular materials[END_REF][START_REF] Nguyen | Local stress analysis in granular materials at a mesoscale[END_REF]. Nguyen et al. [START_REF] Nguyen | Towards a new approach for modeling the behavior of granular materials: A mesoscopic-macroscopic change of scale[END_REF] validated afterwards this new upscaling approach to obtain the macroscopic mechanical behavior of 2D granular materials from that of the meso-domains.

Such a meso-scale holds a definite physical meaning since important information such as the topology of the contact network and the arrangement of particles creating the meso-domains is de facto taken into account. However, the definition of this meso-scale for 2D materials cannot be straightforwardly extended to the 3D case since a 3D granular material has a much more complex micro-structure. Indeed, while each meso-domain for a 2D granular material clearly defines a closed pore which is not interconnected to its neighboring pores, the void space in a 3D granular material is characterized by a set of larger voids interconnected by smaller voids also denoted constrictions. A question that arises here is whether there exists a relevant meso-scale for 3D granular materials at which the features observed for 2D materials can be recovered. Moreover, are only 3D meso-domains composed of particles in mutual contact suited for a mechanical upscaling? Therefore, the definition of a relevant meso-scale for 3D granular materials seems a challenging and a complex task.

This paper presents some methods to define a 3D meso-scale, which can be envisioned as potential answers to the aforementioned questions. These methods come from another research domain where a partition of the void space of a granular material into pores and constrictions was provided. In particular, it was used to study the phenomena of transport and blockage of flowing smaller particles in a granular filter.

A comparison of the obtained sets of 3D meso-domains by these different methods is provided. This study mainly aims to pave the way to define a meso-scale for 3D granular materials. The validation of a given set of mesodomains as a relevant meso-scale in an upscaling approach for studying the mechanical behavior is beyond the scope of this work.

Weighted Delaunay triangulation

The 3D weighted Delaunay triangulation [START_REF] Aurenhammer | Power diagrams: Properties, algorithms and applications[END_REF] provides a partition of a packing of spheres into tetrahedra whose vertices correspond to the centers of four neighboring particles. In this partition, the weight associated to each sphere is equal to the square of its radius. However, the main drawback of this tessellation technique is that it tends to subdivide single pores into several smaller pores and it results in a uniform pore network for which each pore is systematically associated to four exits (constrictions) located on the four faces of each tetrahedron [START_REF] Al-Raoush | Comparison of network generation techniques for unconsolidated porous media[END_REF][START_REF] Reboul | A stastistical analysis of void size distribution in a simulated narrowly graded packing of spheres[END_REF]. Therefore, neighboring tetrahedra need to be merged according to a given criterion to obtain an accurate representation of the void space of a granular assembly [START_REF] Al-Raoush | Comparison of network generation techniques for unconsolidated porous media[END_REF][START_REF] Reboul | A stastistical analysis of void size distribution in a simulated narrowly graded packing of spheres[END_REF][START_REF] Nguyen | Prise en compte d'une échelle mésoscopique dans l'étude du comportement des milieux granulaires[END_REF][START_REF] Sufian | Pore shapes, volume distribution and orientations in monodisperse granular assemblies[END_REF][START_REF] Shire | Constriction size distributions of granular filters: a numerical study[END_REF][START_REF] Seblany | Merging criteria for defining pores and constrictions in numerical packing of spheres[END_REF][START_REF] Van Der Linden | A computational geometry approach to pore network construction for granular packings[END_REF]. In the following, we introduce a concept of sub-domains which are defined as polyhedral volumes composed of one tetrahedron or several merged tetrahedra. Each sub-domain is a candidate for being a meso-domain once the merging process finishes. Different merging criteria have been proposed and they can be classified into three main categories: (i) criteria based on contacts between particles, (ii) criteria based on a possible overlapping between the inscribed void spheres of two adjacent tetrahedra, and (iii) criteria based on the porosity of a face and those of the two adjacent sub-domains that share this common face.

Criteria based on contacts

By drawing an analogy between the 2D and 3D Delaunay networks, one might translate the concept of closed branches for 2D case to 3D case by considering that a triangular face shared by two adjacent tetrahedra is a closed loop if it connects three neighboring particles in mutual contact. In other words, a triangular face is open if it connects, at least, one couple of particles which are not in contact. Two adjacent tetrahedra that share an open face are then merged together. With such a merging criterion, we found that may be generated a very big meso-domain that occupies more than 70% of the total volume of the sample. Consequently, the merging criterion used for 2D case cannot be straightforwardly extended to 3D case.

O'Sullivan et al. [START_REF] O'sullivan | Contact based void partitioning to assess filtration properties in DEM simulations[END_REF] relaxed the above concept of closed loops for 3D materials. According to this method, a 3D closed loop is composed of either one triangular face with three contacts or two adjacent triangular faces with four contacts. Two adjacent pores must share a closed loop, otherwise they are merged together. This concept is appealing; however, the feasibility to define the meso-domains needs to be further investigated. This method is not studied here and the reader can refer to [START_REF] O'sullivan | Contact based void partitioning to assess filtration properties in DEM simulations[END_REF] for more details.

Instead of requiring that the particles of a mesodomain must form closed loops of three or four contacts, we require in the following merging criteria that the particles of a meso-domain are sufficiently close to each other, i.e. the space between them is sufficiently small, so that they can form an appropriate arrangement to resist shear loading.

Criteria based on overlapping inscribed void spheres

These criteria first consist in defining an inscribed void sphere for each tetrahedron and then in merging adjacent tetrahedra depending on the possible overlapping between their inscribed void spheres. Some criteria merge every two adjacent tetrahedra if there exists an overlapping between their inscribed void spheres [START_REF] Reboul | A stastistical analysis of void size distribution in a simulated narrowly graded packing of spheres[END_REF][START_REF] Sufian | Pore shapes, volume distribution and orientations in monodisperse granular assemblies[END_REF][START_REF] Seblany | Merging criteria for defining pores and constrictions in numerical packing of spheres[END_REF]. Other relaxed criteria make use of a threshold value directly related to the intensity of overlapping [START_REF] Al-Raoush | Comparison of network generation techniques for unconsolidated porous media[END_REF][START_REF] Shire | Constriction size distributions of granular filters: a numerical study[END_REF]. For instance, Al-Raoush et al. [START_REF] Al-Raoush | Comparison of network generation techniques for unconsolidated porous media[END_REF] proposed to merge two adjacent tetrahedra if the center of one of the two associated inscribed void sphere is located within the other one.

Reboul et al. [START_REF] Reboul | A stastistical analysis of void size distribution in a simulated narrowly graded packing of spheres[END_REF] adopted a strict criterion to merge tetrahedra in the primary partition given by the weighted Delaunay triangulation (called L 0 ). To avoid over-merging of tetrahedra, the authors introduced a merging level L 1 , for which only the tetrahedron with the largest inscribed void sphere and its direct neighbors are merged at each merging step. Seblany et al. [START_REF] Seblany | Merging criteria for defining pores and constrictions in numerical packing of spheres[END_REF] subsequently modified this merging criterion by introducing a merging level L 0 to eliminate strongly flat tetrahedra issued from the weighted Delaunay triangulation. This level consists in merging every couple of adjacent tetrahedra if the overlapping of the corresponding inscribed void spheres is greater than 99.9999% or if the constriction size is greater than the size of the corresponding inscribed void spheres. By doing so, about 10% of the couples of adjacent tetrahedra are merged whatever the studied packing of spheres (coefficient of uniformity ranging between 1.5 and 6) and density state (loosest or densest state).

The merging criterion L 1 implies that a maximum of five tetrahedra can be merged to form a single mesodomain. It was found that about 50% and 60% of the meso-domains in respective loose and dense assemblies of spheres having a coefficient of uniformity ranging between 1.5 and 6 are represented by a single tetrahedron [START_REF] Reboul | A stastistical analysis of void size distribution in a simulated narrowly graded packing of spheres[END_REF][START_REF] Seblany | Merging criteria for defining pores and constrictions in numerical packing of spheres[END_REF]. It was also proposed to extent this merging to the neighbors of the adjacent tetrahedra (L 2 criterion) but this criterion did not bring a much more complex set of pores [START_REF] Reboul | A stastistical analysis of void size distribution in a simulated narrowly graded packing of spheres[END_REF][START_REF] Seblany | Merging criteria for defining pores and constrictions in numerical packing of spheres[END_REF].

Criteria based on porosities of local particle arrangements

In this approach, a sub-domain is a polyhedral volume composed of one tetrahedron or several adjacent tetrahedra. For each sub-domain, the void space is approximated by an inscribed void sphere. The inscribed void sphere must not intersect adjacent solid particles and its center must be located inside the sub-domain. Starting from the weighted Delaunay triangulation, we proceed iteratively to identify all the sub-domains which do not fulfill this physical requirement and then merge them to their neighbors to form new sub-domains. Thus, this procedure is different from those used by Reboul et al. and Seblany et al. to define their primary set of subdomains (respectively L 0 and L 0 ' levels). It allows us to eliminate efficiently strongly flat tetrahedra issued from the Delaunay triangulation. After this primary merging procedure (PM), further merging is processed by using the two following criteria based on the facial and spacial porosities [START_REF] Nguyen | Prise en compte d'une échelle mésoscopique dans l'étude du comportement des milieux granulaires[END_REF].

Merging criterion based on the facial porosity

Two adjacent sub-domains share a face which can be composed of one or several triangular faces. The porosity n f of a face is defined as n f = A v /A where A v and A are the total void area and the total area of the face. We consider that a face shared by a couple of adjacent sub-domains is open if its porosity n f is bigger than a threshold value n th f . Every two adjacent sub-domains which share an open face are merged together. A similar merging criterion was used by Van der Linden et al. [START_REF] Van Der Linden | A computational geometry approach to pore network construction for granular packings[END_REF].

Merging criterion based on the relative facial porosity

The above criterion to define if a face shared by two adjacent sub-domains is open or closed might be relative: it might be considered to be closed if it is shared by two very porous sub-domains but it might be open if it is shared by two compact sub-domains. In this regard, we define a relative facial porosity n f as n f = n f / min(n 1 , n 2 ) where n 1 and n 2 are the porosities of the two adjacent sub-domains. A face is considered to be open if the relative facial porosity is bigger than a threshold value n th f .

Comparison between different criteria

The different merging techniques presented above are used to partition the loose numerical sample composed of 2500 spherical particles that was simulated with the Discrete Element Method and presented by Seblany et al. [START_REF] Seblany | Merging criteria for defining pores and constrictions in numerical packing of spheres[END_REF]. This sample with periodic boundary conditions has a porosity of 0.38 and a narrowly graded grain size distribution. The criteria based on the absolute and relative facial porosities, respectively abbreviated as FP and RFP, are used after the primary merging procedure (abbreviated as PM) as described in Section 5. These criteria are compared to the L 1 criterion based on overlapping inscribed void spheres (Section 4). The pore size of each sub-domain is characterized by the diameter of an equivalent void sphere whose volume is equal to the void volume of the sub-domain under consideration. Figures 2 and3 show this comparison where different threshold values are used for the FP and RFP criteria. In both figures, the pore size distribution of the tetrahedra obtained from the weighted Delaunay triangulation is represented by the curve L 0 . It can be seen in Figure 2 that the pore size distribution after the primary merging procedure (in Figure 2, the PM

Probability density

Equivalent sphere diameter (mm)

L 1 L 0 PM FP, n th f = 0.26 FP, n th f = 0.3 FP, n th f = 0.4
Fig. 2 Probability density versus the equivalent sphere diameter obtained with the L 0 , L 1 , PM and FP criteria.

curve is almost coincident with the red curve) is significantly lowered with larger pore sizes, compared to the L 0 curve. This means that the PM criterion merges efficiently the flat tetrahedra, which have medium or large pore volumes, giving rise to sub-domains with bigger pore volumes. The FP criterion with a threshold value n th f = 0.4 does not notably change the pore size distribution. This means that the PM criterion allows us to eliminate almost highly porous faces (n f > 0.4). A notable change in the pore size distribution can be observed for n th f = 0.3 that leads to some meso-domains with more than 20 tetrahedra. Van der Linden et al. [START_REF] Van Der Linden | A computational geometry approach to pore network construction for granular packings[END_REF] recommended a threshold value n th f > 0.3 since a lower threshold value tends to over-merge tetrahedra. We can also remark that the pore size distribution given by the FP criterion with n th f = 0.3 is far from that obtained with the L 1 criterion. This is due to the fact that only 14% of the faces meet the FP criterion with n th f = 0.3, much lower than 30% of the faces that meet the L 1 criterion. One can note that even though the threshold value n th f = 0.26 is chosen such that the number of open faces according to the FP criterion is equal to that given by the L 1 criterion, these criteria give different results. Finally, the FP criterion slightly tends to merge more sub-domains with medium or big void volumes to form new sub-domains with bigger void volumes than the L 1 criterion, while the sub-domains obtained with the L 1 criterion are limited in size. In addition, the sub-domains with low porosities are not affected by the FP criterion as they share most likely faces with low porosities.

Figure 3 shows that the RFP criterion with a threshold value n th f = 1.0 does not result in a significant change in the pore size distribution compared to the PM curve. This means that almost all sub-domains obtained after the primary merging technique are less porous than their faces. As the threshold value n th f is reduced, the pore size distribution changes greatly and gets closer to the one obtained with the L 1 criterion. For n th f = 0.75, the pore size distribution curve is quite close to the L 1 curve. We found that about 24% of the faces meet the RFP criterion with n th f = 0.75, which is quite close to the number of 30% given by the L 1 criterion, but this is not the only reason to explain the observed result. We have shown previously that, with almost the same number of open faces, the FP criterion with n th f = 0.26 gives quite a different result from that given by the L 1 criterion. Among the closed faces defined by the FP criterion, some of them might be open according to the RFP criterion as they might be shared by adjacent sub-domains of low porosities (their relative porosity becomes bigger than the absolute one) or vice versa. Indeed, 6% of the closed faces defined by the FP criterion with n th f = 0.26 become open according to the RFP criterion with n th f = 0.75. On the other hand, 24% of the open faces defined by the former become closed according to the latter as they are shared by two sub-domains with high porosities. As a consequence, the RFP criterion merges more sub-domains with small void volumes but less sub-domains with medium or big void volumes than the FP criterion. Finally, in the context of the L 1 merging criterion, it was noted that the overlapping of inscribed void spheres is intrinsically very important in packing of spheres. It means that if two adjacent inscribed void spheres overlap, the opening size of the involved constriction is rather large compared to the pore size of the two corresponding adjacent tetrahedra, leading generally to a high value of the relative porosity n f . It may explain why the RFP criterion with n th f = 0.75 gives a pore size distribution close to that given by the L 1 criterion. By comparing Figures 2 and3, it can be concluded that the RFP criterion provides a better partition of a granular volume into meso-domains than the FP criterion.

Discussion

The L 1 and RFP criteria provide a complete partition of a given 3D granular material into meso-domains. In particular, strongly flat tetrahedra that do not hold any physical meaning are removed from the weighted Delaunay triangulation by using the L 0 level in the L 1 criterion or the primary merging procedure in the RFP criterion. Both criteria give the meso-domains of different sizes and shapes, which may reflect the inherent heterogeneity of the internal state of granular materials. Figure 4 shows that the local porosity of the mesodomains obtained with the L 1 and RFP criteria versus their number of contacts. In this figure, a mean porosity is calculated over all the meso-domains having the same number of contacts. This figure shows that there is a strong correlation between the local porosity and the number of contacts within the meso-domains for both criteria: meso-domains with higher number of contacts tend to occupy zones with higher porosity. This correlation was also observed for 2D closed loops [START_REF] Nguyen | Analysis of structure and strain at the meso-scale in 2D granular materials[END_REF]. This preliminary result suggests that although both criteria are based on the void space between neighboring particles, mechanical properties of the meso-domains related to contacts between particles (those observed for 2D closed loops) might be recovered for 3D materials.

For these reasons, the meso-domains defined by either the L 1 or RFP criterion hold some basic features for being a possible candidate to develop an upscaling approach for the mechanical behavior of 3D granular materials.

Conclusion

In this paper, we proposed two methods to extend the concept of a meso-scale defined initially for 2D granular materials to 3D media. These methods are based on the weighted Delaunay triangulation where adjacent tetrahedra are merged by using criteria holding a physical meaning. Two specific criteria, which take into consideration the void space within local arrangements of particles, were investigated. It was shown that both led to similar partitions of the space into a set of meso-domains. They hold some basic features to be possible candidates for defining a meso-scale on which a mechanical upscaling of a 3D granular material can be applied. Further evidence must be given in this respect and a deeper study will be carried in a future paper.

Compliance with ethical standards

Fig. 1

 1 Fig. 1 Part of a 2D granular assembly is subdivided into closed-loops. Solid and dashed lines depict a closed and an open branch, respectively.

Fig. 3

 3 Fig. 3 Probability density versus the equivalent sphere diameter obtained with the L 0 , L 1 , PM and RFP criteria.

Fig. 4

 4 Fig.4Correlation between the local porosity and the number of contacts for the meso-domains obtained with the L 1 and RFP criteria with n th f = 0.75.

Conflict of Interest: The authors declare that they have no conflict of interest