Madhukar Bhat
email: madhukar.bhat@vitec.com

Jean-Marc Thiesse
email: jean-marc.thiesse@vitec.com

Patrick Le Callet
email: patrick.lecallet@univ-nantes.fr

HVS based perceptual pre-processing for video coding

Keywords: Pre-processing, Contrast sensitivity function, Subjective test, Human Visual system, Visual masking

This paper proposes a perceptually optimized preprocessing technique using state of the art Human Visual System (HVS) model suitable for video compression to reduce bitrate at the same quality. Visual masking models to accurately account HVS has been considered. Frequencies which are visually indistinguishable and need not be encoded are removed within visibility threshold. The scheme is optimized for multiple viewing distances to consider real-world scenarios. Extensive subjective and objective evaluation has been conducted to evaluate proposed pre-processing. Investigation shows significant bit-rate savings compared to a professional real time HEVC video encoder.

I. INTRODUCTION

Perceptual optimization to remove perceptual redundancy and adaptation to real-world scenarios has been an active area in video compression [START_REF] Chen | Perceptual video coding: Challenges and approaches[END_REF]. Perceptual optimizations include modifications in encoding the video according to the scene presented and visual changes to the video outside the encoding loop [START_REF] Chen | Perceptual video coding: Challenges and approaches[END_REF]. These adaptations aim to introduce psycho-visual bitrate savings and minimize resource allocation for the encoding process.

The coding loop is the most critical part of nowadays realtime encoders due to the amount of computations involved in coding decisions and reconstruction together. Thus, it is highly desirable to implement processing outside this loop using the source video to ease encoding. Using pipelining of the various modules finally allow designing the full encoder. Pre-processing is one relevant way to achieve perceptual optimization outside the video coding loop. It aims to remove visually redundant noise and high frequencies from the source video to ease compression and coding efficiency at a given quality level of the video.

Pre-processing can be carried out both in pixel domain or in frequency domain in many ways. There have been attempts of pre-processing as denoising filter [START_REF] Karunaratne | A rate-distortion optimal video preprocessing algorithm[END_REF] whose strengths were controlled by encoding parameters such as motion vectors and energy of the residual [START_REF] Lee | Automatic prefilter control by video encoder statistics[END_REF] [START_REF] Jain | A low-complexity, motion-robust, spatiotemporally adaptive video de-noiser with inloop noise estimation[END_REF]. These works mainly to reduce noisy components in the video introducing unnecessary complexity during encoding. In [START_REF] Shaw | Color difference weighted adaptive residual preprocessing using perceptual modeling for video compression[END_REF] [START_REF] Shao-Ping | Saliency-based fidelity adaptation preprocessing for video coding[END_REF] pre-processing is seen as an adaptive low-pass filter to reduce high-frequency content which are visually not important but potentially challenging to encode. The adaptive low-pass filtering is applied using some more methods like diffusion [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], bilateral filtering [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] and adaptive weighted averaging [START_REF] Ozkan | Adaptive motion-compensated filtering of noisy image sequences[END_REF]. Vidal et al. [START_REF] Vidal | New adaptive filters as perceptual preprocessing for rate-quality performance optimization of video coding[END_REF] introduced two filters combining good features of both bilateral and adaptive weighting average filters. In [11] [12] an anisotropic filter was proposed based on contrast sensitivity estimation for different real world video delivery scenarios such as viewing angle and distance. Just Noticeable Difference models(JND) have been employed in [START_REF] Yang | Motion-compensated residue preprocessing in video coding based on just-noticeable-distortion profile[END_REF] [START_REF] Ding | Video pre-processing with JND-based Gaussian filtering of superpixels[END_REF] to control pre-filters. These JND models remove information within the noticeable threshold from a visual signal to avoid overestimation of removable frequencies. There is still room for improvements in accurate perceptual modeling of frequencies that could be removed combined with condition to be within the visibility threshold (JND) to ease the encoding process.

In contrast to above methods to model removable frequencies, Human Visual System (HVS) models are useful in accurately measuring distortions in content that are not visible to the human eye. In [START_REF] Barten | Evaluation of subjective image quality with the square-root integral method[END_REF] [16] [START_REF] Watson | A standard model for foveal detection of spatial contrast ModelFest experiment[END_REF] [18] computational models of HVS intended to provide an estimate for the visibility difference between a distorted and original signal have been proposed. They include following elements: non-linear sensitivity to calculate absolute luminance, contrast sensitivity to visual frequencies, oblique effect and sub-band decomposition into visual channels, masking effect, Psychometric function and error pooling. These HVS tools can be carefully used to remove or add information from the source video within the visibility threshold [START_REF] Urvoy | Perceptual DFT watermarking with improved detection and robustness to geometrical distortions[END_REF]. Hence, HVS models theoretically perform better than other solutions in accurately measuring redundant information that could be removed from a source within JND to perform pre-processing algorithm.

This paper introduces a new perceptual pre-processing filter named HVSPP which uses such HVS model for Luma and an adaptive low-pass filter for Chroma within JND threshold. HVSPP is built on visual tools such as contrast masking and luminance masking optimized for multiple viewing distance and screen brightness. The contrast sensitivity function (CSF) is calculated in the complex multi-scale steerable pyramid [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] domain of the visual channel model of physical luminance of the source video. The complete model of masking effect is normalized to suit the pre-processing within the JND threshold to get the pre-processed video. Coding efficiency in-terms of bit-rate after application of HVSPP is verified using extensive subjective and objective tests. The tests were conducted on a real time HEVC encoder developed by VITEC, a French company which expertises in real time video encoders.

The paper is organized as follows: Section II explains HVSPP algorithm in detail. Experimental setup for the subjective assessment of the proposed solution is carefully explained in section III. Corresponding results are presented in section IV together with objective evaluation. Conclusions are drawn in section V.

II. FRAMEWORK FOR PROPOSED PERCEPTUAL PRE-PROCESSING FOR VIDEO CODING

In this section, the detail of the proposed HVS based perceptual pre-processing (HVSPP) framework is given. Different HVS tools are applied in the pre-processing step to Luma component. Chroma component is applied with a filter based on morphological operations. The final HVSPP signal is then encoded instead of source. Next subsections describe in Luma and Chroma Pre-processing in detail.

A. Luma pre-processing using HVS tools

As mentioned in the introduction, Luma pre-processing in the proposed algorithm uses HVS models to remove redundant information. The accuracy of the HVS model to account for masking effects in-terms of luminance and contrast with optimization for different real world use cases. This framework considers two viewing distances at 3H (height of the display) which is the recommended minimum distance of the observer from full HD display and 4.5H which is the reference viewing distance for consumers for optimization and evaluation. The luminance of the targeted display is 250cd/m2.

The proposed work-flow of the HVSPP for Luma is shown in Fig 1 . The HVSPP algorithm is carried out in block-wise manner to target hardware implementation. As shown in the Fig. 1 at first, Luma block is converted to physical luminance value. Then, the visual channel model is calculated. This is followed by, the visual channel model decomposed by multiscale representation using the steerable pyramid. The contrast masking is then calculated for each scale and orientation. The calculated contrast masking at each scale and orientation is normalized and applied to multi-scale decomposition of luminance. Finally, the pre-processed block at each scale and orientation is reconstructed back to luminance and Luma. This Luma block is the HVSPP block.

B. Luminance calculation

Generally in video representation and storage gamma encoded values are used for display purpose. The contrast sensitivity functions used in this paper work on the physical luminance in cd/m2. In the proposed model YUV videos are

L(in cd/m2) = M * Y γ + black level (1)
γ value here is 2.2. Y is the normalized Luma component and M is the maximum brightness of the screen in cd/m2. black level is the minimum brightness of the screen.

C. Visual channel model

Visual channel model mimics the human visual system closely when luminance L is fallen on the eye [START_REF] Mantiuk | HDR-VDP-2: a calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. The model accounts for Intra-ocular light scatter which occurs when entering eye, photo receptor sensitivity which are the probability that photo receptors senses a photon, luminance masking which are responses of the eye for light which are usually non-linear and finally achromatic response which is joint cone and rod response. This transformation is represented as P.

D. Multi scale decomposition

The visual cortex of the human visual system presumably goes through decomposition of signals.To mimic that multiscale decomposition is applied to visual channel model P. CSF model are better adapted at localized frequency in steerable pyramid compared to pixel domain, wavelet and cortex transform [START_REF] Mantiuk | HDR-VDP-2: a calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. The multi-scale decomposition of the visual channel model P is achieved using a complex steerable pyramid in the current framework. To avoid reconstruction error introduced in the real-valued pyramid, a complex representation of the steerable pyramid is adapted in this framework. As shown in Fig. 1 there are two decompositions one is for P for calculating contrast masking and other is for L for applying contrast masking.

E. Contrast sensitivity calculation and application

In the proposed framework, the main step is to remove information from a Luma block within the noticeable threshold. Overall information that could be removed in each band is modeled as a sum of signal dependent noise (neural noise {nMask}) and signal independent noise (visual masking {CSF}). Signal independent noise which is CSF, try to remove information such as noise in the block within visibility threshold . The CSF N CSF (ρ) is modeled according to Barten's CSF model [START_REF] Barten | Contrast sensitivity of the human eye and its effects on image quality[END_REF],

N CSF (ρ) = 1 (1 + (p 1 ρ) p 2).(1 -e -(ρ/7) 2 ) p3 , ( 2
)
where ρ is the spatial frequency in cycles per degree and p 1 to p 3 are fitted parameters adapted to 250cd/m2. The signal dependent noise N nM ask try to remove information such as oscillations in texture and structure which are within visibility threshold [START_REF] Mantiuk | HDR-VDP-2: a calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. For f-th spatial frequency band and o-th orientation of the complex steerable pyramid, overall contrast masking(CM) is expressed as,

CM (f, o) = 1 N 2 nM ask + N 2p nCSF (3)
where, p is the gain control function which controls the shape of the masking function. The typical range of p is between 1.3 (contrast discrimination) up to 4 (contrast detection). In case of pre-processing the objective is to remove as much information as possible from the scene which are invisible i.e within JND.

To make sure contrast masking (CM) is within JND threshold at scale f and orientation o of steerable pyramid, psychometric function can be calculated. This gives probability of detection at that band. It is expressed as,

Ψ(f, o) = 1 -e log(0.5) |B P P P (f,o)-Borg (f,o)| CM (f,o) p (4)
where, B {PPPB} is the pre-processed band, B {org} is original band. The value of p is optimized in such a way that Ψ(f, o) should be less than 0.5. Another objective attached for p is coding efficiency after encoding. Using these two conditions p has been set to 1.5 for high-pass and 3 for bandpass after careful analysis. The total neural noise CM(f,o) is in the domain of multiscale decomposition of visual channel model(P). To have a reconstructible solution for application of neural noise to remove information from each band following normalization is performed.

Bt normalized (f, o) = CM (f, o) max(CM f, o) (5)
The information that can be removed from each band within probability of noticeable visual change is the normalized neural noise. This removal is performed as follows,

B P P P B (f, o) = B org (f, o) -[B org (f, o) * Bt normalized (f, o)] (6)
Then reconstruction of pre-processed steerable pyramid is carried out. Finally the physical luminance is converted back to Luma representation to give pre-processed block in the pixel domain. Fig. 2b shows the exaggerated absolute difference between source frame and pre-processed frame. HVSPP removes un-necessary noise in the flat region like sky along with visually redundant oscillations in texture and edges.

F. Chroma pre-processing

Pre-processing of Chroma follows a morphological analysis based pixel domain filter process controlled by quantization level of video to be encoded. According to encoding QP a threshold of how much a pixel can be changed is set. Then 3x3 and 5x5 Gaussian filters are applied to image separately and also gradient of the image with 3x3 and 5x5 kernel is found. Then a particular pixel is changed according to its gradient. If a pixel is not a gradient pixel with 5x5 gradient filter then pixel value of Gaussian filtered image with 5x5 kernel is applied.

Then the same step is repeated for 3x3 kernel. In this way strong edges are kept and flat areas are modified. 5x5 kernel is used see if the pixel is just a noise in flat area more robustly.

III. EXPERIMENTAL SETUP

The performance of the perceptual pre-processing framework is evaluated on four 1080p video sequences. Extensive paired comparison subjective tests and in addition objective tests have been conducted to validate performances of HVSPP in the context of video compression using real time HEVC encoder developed by VITEC. Subjective tests incorporated to evaluate the proposed framework follows ITU-T P.910 [START_REF]Subjective video quality assessment methods for multimedia applications[END_REF]. The subjective test protocols is explained in the following paragraphs.

A. Content and Stimuli selection

Video sequences used for evaluating perceptual preprocessing algorithm is selected containing a wide range of real-world scenarios. Four 1080p videos at 50fps is selected for the subjective test. Two videos are JCT-VC common test condition videos [START_REF] Jctvc-L | Common test conditions and software reference configurations[END_REF]. Another two are from SVT test videos [START_REF] Haglund | The SVT high definition multi format test set[END_REF]. Selected videos are listed in

B. Test environment

The testing environment incorporated here is based on Recommendation ITU-T P.910 [START_REF]Subjective video quality assessment methods for multimedia applications[END_REF]. Two viewing distances corresponding to the height of the display (H) at 3H and 4.5H are used. The luminance of the screen used for the test is 250cd/m2.

C. Paired comparison test design

Paired comparison test conducted in this work uses squared design [START_REF] Dykstra | Rank analysis of incomplete block designs: A method of paired comparisons employing unequal repetitions on pairs[END_REF] with 9 Hypothetical Reference Circuits (HRCs) and 18 pairs to compare for each sequence. HRC includes source and eight videos encoded at different QPs from HVSPP and original video. 30 naive observers took part in the test for each viewing distance.

IV. RESULTS

In this section, performance of HVSPP is discussed for the setup explained in the last section for video compression. First, BD-rate is calculated using objective metrics and then subjective test data is carefully analyzed to measure performance of HVSPP. Objective metrics used for the BD-rate measurement are PSNR, MS-SSIM, VMAF and HDR-VDP. TABLE II and TABLE III shows BD-Rate for HVSPP against original at 3H and 4.5H respectively. The performance of HVSPP is almost the same as original at 3H and worse at 4.5H for PSNR. However, this is expected as PSNR does not perform well to capture perceptual improvements [START_REF] Huynh-Thu | Scope of validity of PSNR in image/video quality assessment[END_REF]. Other metrics show BD-rate gain for HVSPP compared to original at the same quality. At 4.5H, average bit-rate savings at the same quality increase compared to 3H for all perceptual metrics except PSNR. Objective metrics except HDR-VDP are not optimized for different viewing distances. Hence, the more suitable way of accurately quantifying improvements of perceptual algorithms like HVSPP is subjective test [START_REF] Hanhart | How to benchmark objective quality metrics from paired comparison data?[END_REF]. As explained in section III, an extensive subjective test is conducted to measure performance of HVSPP. Paired comparison data obtained from subjective test is processed using Bradley-Terry model [START_REF] Bradley | 14 paired comparisons: Some basic procedures and examples[END_REF] to get relative score (D) for each sequence with 95% confidence intervals(CI). It is important to note for Bradley-Terry scores that, BD-quality at same rate for different sequences does not have any meaning [START_REF] Hanhart | How to benchmark objective quality metrics from paired comparison data?[END_REF]. This is because subjective scores are calculated relatively on persequence basis with 95% CI and are not comparable to other sequences. Hence, only BD-rate at relative subjective scores can be calculated using subjective comparison of encoder based on fitted curve [START_REF] Hanhart | Calculation of average coding efficiency based on subjective quality scores[END_REF]. This model takes 95% CI attached to relative scores from Bradley-Terry scores (D) into account. Mean, maximum and minimum bit-rate gain for both the distances are given in TABLE IV. They are calculated as,

∆R = ψ(r 1 (D), r2 (D), D L , D H ) (7)
∆R min = min{ψ(r - 1 (D), r+ 2 (D), D L , D H ), ψ(r 1 (D), r2 (D), D L , D H )} (8)
∆R max = max{ψ(r - 1 (D), r+ 2 (D), D L , D H ), ψ(r 1 (D), r2 (D), D L , D H )} (9)
Where, ∆R is the mean bit-rate difference between original and HVSPP using mean relative subjective scores using ψ fitting function [START_REF] Hanhart | Calculation of average coding efficiency based on subjective quality scores[END_REF]. r1 (D) indicates rate at particular relative scores using mean scores (mean grade series from D L to D H ). r+/- In most cases, the stimuli pair at the same QP are not significantly different, as shown in Fig. 3. This means, there is no significant subjective quality difference between HVSPP and original video at same QP for most of the cases with 95% of CI. This implies that, statistically bit-rate reduction produced by HVSPP at a particular QP can be considered as bit-rate saving at that QP for HVSPP. The average bit-rate reduction at QP's from Section III for different sequences are shown in Fig. 4. Bit-rate reduction with HVSPP compared to original reduces with an increase QP, because overall bit-rate needed to encode the video also reduces. It can be concluded that, subjectively HVSPP can produce a bit-rate gain of up-to 36% depending on QP with 95% CI.

V. CONCLUSION

This paper introduced perceptual pre-processing to improve video compression efficiency. In the proposed framework, various HVS tools are used to remove information from the source video within the JND threshold. The performance of HVSPP has been carefully evaluated using objective metrics and a subjective test. The analysis of results demonstrates significant benefits of using this perceptual pre-processing. HVSPP demonstrates very promising subjective BD-rate gains of 5.75% at 3H and 8.5% at 4.5H. This performance is captured by most of the tested objective metrics except PSNR. Above all, the subjective assessment proves that at same QP, the difference between HVSPP and original videos are not significant in most of the cases with HVSPP requiring significantly less bit-rate to encode.

Fig. 2 :

 2 Fig. 2: (a) For crowdrun video sequence, Absolute difference (original-HVSPP) for pre-processed frame(exaggerated by 100 times).

Fig. 3 :Fig. 4 :

 34 Fig. 3: Significance of difference between HVSPP and original encoded at same QP

TABLE I .

 I Selection of stimuli for subjective test is based on HVSPP and original video encoded at 4 QP's. TABLE I shows QP's selected for each sequence. The selection of QP is based on their RD-Curve.

	Sequences	Selected QP's
	Crowdrun	22,27,32,38
	BasketBallDrive	20,27,35,40
	BQTerrace	21,25,32,38
	OldtownCross	22,27,32,37

TABLE I :

 I Selected sequences and stimuli

TABLE II

 II

	: BD-rate at same quality for observer at 3H the
	height of screen
	Bit-rate saving/	PSNR	VMAF	MS-SSIM HDR-VDP
	/Sequence
	Crowdrun	1.87%	-2.52%	-0.27%	-1.19%
	Basketball	4.67%	-3.42%	-1.92%	-5.02%
	BQ terrace	2.27%	-5.25%	-2.13%	-5.07%
	oldtowncross	0.67%	-3.23%	-3.93%	-5.09%
	Overall	2.37%	-3.61%	-2.06%	-4.09%

TABLE III :

 III BD-rate at same quality for observer at 4.5H the height of screen

  indicates rate at relative score including confidence interval (D ± CI, max/min grade series). ∆R min and ∆R max indicates minimum and maximum bit-rate difference, hence indicating maximum and minimum bit-rate gain respectively[START_REF] Hanhart | Calculation of average coding efficiency based on subjective quality scores[END_REF]. Overall mean subjective bit-rate savings are 5.76% at 3H and 8.5% at 4.5H with maximum subjective bit-rate saving going up to 17.41% and 19.86% for respective distances.

			Mean	Minimum	Maximum
		Bit-rate gain	Bit-rate gain	Bit-rate gain
		(∆R)	(∆Rmax)	(∆R min )
	Sequences	3H	4.5H	3H	4.5H	3H	4.5H
	Crowdrun	-6.74%	-4.71%	4.94% 4.22%	-13.07% -15.64%
	Basketball	-5.02%	-9.01%	4.23% 4.41%	-15.64% -21.82%
	BQ terrace	-6.05%	-2.53%	14.5% 14.09% -20.59% -14.93%
	oldtowncross -5.24%	-17.74%	4.06% -5.00%	-20.32% -27.08%
	Overall	-5.76%	-8.5%	6.93% 4.43%	-17.41% -19.86%

1 (D)

TABLE IV

 IV

		: BD-rate at same quality for observer at 3H and
	4.5H. (A negative value indicates decrease in bit-rate (gain)
	for HVSPP at same subjective score compared to original)
		Pre-Processing Vs Reference
	75.00%	70.00%
					Significantly better
					No significant differrence
	15.00%	10.00%	20.00%	10.00%	Significantly worse
	3H			4.5H


Louise Fouqueau

Denis Roze
email: roze@sb-roscoff.fr

The evolution of sex along an environmental gradient

Keywords: Evolutionary quantitative genetics, geographic parthenogenesis, linkage disequilibrium, range margins, recombination

While temporally changing environments generally favor sex and recombination, the effects of spatial environmental heterogeneity have been less explored. In this paper, we use a classical model of adaptation along an environmental gradient to study the selective forces acting on reproductive mode evolution in the central and marginal parts of the distribution range of a species. The model considers a polygenic trait under stabilizing selection (the optimal trait value changing across space) and includes a demographic component imposing range limits. The results show that in the central part of the range (where populations are well adapted), recombination tends to increase the mean fitness of offspring in regimes where drift is sufficiently strong (generating a benefit for sex), while it has the opposite effect when the effect of drift stays negligible. However, these effects remain weak and are easily overwhelmed by slight intrinsic fitness differences between sexuals and asexuals. In agreement with previous results, asexuality may be favored in marginal populations, as it can preserve adaptation to extreme conditions. However, a substantial advantage of asexuality is possible only in conditions maintaining a strong maladaptation of sexuals at range limits (high effective environmental gradient, weak selection at loci coding for the trait).

INTRODUCTION

Despite the many costs associated with sexual reproduction, exclusive asexuality is rare in the eukaryotic kingdom, at least among multicellular plants and animals: the proportion of obligate asexuals is thought to be of the order of 0.1% in animals, and less than 1% in angiosperms [START_REF] Vrijenhoek | Animal clones and diversity[END_REF][START_REF] Whitton | The dynamic nature of apomixis in the angiosperms[END_REF]. Interestingly, asexual lineages often present different geographical distributions than their closest sexual relatives, a phenomenon described as "geographic parthenogenesis" by [START_REF] Vandel | La parthénogenèse géographique. Contribution à l'étude biologique et cytologique de la parthénogenèse naturelle[END_REF]. In a number of cases, asexuals are found at higher latitudes or altitudes, or in environments qualified as extreme or marginal [START_REF] Glesener | Sexuality and the components of environmental uncertainty: clues from geographic parthenogenesis in terrestrial animals[END_REF][START_REF] Lynch | Destabilizing hybridization, general purpose genotypes and geographic parthenogenesis[END_REF][START_REF] Bierzychudek | Patterns in plant parthenogenesis[END_REF][START_REF] Haag | A new hypothesis to explain geographic parthenogenesis[END_REF][START_REF] Hörandl | Geographical parthenogenesis: opportunities for asexuality[END_REF]. A large variety of patterns are observed, however, asexuals sometimes occupying a much wider geographic range than their sexual counterparts, possibly coexisting with them in part of the range (e.g., [START_REF] Van Dijk | Ecological and evolutionary opportunities of apomixis: insights from Taraxacum and Chondrilla[END_REF][START_REF] Schmit | Linking present environment and the segregation of reproductive modes (geographical parthenogenesis) in Eucypris virens (Crustacea: Ostracoda)[END_REF]. Several hypotheses have been proposed to explain these different distributions [START_REF] Haag | A new hypothesis to explain geographic parthenogenesis[END_REF][START_REF] Tilquin | What does the geography of parthenogenesis teach us about sex? Phil[END_REF]. Some are based on a direct selective advantage of uniparental reproduction, such as a better colonizing ability [START_REF] Baker | Self-compatibility and establishment after 'long-distance' dispersal[END_REF][START_REF] Cuellar | Biogeography of parthenogenetic animals[END_REF]. Other involve spatial variation in the indirect benefits associated with sexual recombination, due to environmental heterogeneity: in particular, sex may be less advantageous in areas facing high rates of maladapted gene flow [START_REF] Antonovics | Evolution in closely adjacent plant populations V. Evolution of self-fertility[END_REF][START_REF] Peck | Explaining the geographic distributions of sexual and asexual populations[END_REF], or where inbreeding is high due to frequent population bottlenecks [START_REF] Haag | A new hypothesis to explain geographic parthenogenesis[END_REF], or where parasites or predators are less abundant [START_REF] Glesener | Sexuality and the components of environmental uncertainty: clues from geographic parthenogenesis in terrestrial animals[END_REF]. Although these theories have hardly ever been tested empirically, the different distributions of sexual and asexual species represent an interesting avenue for exploring the possible benefits of sexual reproduction [START_REF] Bell | The Masterpiece of Nature[END_REF][START_REF] Tilquin | What does the geography of parthenogenesis teach us about sex? Phil[END_REF].

Most theories on the evolutionary advantage of sex and recombination assume spatially homogeneous environments [START_REF] Agrawal | Evolution of sex: why do organisms shuffle their genotypes?[END_REF][START_REF] Otto | The evolutionary enigma of sex[END_REF]. Several models showed that recombination may be favoured because it increases genetic variation and therefore the rate of response to directional selection, whenever negative genetic associations (linkage disequilibria) exist within populations (that is, when alleles increasing fitness tend to be associated with alleles decreasing fitness at other loci). Such negative associations may be generated by negative epistasis (a negative curvature of the fitness function, occurring in particular in most models of optimizing selection on quantitative traits, e.g., [START_REF] Smith | Selection for recombination in a polygenic model[END_REF], 1988;[START_REF] Charlesworth | Directional selection and the evolution of sex and recombination[END_REF][START_REF] Barton | A general model for the evolution of recombination[END_REF][START_REF] Kondrashov | Evolution of amphimixis and recombination under fluctuating selection in one and many traits[END_REF][START_REF] Vanhoenacker | Stabilizing selection, mutational bias, and the evolution of sex[END_REF]. In that case, recombination tends to decrease the mean fitness of offspring in the short term, although it increases the rate of adaptation in a longer term [START_REF] Barton | A general model for the evolution of recombination[END_REF][START_REF] Otto | The evolution of recombination in changing environments[END_REF][START_REF] Otto | Resolving the paradox of sex and recombination[END_REF]. Negative associations are also generated by stochastic effects occurring in finite populations [START_REF] Hill | The effect of linkage on limits to artificial selection[END_REF][START_REF] Felsenstein | The evolutionary advantage of recombination[END_REF], favoring sex and recombination in the absence of epistasis (e.g., [START_REF] Otto | The evolutionary enigma of sex[END_REF]Barton, 1997, 2001;[START_REF] Barton | Evolution of recombination due to random drift[END_REF][START_REF] Keightley | Interference among deleterious mutations favours sex and recombination in finite populations[END_REF][START_REF] Roze | Selection for sex in finite populations[END_REF].

In spatially heterogeneous environments, another possible source of linkage disequilibrium corresponds to the mixing between populations with different allele frequencies -for instance, due to spatial differences in selection [START_REF] Barton | Genetic analysis of hybrid zones[END_REF].

In particular, if alleles A and B at two loci are more frequent in a given population while alleles a and b are more frequent in another, dispersal will tend to produce an excess of AB, ab combinations in both populations (relative to Ab, aB combinations).

The effects of spatial heterogeneity on selection for recombination have been explored by [START_REF] Pylkov | Migration versus mutation in the evolution of recombination under multilocus selection[END_REF] and [START_REF] Lenormand | The evolution of recombination in a heterogeneous environment[END_REF] (see also [START_REF] Charlesworth | Selection on recombination in clines[END_REF], using models representing discrete demes connected by dispersal.

They showed that, with spatial differences in selection, recombination may increase the mean fitness of offspring (for example in the simple scenario mentioned above, if Ab, aB genotypes are fitter on average than AB, ab genotypes in both environments). As a consequence, recombination may be favored under wider ranges of parameters than in the absence of spatial heterogeneity, in particular when the strength of selection differs among demes [START_REF] Lenormand | The evolution of recombination in a heterogeneous environment[END_REF][START_REF] Otto | Resolving the paradox of sex and recombination[END_REF][START_REF] Otto | The evolutionary enigma of sex[END_REF]; see Agrawal, 2009 for similar results on the benefits of segregation in diploids).

These models have lead to very important insights, but are not ideally designed to explore the spread of asexual mutants in sexual populations living along an environmental continuum, as they generally assume recombination modifier loci with small effects and consider either a two-deme population or the island model of population structure.

Models of continuous populations living along environmental gradients have been used to explore the causes of species' range limits [START_REF] Kirkpatrick | Evolution of a species' range[END_REF]Barton, 2001;[START_REF] Bridle | Why is adaptation prevented at ecological margins? New insights from individual-based simulations[END_REF]Polechová and Barton, 2015;[START_REF] Polechová | Is sky the limit? On the expansion threshold of a species' range[END_REF]. The general framework is based on models representing the evolution of a polygenic, quantitative trait whose optimal value changes across space (Felsenstein, 1977;[START_REF] Slatkin | Spatial patterns in the distributions of polygenic characters[END_REF]Barton, 1999), to which a demographic component is added. The most commonly used scenario considers a one-dimensional habitat, with Gaussian stabilizing selection around a phenotypic optimum whose value changes linearly across space. In the case where all parameters (except the phenotypic optimum) stay constant over space, and neglecting the effects of drift and mutation, one can show that the loci coding for the trait should be at linkage equilibrium: stabilizing selection tends to produce negative genetic associations, whereas migration generates positive associations (as each location receives an excess of alleles increasing the trait from one direction, and an excess of alleles decreasing the trait from the other), and the two effects compensate exactly (Felsenstein, 1977;[START_REF] Slatkin | Spatial patterns in the distributions of polygenic characters[END_REF]Polechová and Barton, 2015). Therefore, one would expect that, in the absence of linkage disequilibrium, recombination should become neutral in this type of model. However, Felsenstein's result strictly holds in the case of an infinite, linear gradient, Gaussian stabilizing selection and no drift or mutation, and it is unclear how it may be affected by deviations from these idealized conditions.

In particular, selective pressures acting on recombination may differ at the range margins, due to the asymmetric gene flow coming from central parts of the population and to stronger drift effects, preventing adaptation to marginal conditions. Under the assumptions of the infinitesimal model of quantitative trait variation, [START_REF] Peck | Explaining the geographic distributions of sexual and asexual populations[END_REF] showed that asexual reproduction may be favored at the margins, since asexuality may maintain adaptation to extreme conditions in the face of a maladapted gene flow: this corresponds to one of the proposed hypotheses to explain geographic parthenogenesis.

Ecological, demographic and genetic parameters influencing adaptation at range limits (e.g., slope of the environmental gradient, dispersal distance of individuals, strength of selection relative to drift at loci involved in local adaptation) should affect the benefit of asexuals at range margins, but this has not been formally explored.

In this paper, we extend previous models of quantitative trait evolution along environmental gradients to explore how environmental heterogeneity may affect reproductive mode evolution in central and marginal parts of the distribution range of a species. We first consider the cases of fully sexual and fully asexual populations, in order to explore how the reproductive mode of organisms may affect the level of adaptation and components of genetic variance maintained within populations at equilibrium. In a second step, we explore how the different parameters of the model affect the spread of asexual mutants occurring in sexual populations. The results show that either sexual or asexual reproduction may be favored in the central part of the range depending on parameter values; however, the strength of indirect selection stays weak and is quickly overwhelmed by any direct fitness effect associated with the reproductive mode. By contrast, maladaptation at range margins tends to favor the local spread of asexual clones with fitter genotypes. This benefit of asexuals is stronger for parameter values leading to an important migration load at range limits (in particular, when the effect of drift is important at loci coding for the trait), and can maintain a pattern of geographical parthenogenesis when other indirect or direct forces maintain sex in the central part of the range. However, our results show that in this case, asexuals are generally restricted to small areas at the range margins.

METHODS

General framework. Our multilocus individual-based simulation program (written in C++, and available from Zenodo) is based on Polechová and Barton's (2015) model of adaptation along an environmental gradient, to which variability in reproductive mode is added. We provide in this section a general description of the program (whose main parameters are summarized in Table 1), while more details can be found in the Supplementary Material. coded by a finite number L of biallelic loci with two possible alleles (denoted 0 and 1), allele 1 increasing the value of the trait by an amount α: the trait value is given by z = α L i=1 X i , where X i equals 1 if the individual carries allele 1 at locus i, and 0 otherwise. Mutation occurs at a rate u per generation between the two alleles at each locus (in both directions); the overall mutation rate on the trait is thus U = uL.

The optimal trait value increases linearly across space: θ = bj + θ 0 , where j is the deme number (from 1 to n d ) and b the slope of the environmental gradient, while θ 0 is set so that the optimal trait value is αL/2 in the mid-point of the 1D-habitat (and is thus achieved in individuals carrying allele 1 at half of the loci coding for the trait).

Individuals are haploid and produce gametes by mitosis; female gametes are fertilized by male gametes produced by individuals from the same deme, the contribution of each individual to the pool of male gametes being proportional to its fitness W (densitydependence and local adaptation thus affect the overall reproductive output of each individual, through both male and female components). Gamete fusion occurs at random, including the possibility of selfing. Fertilization is immediately followed by meiosis to produce haploid juveniles (one per diploid zygote); we assume that the L loci coding for the trait are evenly spaced along a linear genome with map length R (average number of cross-overs at meiosis). Finally, juveniles disperse to neighboring demes according to a Gaussian dispersal kernel (discretized and truncated as explained in Appendix 3 of [START_REF] Polechová | Speciation through competition: a critical review[END_REF], where σ 2 measures the variance of the distance between parent and offspring (individuals move to demes 1 or n d if their dispersal distance would lead them to move beyond the metapopulation). We also considered a different life cycle in which dispersal occurs at the gametic stage, the different events thus occurring in the order selection → dispersal → recombination (rather than selection → recombination → dispersal), but both models yielded very similar results.

Previous results and parameterization. Some of the main predictions of the model in the case of a sexual population can be found in Barton (2001) and Polechová and Barton (2015), and will be only briefly summarized here. When the environmental gradient is not too steep and genetic drift is not too strong (relative to selection) at loci coding for the trait, the population can adapt to the gradient and increase its range. In that case, a succession of clines form along space at loci coding for the trait, so that the mean population phenotype z matches the optimal phenotype in all locations. At each locus, selection pushes towards the fixation of the most common allele with a strength s = α 2 / (2V s ) while dispersal tends to homogenize allele frequencies between demes, the balance between these two forces maintaining clines of width 4σ/ √ 2s (Barton, 1999(Barton, , 2001)). Neglecting the effects of drift and mutation, one predicts that the different loci coding for the trait should be at linkage equilibrium (as the effect of stabilizing selection generating negative associations exactly compensates the effect of dispersal generating positive associations, e.g., Felsenstein, 1977;Polechová and Barton, 2015), and that the genetic variance V g (the variance of z) equilibrates to V g = bσ √ V s in all places. This genetic variance corresponds to a "migration load" (reduction in mean fitness caused by maladapted immigrants) reducing the equilibrium deme size, given by:

N eq = K 1 - V g 2r m V s = K 1 - b σ 2r m √ V s . (1)
When the parameters stay constant over space, the population may expand indefinitely, as long as genetic variation is available (Barton, 2001). A range limit may be generated if the effective environmental gradient steepens (in particular, if b or σ increase), or if the strength of drift relative to selection increases (in particular, if K or r m decrease, Polechová and Barton, 2015). In our model, we generate range limits by imposing a decrease in carrying capacity K at both extremities of the metapopulation:

K stays equal to K max in the central part, and decreases linearly towards zero on both sides, over a distance corresponding to ω demes (see Figure 1). This allows us to compare evolutionary patterns in the central part of the metapopulation (where K = K max ) and in the margins (areas where K declines), in which gene flow is asymmetric (from more densely to less densely populated demes). We also used a model in which the range limits were generated by decreasing the per capita growth rate r m , but the results obtained were similar and will not be shown here.

The main parameters of the model are the slope of the environmental gradient b, the strength of stabilizing selection on the trait V s , the strength of selection acting at each locus s when z = θ (the effect of allelic change α being set to √ 2V s s), the variance of the distance between parent and offspring σ 2 and the carrying capacity in the central part of the population K max (together with ω that determines the slope of K at the margins). Some of the other parameters are expected not to affect the evolutionary outcome, at least as long as they are sufficiently large: this is the case in particular for the total number of demes n d , as long as it is large enough so that the amount of genetic variation maintained in the central part of the population is not affected by edge effects. For this, the program ensures that the overall range (number of demes) of the population is always at least 3 cline widths (see Supplementary Methods), while the number of demes in the central part (with K = K max ) is at least 100. The program also adjusts the total number of loci affecting the trait (L) so that genotypes coding for the optimal trait value are possible in all demes. Because the average strength of linkage disequilibria is expected to depend on the harmonic mean recombination rate ρ h between pairs of loci (e.g., [START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF][START_REF] Charlesworth | Directional selection and the evolution of sex and recombination[END_REF], the genome map length R is adjusted in order to ensure a given value of ρ h (generally 0.35) using equation A1 in [START_REF] Blanckaert | Epistasis, pleiotropy and the mutation load in sexual and asexual populations[END_REF]. Initialization is done as in Polechová and Barton (2015): the population initially occupies the central zone (where

K = K max )
and is well adapted in this area; it is then allowed to expand until it reaches its range limits (see Supplementary Material for more details).

Components of the genetic variance and migration load. We first considered separately the cases of obligately sexual and obligately asexual populations. Asexual populations follow the same life cycle as described above, except that haploid adults produce new individuals mitotically (apomixis). In both cases, the program records the mean trait value z and genetic variance V g in all demes. The genetic variance can be decomposed into two terms:

V g = α 2 L i=1 p i q i + α 2 i =j D ij (2)
(e.g., [START_REF] Lynch | Genetics and Analysis of Quantitative Traits[END_REF], where the first term depends on the frequencies of alleles 1 and 0 at each locus i (denoted p i and q i ) and corresponds to the "genic" variance (denoted V g,0 thereafter), while the second (denoted D thereafter) corresponds to the overall effect of linkage disequilibria between loci (D ij ) on the variance. Allele frequencies and V g are measured at the different steps of the life cycle in the different demes, in order to quantify V g,0 and D under the two different reproductive modes.

The genetic load corresponds to the reduction in the mean fitness of a population due to the presence of suboptimal genotypes, and is classically defined as

L = 1 -W /W max ,
where W is the mean fitness and W max the fitness of an optimal genotype (e.g., [START_REF] Crow | Genetic loads and the cost of natural selection[END_REF]. In our model, L is generated by recurrent mutation creating new genotypes (mutation load) and by the immigration of maladapted individuals from other demes (migration load); however, the contribution of mutation should stay negligible as long as U is small, and the migration load may thus be measured as:

L m = 1 -E exp - (z -θ) 2 2V s (3)
where E stands for the average over all individuals in a deme. When selection is

weak (V g V s ), the average in equation 3 is approximately equal to exp[-V g / (2V s )] × exp -(z -θ) 2 / (2V s ) ,
showing that the load may be decomposed into two components: the effect of genetic variance V g , and the effect of displacements of the mean phenotype from the optimum (e.g., [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF]. The migration load and its two components were measured in sexual and asexual populations at equilibrium, contrasting central and marginal populations.

In a second step, we considered the case of an initially sexual population, in which asexual mutants occur at a rate µ per generation (fixed to 10 -3 ; the program also includes back mutations restoring sex, also at a rate µ -see Supplementary Methods).

A direct cost (or benefit) of sex is introduced by multiplying the number of offspring produced by asexual females by a parameter c: c = 1 corresponds to the situation where sexuals and asexuals have the same baseline fitness (no cost of sex), while asexual mutants have a direct advantage when c > 1, and a disadvantage when c < 1.

We assume that asexuals do not produce any male gamete (no contagious asexuality), and only produce offspring by apomixis. As for the load, the change in the frequency of sexuals (denoted p sex ) due to indirect selection can be decomposed into two terms, representing the effects of differences in genetic variance V g and differences in mean phenotype z between sexuals and asexuals (sometimes called the "short-term" and "long-term" effect of recombination, e.g., [START_REF] Charlesworth | Directional selection and the evolution of sex and recombination[END_REF][START_REF] Barton | A general model for the evolution of recombination[END_REF][START_REF] Otto | The evolution of recombination in changing environments[END_REF][START_REF] Agrawal | Evolution of sex: why do organisms shuffle their genotypes?[END_REF][START_REF] Vanhoenacker | Stabilizing selection, mutational bias, and the evolution of sex[END_REF]. Based on previous works [START_REF] Barton | A general model for the evolution of recombination[END_REF][START_REF] Vanhoenacker | Stabilizing selection, mutational bias, and the evolution of sex[END_REF], these two components can be quantified as:

S Vg ≈ - 1 2 (V g + V s ) E (X sex -p sex ) (z -z) 2 (4) S z ≈ - z -θ V g + V s E [(X sex -p sex ) (z -z)] (5)
where X sex equals 0 and 1 in asexual and sexual individuals, respectively. The program measures the moments appearing on the right-hand sides of equations 4 and 5 in the different demes, thus providing estimates of the relative strengths of the two components of indirect selection for sex (S Vg and S z ).

RESULTS

Adaptation to environmental heterogeneity in sexual and asexual populations. We first consider the case of obligately sexual or asexual populations (without letting the reproductive mode evolve). In regimes where the population can adapt to the environmental gradient, the amount of genetic variation maintained at each location is predicted to increase with the slope of the gradient and with the dispersal distance of individuals, and to decrease with the strength of stabilizing selection (e.g., Barton, 1999). Assuming that the mean phenotype matches the optimum everywhere and that selection acting on the trait is weak (V g V s , so that differences in fitness between individuals remain small), the per-generation change in genetic variance Barton, 1999). If linkage disequilibria stay negligible in sexual populations (so that the effect of recombination can be neglected), and if the mutational variance V m = α 2 U is small relative to (bσ) 2 , the equilibrium genetic variance is given by ∆

due to stabilizing selection is ∆ sel V g ≈ -V g 2 /V s , while the change due to dispersal is ∆ disp V g ≈ (bσ) 2 (
sel V g ≈ -∆ disp V g , leading to V g ≈ bσ √ V s .
This shows that gene flow may substantially inflate the genetic variance, so that the hypothesis that V g V s may not necessarily hold. In that case, more accurate expressions can be obtained by assuming that the distribution of phenotypes in each location is approximately Gaussian: the [START_REF] Bulmer | The Mathematical Theory of Quantitative Genetics[END_REF], so that the equilibrium value of V g after dispersal, before selection is given by:

change in V g caused by selection then becomes ∆ sel V g ≈ -V g 2 / (V g + V s ) (e.g.,
V g, aft disp ≈ 1 2 (bσ) 2 + bσ (bσ) 2 + 4V s (6)
while the variance after selection, before dispersal is:

V g, aft sel ≈ V g, aft disp -(bσ) 2 , (7)
the average between both quantities being (bσ/2) (bσ) 2 + 4V s (and thus converging to bσ √ V s when (bσ) 2 V s , so that V g V s ). Figure 2 shows that the predictions from equations 6 and 7 match our simulation results for all tested values of b, σ and V g (this was verified for both sexual and asexual populations, as well as for the selection → dispersal → recombination life cycle). The same results are obtained for different strengths of selection per locus s, including values such that N eq s 1 (so that the effect of drift is strong relative to selection at each locus) and N eq s 1.

Despite the fact that sexual and asexual populations display approximately the same amount of genetic variance at equilibrium, slight differences are observed.

In particular, while the average variance usually matches very well the theoretical prediction in asexuals (except for high values of σ), sexuals may display slightly higher or lower V g (Figure S1). These differences are due to the effect of recombination:

although linkage disequilibria are small at the time when recombination occurs, they may be slightly positive or negative depending on parameter values (Figure 3). When linkage disequilibria are positive (D > 0), recombination tends to decrease V g by reducing the second term of equation 2, while recombination has the opposite effect when D < 0. Figure 3A and 3B indicate that positive D arises when the per locus strength of selection s is small relative to the effect of drift, while D becomes negative as deme size or s increases. Although we could not obtain an analytical expression for D at migration-selection-drift equilibrium, insights can be gained from quasi-linkage equilibrium results under weak selection: indeed, when D remains small and V g V s ,

the change in D due to selection is approximately -(V g,0 ) 2 /V s (e.g., [START_REF] Turelli | Dynamics of polygenic characters under selection[END_REF], while the change due to dispersal is (bσ) 2 . As long as D remains small, and ignoring the effects of mutation and drift, V g,0 ≈ V g ≈ bσ √ V s , so that the two effects compensate exactly (Felsenstein, 1977;Polechová and Barton, 2015). Mutation tends to increase the genic variance V g,0 , so that the effect of selection becomes stronger than the effect of dispersal, generating negative D (Slatkin, 1978, Figure 3F). By contrast, drift tends to decrease V g,0 , generating positive D. As shown by Figure 3, σ and V s have only limited effects on the contribution of linkage disequilibria to V g , while increasing b tends to increase the relative importance of positive D in the weak selection regime.

In asexuals, simulations indicate that while the genetic variance V g remains constant over time, D tends to become more and more negative (while V g,0 becomes more and more positive), due to the constant appearance of new genotypes carrying different combinations of compensatory mutations (the same effect is predicted to occur in the case of fully selfing populations, e.g., [START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF][START_REF] Lande | Maintenance of quantitative genetic variance under partial self-fertilization, with implications for the evolution of selfing[END_REF]. Although the observed V g matches well the theoretical prediction (equation 6)

when averaged over demes (Figure S1), we observed that a stable pattern of periodic oscillations of V g and z-θ sets up along the metapopulation for relatively low values of s and/or high values of σ (e.g., for s = 0.001 and σ = 6 in Figure 4). This corresponds to a situation in which optimal genotypes are predominant only in some demes separated by a fixed distance: in those demes, z = θ and V g reaches a minimal value, while the demes in between contain mixtures of maladapted immigrant genotypes, causing V g to increase and z to deviate from θ (in this situation, D is initially positive and becomes negative over time, as shown by Figure S2). As shown by Figure S3, the amplitude of these oscillations increases as deme size decreases. This clustering of the phenotypic distribution of asexuals along an environmental gradient has already been shown by [START_REF] Polechová | Speciation through competition: a critical review[END_REF] using a deterministic version of the same model (see also [START_REF] Champagnat | Invasion and adaptive evolution for individualbased spatially structured populations[END_REF], and is due to a structural instability of the smooth solution with constant V g and z = θ when the mutational variance is small (in other words, two asexual clones with very similar phenotypes cannot coexist in adjacent demes). Interestingly, in this regime the phenotypic distribution is multimodal and thus non-Gaussian in most of the demes, but the Gaussian approximation still provides a correct prediction of the average genetic variance over demes (Figures 2,S1).

As shown by Figure 4, periodic fluctuations do not occur in the case of sexual populations: V g stays approximately constant over space, for all tested parameter values. In both sexuals and asexuals, mean phenotypes usually stay close to the optimum in the central part of the population (although some deviations arise in the oscillating regime just described), so that the migration load is mostly generated by V g . Therefore, sexuals have a slightly higher load than asexuals in the regime where D < 0 in sexuals (recombination increases V g ) and a slightly lower load in the regime where D > 0 (recombination decreases V g , see Figure 5). However, the oscillations of V g occurring in asexuals in the high σ / low s regime generate oscillations of the load, as shown by Figure S4.

At range margins, the lower population size causes a decrease in the genic variance V g,0 in sexual populations. Interestingly, we generally observed that this decrease in V g,0 is almost exactly compensated by an increase in D, so that V g stays approximately constant in most of the range, and remains well predicted by equation 6 (see Figure S5 for an example). Moreover, Figure 6 shows that the mean phenotype tends to depart from the optimum at range limits (due to drift and to asymmetric gene flow from more densely populated demes), becoming an important component of the load (in particular when selection per locus s is weak, as shown by Figures 6, S6 and S7). Maladaptation at range limits is more important in the case of asexual populations (Figures 6, S6, S7), probably due to the increased effects of drift caused by interference between selected loci in the absence of recombination (e.g., Barton, 2010).

As a consequence, the species range tends to be smaller in asexuals than in sexuals (Figure S8), the difference being again mostly marked when s is small. Despite the fact that purely asexual populations tend to suffer more from migration load at range limits, we will see that asexual mutants can nevertheless be favored at the range limits of sexual populations, since they are not affected by migration load in this case.

Evolution of reproductive mode. When asexual mutants are introduced within sexual populations, we observed that in the absence of any direct fitness effect of the reproductive mode (c = 1), sexual reproduction stays predominant in the central part of the range for parameter values leading to a lower load in sexuals than in asexuals, while asexuality prevails for parameter values leading to a lower load in asexuals. In particular, Figure 7A shows that the frequency of sexuals (p sex ) stays high in the central part of the range when the per locus strength of selection s is small (leading to positive linkage disequilibria, Figure 3), while asexual mutants invade at higher values of s (Figure S9 shows that the effects of b, σ and V s on the frequency of sexuals can also be predicted from their effects on the relative migration load of sexuals and asexuals at equilibrium). However, in the parameter region where sex is favored (low s), asexuals are maintained at non-negligible frequencies (e.g., p sex ≈ 0.8 for s = 0.001). This may be due to the fact that differences in mean fitness between sexual and asexual populations are weak, but also to the oscillations occurring in asexual populations in the low s regime (Figure 4) which may favor asexuals over sexuals in certain demes.

Figure 7B confirms that selection on the reproductive mode is mainly driven by the difference in genetic variance V g between sexual and asexual populations (quantified by S Vg , given by equation 4), while the mean phenotype stays close to the optimum (on average) for both reproductive modes (leading to S z ≈ 0). Differences in V g and in the migration load between sexual and asexual population stay relatively weak (Figure 5), so that any direct fitness effect of the reproductive mode quickly prevails and leads to the near fixation of sexuals when c < 1 (direct cost of asexuality), and to the near fixation of asexuals when c > 1 (direct cost of sex, see Figure 8).

By contrast, the evolution of reproductive mode at the range limits cannot be predicted from the migration load of purely sexual and purely asexual populations at equilibrium. Indeed, while purely asexual populations always suffer from a similar or higher load than sexual populations at range limits (Figure S6), we observed that asexual clones may often spread in peripheral populations. When sex is maintained in the central part of the range (either due to a higher load in asexuals when c ∼ 1 or to an intrinsic advantage of sexuals, i.e., c < 1), the asexual population maintained at range limits does not suffer from migration load (since immigrants coming from the central area are mostly sexual), leading to a stable pattern of geographic parthenogenesis. The spread of asexuals at range margins is due to the maladaptation of sexual populations in these areas, generating an advantage for locally adapted mutant clones: indeed, asexuals are favored because they can maintain phenotypes that are closer to the optimum (S z < 0, see Figure S10). This advantage of asexuals increases with the degree of maladaptation of sexuals, which (as shown by Polechová and Barton, 2015) depends on the effective environmental gradient (proportional to bσ/ √ V s ) and on the effect of drift relative to selection at loci coding for the trait (see Figure 9). For sufficiently strong values of bσ/ √ V s and sufficiently weak selection relative to drift at loci coding for the trait (low s), asexuality may thus be maintained at range limits despite an intrinsic advantage of sexuals (c < 1, see Figures S11 -S16 showing the abundance of sexuals and asexuals along space for different parameter values). In general, the maintenance of geographic parthenogenesis in the present model is only possible for sufficiently weak selection at loci coding for the trait: under stronger selection (e.g., s = 0.05 in Figures S15,S16), the benefit of asexuals at range limits stays weak and is easily overwhelmed by any advantage associated with sex (represented here by c < 1), while such an advantage of sex is necessary to maintain sexuals in the central part.

Finally, Figures S11 -S16 show that the area in which asexuals are predominant at range limits (under conditions favoring sex in the central area) generally stays rather small. As shown by Figure 10A, decreasing the slope of the change in carrying capacity over space (by increasing ω) has little effect on the size of this area. Although we have seen that the advantage of asexuals at range limits increases with bσ/ √ V s (Figure 9), Figure 10B shows that asexuals tend to be confined in a smaller area (measured in units of dispersal distance σ) as bσ/ √ V s increases. This is likely due to the fact that although a locally adapted clone may benefit from a stronger advantage when bσ/ √ V s is higher (due to the stronger maladaptation of sexuals), this benefit is quickly lost as individuals move to neighboring demes (since the effective environmental gradient is stronger).

DISCUSSION

Temporal environmental change generally favors sex and recombination due to the fact that recombination increases the amount of genetic variation upon which direc-tional selection can act [START_REF] Charlesworth | Directional selection and the evolution of sex and recombination[END_REF][START_REF] Barton | A general model for the evolution of recombination[END_REF][START_REF] Otto | The evolution of recombination in changing environments[END_REF][START_REF] Otto | The evolutionary enigma of sex[END_REF]Becks and Agrawal, 2012). By contrast, the effect of spatial environmental change on selection for sex has been somewhat less studied. [START_REF] Pylkov | Migration versus mutation in the evolution of recombination under multilocus selection[END_REF] and [START_REF] Lenormand | The evolution of recombination in a heterogeneous environment[END_REF] demonstrated that spatial variation in the direction or strength of selection may sometimes generate a short-term benefit for recombination (increase in mean fitness of recombinant offspring), while [START_REF] Peck | Explaining the geographic distributions of sexual and asexual populations[END_REF] showed that asexuality may be favored at species range margins as it can preserve locally adapted genotypes. In this study, we used a classical model of evolution along an environmental gradient (Felsenstein, 1977;[START_REF] Slatkin | Spatial patterns in the distributions of polygenic characters[END_REF]Barton, 1999;Polechová and Barton, 2015) to further explore the effects of spatial heterogeneity on reproductive mode evolution. We found that in situations where demes are locally adapted (i.e., when the average phenotype matches the local optimum), recombination either increases or decreases the mean fitness of offspring depending on the relative strengths of selection, mutation and drift. While the effect of stabilizing selection generating negative LD exactly compensates the effect of dispersal generating positive LD when mutation and drift are neglected (Felsenstein, 1977;Polechová and Barton, 2015), mutation tends to produce negative LD [START_REF] Slatkin | Spatial patterns in the distributions of polygenic characters[END_REF] while drift generates positive LD (recombination being beneficial under positive LD as it tends to increase the proportion of locally adapted genotypes, and disadvantaged under negative LD). This leads to a benefit for sexual reproduction when deme size is sufficiently small, and when selection and mutation are sufficiently weak. However, in most cases the difference in mean fitness between sexuals and asexuals stays rather weak, so that indirect selection acting on the reproductive mode is easily overwhelmed by even slight direct costs or benefits associated with asexuality. This model therefore does not represent a realistic scenario for the evolutionary maintenance of costly sex within populations. However, it is interesting to note that while stabilizing selection generally favors asexuality in equilibrium populations living in homogeneous habitats (in the absence of mutational bias, e.g., [START_REF] Charlesworth | Directional selection and the evolution of sex and recombination[END_REF][START_REF] Vanhoenacker | Stabilizing selection, mutational bias, and the evolution of sex[END_REF], a spatial environmental gradient may erase this advantage of asexuals (at least in the central part of the species range), making it easier for sex to be favored through other mechanisms. Exploring the effects of selection, drift and dispersal on the sign and magnitude of linkage disequilibria (and the resulting strength of selection for sex and recombination) under more general scenarios (including different forms of fitness function or non-linear changes in phenotypic optimum over space) would be of interest, as they may generate stronger selection on recombination.

From an empirical perspective, a better understanding of the selective forces affecting the evolution of reproductive modes could be gained by measuring the effect of recombination on the mean and variance in fitness among offspring (e.g., [START_REF] Sharp | Evolution of sex: using experimental genomics to select among competing theories[END_REF]. This type of data remains scarce, however [START_REF] Peters | Liberating genetic variance through sex[END_REF][START_REF] Becks | The effect of sex on the mean and variance of fitness in facultatively sexual rotifers[END_REF] and references therein) and to our knowledge, only one study used fitness measures obtained under natural conditions [START_REF] Kelley | A test of the short-term advantage of sexual reproduction[END_REF]. Based on our results and those of previous models (Felsenstein, 1977;[START_REF] Lenormand | The evolution of recombination in a heterogeneous environment[END_REF], it would be of particular interest to compare the effects of recombination on offspring fitness in populations subject to different forms of environmental heterogeneity. In particular, while recombination is expected to increase the variance in fitness and decrease mean fitness under spatially homogeneous stabilizing selection, these effects may vanish or even reverse in the case of populations living along environmental gradients. This could be done in principle by comparing the fitness distributions of sexually and asexually produced offspring in species using both reproductive modes [START_REF] Kelley | A test of the short-term advantage of sexual reproduction[END_REF]Becks and Agrawal, 2012), although the fact that sexually and asexually produced propagules are often not ecologically or physiologically equivalent represents a challenge for such experiments.

Our simulations of obligately asexual populations showed that in some regimes, adaptation cannot be maintained at all geographic locations (a result that was already obtained in previous numerical explorations of similar models, in particular Polechová and Barton, 2005;[START_REF] Champagnat | Invasion and adaptive evolution for individualbased spatially structured populations[END_REF]. In such situations, the population However, our results show that geographic parthenogenesis is difficult to maintain in this type of model when the strength of selection acting at loci involved in local adaptation is not sufficiently small. Furthermore, asexuality usually stays restricted to a small area at range limits where density is low. While this may be compatible with some of the reported cases of geographic parthenogenesis (although estimates of σ are rarely available), situations in which asexuals occupy a very large area (e.g., [START_REF] Van Dijk | Ecological and evolutionary opportunities of apomixis: insights from Taraxacum and Chondrilla[END_REF][START_REF] Schmit | Linking present environment and the segregation of reproductive modes (geographical parthenogenesis) in Eucypris virens (Crustacea: Ostracoda)[END_REF] seem more difficult to explain with this type of hypothesis.

Demographic advantages associated with uniparental reproduction (better colonizing ability, [START_REF] Baker | Self-compatibility and establishment after 'long-distance' dispersal[END_REF] may possibly explain such patterns in the case of recent geographic expansions: such benefits of asexuality could be incorporated in the model by considering dioecious or self-incompatible sexual organisms.

While changes in species ranges and evolutionary dynamics at population expansion fronts have been the subject of increased attention in a context of accelerating global change (e.g., [START_REF] Gilbert | Local adaptation interacts with expansion load during range expansion: maladaptation reduces expansion load[END_REF][START_REF] Peischl | Expansion load and the evolutionary dynamics of a species range[END_REF][START_REF] Williams | How evolution modifies the variability of range expansion[END_REF], the selective pressures acting on the evolution of recombination or reproductive systems during range expansion remain little explored. Yet, shifts in reproductive system have profound demographic and genetic consequences that should in turn affect population expansion [START_REF] Barrett | Plant reproductive systems and evolution during biological invasion[END_REF]. Preliminary results obtained using our simulation models indicate that asexual mutants may spread at expansion fronts due to the im- space occurs in asexual populations when σ is high and s is small. population density quickly drops to zero in those demes; they can be later recolonized as the population expands.

� / � � � ��� ��� ��� ��� ��� ��� � ���� ���� ���� ���� � / � � � � � � � � � σ -���� ���� ���� ���� ���� ���� � / � � � �� �� �� �� �� � � -���� ���� ���� ���� � / � � � �� -� �� -� �� -� � -���� -���� -���� ���� ���� � / � � s=0.05 s=0.001
� ����� ���� ��� � -����� ����� ����� ����� ����� � � ���� -� � ��� � ��� ��� ��� ��� ��� ��� ��� � ��� ����� ����� ����� ����� � � ���� -� � ��� � ��� ��� ��� ��� ��� ��� � ����� ����� ����� ����� ����� � � ���� -� � ��� � � � � � � � σ -����� ����� ����� � � ���� -� � ��� � �� �� �� �� �� � � -����� -����� ����� ����� ����� ����� ����� � � ���� -� � ��� � �� -� �� -� �� -� � -����� -����� ����� ����� � � ���� -� � ���
Measures and numbers of generations. Genetic and genic variances are measured every 1,000 generations within each deme, after each step of the life cycle (selection, recombination and dispersal). The genetic variance V g corresponds to the variance of z within the deme, while the genic variance is given by V g,0 = α 2 L i=1 p i q i , where p i and q i are the frequencies of the two alleles at locus i. Values of V g and V g,0 after selection (before recombination) are obtained by weighting each individual by its relative fitness. The program also measures the mean phenotype z and migration load, given by equation 3. Except for the simulation results shown on Figure 9, a single replicate was performed for each set of parameters, over a large enough number of generations to ensure that an equilibrium was reached long before the end of the simulation. In the case of obligately sexual populations (without letting the reproductive mode evolve), simulations generally lasted 50,000 generations, and averages were performed over the last 40,000 generations (we checked by visual inspection that equilibrium was reached during the first 10,000 generations). For most simulations of obligately asexual populations, only phenotypic moments were measured and we used a modified program in which only the total number of alleles 1 in the genome of each individual is represented (without representing explicitly the L loci, in order to increase execution speed); however, measures of V g,0 and D (Figure S2) were obtained for some parameter values using the standard program. Because asexual populations took longer to reach their range limit, simulations were run over a larger number of generations (500,000 in the case of Figures 5,6, S6-S8). Simulations in which we let the reproductive mode evolve generally lasted 90,000 generations, asexual mutants being introduced after 60,000 preliminary generations (each sexual individual becoming asexual with probability µ = 10 -3 ).

Selection -dispersal -recombination. A different program represented gametic dispersal, so that recombination occurs just after dispersal, before selection (instead of occurring after selection, before dispersal as in our standard program). In that case, female gametes migrate first according to the dispersal kernel. The probability that a female gamete present in deme i is fertilized by a male gamete immigrating from deme j is then given by W ♂j m ij / k W ♂k m ik where W ♂j is the average male fitness in deme j, and m ij the probability that a gamete produced in deme j migrates to deme i. Once the deme of origin of the male gamete has been drawn, a father is sampled according to its male fitness as in the standard program. The results obtained using this modified program were very similar to those obtained under the selectionrecombination -dispersal life cycle (not shown). . Fluctuations of the genetic variance V g over space in asexual populations, at generation 500,000 (same parameter values and colors as in Figure S2). The migration load is averaged over the demes located at both edges of the range, in which deme size is comprised between 10 and 50 individuals.

SUPPLEMENTARY FIGURES

� ����� ���� ��� � -����� -����� ����� ����� (� � -� � �� )/ (�� � ) ���� ��� � ��� ��� ��� ��� ��� ��� ��� � ��� -����� -����� -����� -����� -����� ����� ����� (� � -� � �� )/ (�� � ) � ��� ��� ��� ��� ��� ��� � -����� -����� -����� ����� ����� (� � -� � �� )/ (�� � ) � � � � � � � σ -����� -����� -����� -����� ����� ����� (� � -� � �� )/ (�� � ) � �� �� �� �� �� � � -����� -����� -����� -����� -����� ����� ����� (� � -� � �� )/ (�� � ) � �� -� �� -� �� -� � ����� ����� ����� ����� (� � -� � �� )/ (�� � ) asex, s = 0.
� � ���� ���� -� � ���� ��� � � � � � � � σ ��� ��� ��� ��� ��� � � ���� ���� -� � ���� ��� � �� �� �� �� �� � � ���� ���� ���� � � ���� ���� -� � ���� ���
� ����� ���� ��� � -���� -���� -���� -���� ���� (� � ��� -� � �� )/ (�� � ) � ��� ��� ��� ��� ��� ��� � -��� -��� -��� -��� -��� ��� ��� (� � ��� -� � �� )/ (�� � ) � � � � � � � σ -��� -��� ��� (� � ��� -� � �� )/ (�� � ) � �� �� �� �� �� � � -���� -���� ���� (� � ��� -� � �� )/ (�� � )

  The model represents a population inhabiting a onedimensional habitat consisting of n d demes connected by dispersal. Individuals are haploid and hermaphrodite; generations are discrete. Each generation, the number of female gametes produced by an individual is sampled from a Poisson distribution with parameter W = e r , where r represents the Malthusian fitness of the individual and is the sum of two components: a density-dependent term r N = r m 1 -N K , where N is the size of the local deme, K its carrying capacity and r m the maximum per capita growth rate, and a local adaptation component r z = -(z -θ) 2 / (2V s ), where z represents the value of a phenotypic trait of the individual, θ the optimal value of the trait and V s the strength of stabilizing selection (r = r N + r z ). The trait z is

  consists of different subtypes adapted to specific locations separated by a more or less constant distance (which increases with the dispersal distance of individuals). As a consequence, the genetic variance and level of adaptation fluctuate along space in a periodic fashion (Figure4). Comparing the patterns of adaptation along environmental gradients in sexual and asexual populations represents another possible direction for empirical studies, either by comparative analyses of natural populations or by using experimental evolution approaches. At range limits, our simulations indicate that obligately asexual populations generally suffer from a stronger migration load (caused by migration bias from more densely populated habitats) than obligately sexual populations, in particular when the strength of selection per locus is small. Similarly,[START_REF] Vanhoenacker | Stabilizing selection, mutational bias, and the evolution of sex[END_REF] showed that mutational biases tend to have a stronger impact on the mutation load of asexual populations, due to the reduced efficiency of directional selection in asexuals. Extending our model to include a migration bias in the central part of the population (that may be caused by winds or ocean currents in natural habitats) could thus be of interest, as it may increase the strength of selection for sex.Finally, our results confirm that the constant maladaptation of peripheral populations can favor the spread of asexual mutants adapted to these extreme environments, which may lead to a pattern of geographic parthenogenesis when sex is maintained in the central part of the range (either by direct or indirect selection). The strength of selection for asexuality correlates positively with the magnitude of the migration load in peripheral areas, which itself increases with the steepness of the environmental gradient (b), the variance of the dispersal distance (σ 2 ) and the strength of stabilizing selection (1/V s ), but decreases with the strength of selection per locus (s). This is consistent with the results obtained by[START_REF] Peck | Explaining the geographic distributions of sexual and asexual populations[END_REF] using the infinitesimal model (extremely small strength of selection per locus), also showing patterns of geographic parthenogenesis in the presence of a direct fitness cost associated with asexuality.

Figure 1 .�Figure 2 .

 12 Figure 1. Variation of the carrying capacity K across space: K equals K max in the

Figure 3 .Figure 4 .

 34 Figure 3. Contribution of linkage disequilibria to the genetic variance D/V g (where

Figure 5 .

 5 Figure 5. Difference in migration load L m (measured by equation 3) between asexual

Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .

 678910 Figure 6. Squared difference between the mean phenotype and the optimum, scaled

Figure S2 .

 S2 Figure S1. Difference between the average genetic variance V g per deme observed in

Figure S4 .Figure S5 .

 S4S5 Figure S4. Variation of the migration load along space in sexuals (black) and asexuals

Figure S6 .

 S6 Figure S6. Difference in migration load L m between asexual and sexual populations

Figure S7 .Figure S8 .Figure S9 .SFigure S11 .Figure S12 .

 S7S8S9S11S12 Figure S7. Difference between the average genetic variance V g per deme observed in

Figure S13 .Figure S14 .Figure S15 .Figure S16 .

 S13S14S15S16 Figure S13. Number of sexual (black) and asexual (grey) individuals per deme at

Community (project MARFOR Biodiversa/004/2015) and Region Bretagne (ARED 2017 REEALG) for funding, and the Ecole Doctorale "Science de la Nature et de l'Homme, Ecologie et Evolution" for a Transhumance mobility grant to LF. Data availability: The code for the C++ simulation program is available on Zenodo (DOI 10.5281/zenodo.4643637).

SUPPLEMENTARY METHODS

Initialization. All demes are initially populated with N eq individuals, given by equation 1 (r m is readjusted if needed to ensure that N eq ≥ 30 in the central part of the population). As in Polechová and Barton (2015), a series of clines is set up at loci coding for the quantitative trait in the central part of the population, in order to ensure that demes in the central part are already well adapted. At migration-selection equilibrium, the equation for each cline is given by:

where p is the frequency of allele 1, x the spatial position and x 0 the position of the center of the cline (at which p = 0.5), and where the cline width is given by