Audrey Lyrata Audrey Le Veve

Nicolas Burghgraeve

Mathieu Genete

Christelle Lepers-Blassiau

Margarita Takou

Juliette De Meaux

Barbara K Mable

Eléonore Durand

Xavier Vekemans

Vincent Castric
email: vincent.castric@univ-lille.fr

Long-term balancing selection and the genetic load linked to the self-incompatibility locus

Keywords: balancing selection, linked selection, sheltered genetic load, deleterious mutations, S-locus, Arabidopsis, self-incompatibility

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

Introduction

Balancing selection refers to a variety of selective regimes maintaining advantageous genetic diversity within populations [START_REF] Delph | On the importance of balancing selection in plants[END_REF]. Notable examples include heterozygote advantage, negative frequency-dependent selection, and spatial heterogeneity. The implication of balancing selection in maintaining genetic diversity across the genome has been widely debated (e.g. [START_REF] Asthana | A limited role for balancing selection[END_REF]. In contrast to genetic linkage to genomic sites subject to either positive or negative selection that generally tends to eliminate surrounding genetic variation (Smith and Haigh 1974, for hitchhiking effect;[START_REF] Charlesworth | The effect of deleterious mutations on neutral molecular variation[END_REF][START_REF] Loewe | Background selection in single genes may explain patterns of codon bias[END_REF]Charlesworth 2007, for background selection effects), linkage to loci under balancing selection is expected to locally promote the long-term persistence of variation in surrounding sites [START_REF] Charlesworth | Balancing selection and its effects on sequences in nearby genome regions[END_REF]. Theoretical studies by [START_REF] Takahata | Footprints of intragenic recombination at HLA loci[END_REF], [START_REF] Schierup | The effect of hitch-hiking on genes linked to a balanced polymorphism in a subdivided population[END_REF] and [START_REF] Wiuf | The probability and chromosomal extent of trans-specific polymorphism[END_REF] showed that besides the strength of balancing selection and the local rate of recombination, the magnitude of the local diversity increase and its extent along the chromosome critically depend on details of the exact form of balancing selection [START_REF] Llaurens | Genetic architecture and balancing selection: the life and death of differentiated variants[END_REF] for a review). An important result from these studies is that the extended time over which the balanced allelic lineages are maintained [START_REF] Takahata | Allelic genealogy under overdominant and frequency-dependent selection and polymorphism of major histocompatibility complex loci[END_REF][START_REF] Vekemans | Gene and allelic genealogies at a gametophytic self-incompatibility locus[END_REF]) also means more time for recombination to decouple them from their linked sites, such that the extent of the region affected may end up being quite narrow [START_REF] Kaplan | The coalescent process in models with selection and recombination[END_REF][START_REF] Schierup | Recombination, balancing selection and phylogenies in MHC and self-incompatibility genes[END_REF]). In addition to the sheer increase of polymorphism, several balancing selection processes promote heterozygosity. They are thus expected to mask recessive deleterious mutations [START_REF] Maruyama | Genetic variability maintained by mutation and overdominant selection in finite populations[END_REF], such that linkage to a locus under balancing selection can negatively interfere with purifying selection, diminishing its efficacy and facilitating the local accumulation of a potentially strong genetic load, referred to as the "sheltered load" [START_REF] Uyenoyama | Genealogical structure among alleles regulating self-incompatibility in natural populations of flowering plants[END_REF][START_REF] Uyenoyama | Evolution under tight linkage to mating type[END_REF][START_REF] Hartfield | Recombination and hitchhiking of deleterious alleles[END_REF]. This phenomenon has been considered as the "evolutionary cost" of balancing selection [START_REF] Van Oosterhout | A new theory of MHC evolution: beyond selection on the immune genes[END_REF][START_REF] Lenz | Excess of deleterious mutations around HLA genes reveals evolutionary cost of balancing selection[END_REF], and in humans a large number of diseases are indeed associated with variants at genes linked to one of the classical examples of balancing selection in the human genome, the Major Histocompatibility Complex (MHC; e.g. [START_REF] Lenz | Excess of deleterious mutations around HLA genes reveals evolutionary cost of balancing selection[END_REF][START_REF] Matzaraki | The MHC locus and genetic susceptibility to autoimmune and infectious diseases[END_REF].

Evaluating the importance of balancing selection and determining its evolutionary consequences has been the focus of sustained interest in the field [START_REF] Llaurens | Genetic architecture and balancing selection: the life and death of differentiated variants[END_REF]. Genome resequencing studies have revealed that balancing selection can be a potent force throughout the genome [START_REF] Degiorgio | A model-based approach for identifying signatures of ancient balancing selection in genetic data[END_REF], but because of the inherent technical challenges related to the high levels of polymorphism in the genomic regions affected [START_REF] Vekemans | Whole-genome sequencing and genome regions of special interest: Lessons from major histocompatibility complex, sex determination, and plant self-incompatibility[END_REF] it is still unclear how widespread the various forms of balancing selection actually are (see e.g. [START_REF] Fijarczyk | Detecting balancing selection in genomes: limits and prospects[END_REF]. In support of the sheltered load hypothesis, [START_REF] Lenz | Excess of deleterious mutations around HLA genes reveals evolutionary cost of balancing selection[END_REF] observed a specific accumulation of putatively deleterious mutations (missense variants) in genes that are located inside the human MHC region but have no function in immunity and just happen to be linked to the MHC alleles. Interestingly, this sheltered load was mostly due to an increase in the mean population frequency of deleterious mutations as compared to genes in a series of "control" regions, but not to an elevation of their overall number, suggesting that the balancing selection process at play for the human MHC region elevates polymorphism locally by distorting the frequency of deleterious mutations, rather than by increasing their density. Whether this observation can be generalised to other biological systems under balancing selection is not known.

Self-incompatibility (SI) in plants is perhaps the best understood case of long-term balancing selection [START_REF] Castric | Plant self-incompatibility in natural populations: a critical assessment of recent theoretical and empirical advances[END_REF]. SI is a genetic mechanism allowing recognition and rejection of self-pollen, thereby preventing inbreeding and promoting outcrossing in hermaphroditic plants [START_REF] Nettancourt | Incompatibility and incongruity in wild and cultivated plants[END_REF]. Pollination between partners expressing identical haplotypes at the S-locus leads to rejection of the pollen. This genetic system enforces outcrossing and promotes higher heterozygosity than expected under random mating. In addition, as noted by [START_REF] Wright | The distribution of self-sterility alleles in populations[END_REF], pollen produced by individuals carrying rare S-alleles will more rarely land on incompatible pistils than pollen produced by individuals carrying S-alleles that are more frequent. The action of natural selection on the S-locus is thus well elucidated and corresponds to an intense form of negative frequency-dependent selection that allows the stable maintenance of a large number of S-alleles within populations. The S-alleles are maintained over very long evolutionary times [START_REF] Vekemans | Gene and allelic genealogies at a gametophytic self-incompatibility locus[END_REF], and theoretical models predict that a local increase of nucleotide polymorphism should be observed in the linked genomic region [START_REF] Uyenoyama | Genealogical structure among alleles regulating self-incompatibility in natural populations of flowering plants[END_REF][START_REF] Schierup | The effect of hitch-hiking on genes linked to a balanced polymorphism in a subdivided population[END_REF]. In gametophytic SI (GSI), pollen SI specificity is determined by its own haploid genome (as found e.g. in Solanaceae), whereas in sporophytic SI (SSI), the pollen recognition phenotype is determined by the male diploid parent (as found e.g. in Brassicaceae). In the Brassicaceae, SSI is controlled by a single genomic region, the S-locus [START_REF] Schopfer | The male determinant of self-incompatibility in Brassica[END_REF][START_REF] Kusaba | Self-incompatibility in the genus Arabidopsis: characterization of the S locus in the outcrossing A. lyrata and its autogamous relative A. thaliana[END_REF]), composed of two linked genes, SCR (encoding the S-locus cysteine-rich protein) and SRK (encoding the S-locus receptor kinase protein), encoding the male and female specificity determinants, respectively.

The phenotypic effect of the sheltered genetic load linked to the S-locus can be revealed by controlled crosses to experimentally enforce homozygosity at the S-locus and isolate the specific effect of this homozygosity on proxies of fitness. To the best of our knowledge, such experiments have been performed in five plant species only: Papaver rhoeas (Lane and Lawrence 1995), Solanum carolinense [START_REF] Stone | Sheltered load associated with S-alleles in Solanum carolinense[END_REF]), Arabidopsis halleri [START_REF] Llaurens | The sheltered genetic load linked to the S locus in plants: new insights from theoretical and empirical approaches in sporophytic self-incompatibility[END_REF], A. lyrata [START_REF] Stift | Inbreeding depression in self-incompatible North-American Arabidopsis lyrata: disentangling genomic and S-locus-specific genetic load[END_REF] and Rosa [START_REF] Vieira | The identification of the Rosa S-locus and implications on the evolution of the Rosaceae gametophytic self-incompatibility systems[END_REF]. In these species, a detectable genetic load linked to the S-locus could be revealed, although its magnitude varied among the S-alleles that were brought to the homozygous state. The fact that this load was detectable at the phenotypic level in spite of the inherently limited experimental power of these studies, suggests that the S-locus does indeed shelter a substantial load of deleterious mutations. These studies focused on phenotypic characterization of the load, and thus provided no indication about its genomic architecture. At this stage, the nature of this load therefore remains elusive. In A. halleri, the S-locus has been sequenced entirely in multiple haplotypes, revealing that the non-recombining S-locus region contains no protein-coding genes besides the ones controlling the SI machinery itself (SCR and SRK; [START_REF] Goubet | Contrasted patterns of molecular evolution in dominant and recessive self-incompatibility haplotypes in Arabidopsis[END_REF]. The load detected phenotypically is therefore likely caused by mutations in the partially linked flanking regions rather than in the non-recombining S-locus region itself. A series of studies have set out to determine the genomic extent of the flanking region over which polymorphism was altered by linkage to the S-locus in A. halleri and A. lyrata [START_REF] Kamau | Balancing selection and low recombination affect diversity near the self-incompatibility loci of the plant Arabidopsis lyrata[END_REF][START_REF] Kamau | Linkage disequilibrium and recombination rate estimates in the self-incompatibility region of Arabidopsis lyrata[END_REF][START_REF] Ruggiero | Hitch-hiking to a locus under balancing selection: high sequence diversity and low population subdivision at the S-locus genomic region in Arabidopsis halleri[END_REF][START_REF] Roux | Recent and ancient signature of balancing selection around the S-locus in Arabidopsis halleri and A. lyrata[END_REF]. These studies sequenced short fragments of a small subset of the genes immediately flanking the S-locus, as well as more distant genes, and compared their polymorphism to that of a handful of "control" coding sequences from across the genome. In both species, the increase of polymorphism was limited to the genes immediately flanking the S-locus only, but the very sparse sampling of genes and the sequencing of small gene fragments only did not allow these previous studies to reach solid conclusions on the true genomic extent of this increase, and to precisely quantify the accumulation of deleterious mutations.

In this study, we combined whole genome sequencing data with a targeted resequencing approach to comprehensively analyze all genes and intergenic sequences within 75kb on either side of the S-locus in three sample sets each of A. halleri and A. lyrata. We compared the observed patterns of polymorphism in these regions with those of 100 unlinked randomly chosen regions used as genomic controls. The use of internal genomic controls provides a powerful way to factor out differences in demographic histories and/or sample structure. We consistently observed an increase of polymorphism within the first 25-kb region immediately flanking the S-locus only, with no detectable effect further along the chromosome. Contrary to predictions from classical models of sheltered genetic load, the putatively deleterious mutations that this narrow region carries do not segregate at higher population frequencies than the overall genomic background, and the relative rate of accumulation of non-synonymous to synonymous variants is also not elevated. These patterns are remarkably consistent across the different sequencing methods we employed and also across the different sample sets we studied in spite of differences in their specific demographic histories. Hence, our data suggest that linkage to one of the strongest known balanced polymorphisms does indeed result in elevated polymorphism, but is not associated with a detectable reduction of selection efficacy.

Results

Sequencing the S-locus flanking regions and control regions in large sample sets

To evaluate the impact of balancing selection on the genomic regions flanking the S-locus, we focused on the region where previous studies indicated that the signature of balancing selection was most likely encompassed; i.e. over a maximum of 75 kb on each side of the S-locus [START_REF] Kamau | Balancing selection and low recombination affect diversity near the self-incompatibility loci of the plant Arabidopsis lyrata[END_REF][START_REF] Kamau | Linkage disequilibrium and recombination rate estimates in the self-incompatibility region of Arabidopsis lyrata[END_REF][START_REF] Ruggiero | Hitch-hiking to a locus under balancing selection: high sequence diversity and low population subdivision at the S-locus genomic region in Arabidopsis halleri[END_REF][START_REF] Roux | Recent and ancient signature of balancing selection around the S-locus in Arabidopsis halleri and A. lyrata[END_REF]. We divided this region in three consecutive non-overlapping windows of 25kb (-25, -50 and -75 kb one one side and +25, +50 and +75 kb on the other side; Fig. 1). Together, the two 25kb windows closest to the S-locus contain a total of 11 annotated genes in the A. lyrata genome, the next upstream and downstream 25-50kb windows together contain 9 genes, and the most distant 25kb regions contain 13 genes (Hu et al. 2011, Fig 1). To compare these regions to the background level of nucleotide polymorphism, we also included in the analysis one hundred 25kb "control" regions unlinked to the S-locus. These control regions were randomly chosen across the A. halleri genome and selected to closely match the density of protein-coding sequences and transposable elements found at the S-locus flanking regions (proportion of CDS within the interval = 0.23 +/-0.0023; proportion of TEs = 0.28 +/-0.0028, Fig. 1). Because the extreme level of sequence divergence of the non-recombining interval containing the S-locus itself precludes mapping of short reads among S-haplotypes [START_REF] Goubet | Contrasted patterns of molecular evolution in dominant and recessive self-incompatibility haplotypes in Arabidopsis[END_REF], we excluded this region from further analysis and focused on the flanking regions only (Fig. 1).

Figure 1: Schematic representation of the S-locus and its flanking regions.

Protein-coding genes in the flanking regions are represented as filled blue rectangles. The genomic regions studied are distributed between positions 9,264,458 and 9,451,731 along chromosome 7 of the A. lyrata genome assembly [START_REF] Hu | The Arabidopsis lyrata genome sequence and the basis of rapid genome size change[END_REF]. The S-locus (in grey) contains two protein-coding genes only, SRK and SCR (white rectangles) and is flanked by the ARK3 and Ubox genes (in the 5 ' and 3' directions, respectively). The S-locus region itself was not analysed in the present study. The percentage of CDS in each 25kb window is given on top of the figure.

To provide a comprehensive picture of the indirect effects of balancing selection, we analysed two closely related species that share the same orthologous SI system and show extensive trans-specific polymorphism at the S-locus, A. halleri and A. lyrata [START_REF] Castric | Repeated adaptive introgression at a gene under multiallelic balancing selection[END_REF]. In order to evaluate the robustness of our conclusions to different demographic histories, we analysed nucleotide sequence polymorphism data from two natural populations of A. halleri (Nivelle, n=25 and Mortagne, n=27 that have been recently introduced in the North of France in association with industrial activities) and two natural populations of A. lyrata (Plech, n=18 from the core of the species range and Spiterstulen, n=23 from the edge of the species range, Table S1). To evaluate the robustness of our conclusions to different sampling strategies, we also included samples from more extended geographic regions of A. lyrata (North America, n=27 distributed across three distinct populations) and A. halleri ssp. gemmifera (Japan, n=47 distributed across six distinct populations). For the Nivelle, Mortagne and North American samples, we developed a dedicated sequence capture protocol specifically targeting the control and S-locus flanking regions. For the Japan, Plech and Spiterstulen samples, we took advantage of published whole-genome resequencing datasets, but analysed only polymorphism of the regions included in the capture protocol.

We obtained an average of 59 million reads mapped for the samples sequenced by sequence capture and 1,310 million reads for the WGS samples. After stringent filtering, we were able to interrogate with confidence an average of 960,368 positions in control regions, 28,432 of which were variable and biallelic (3%). As expected, the number of variable sites across the control regions differed among sample sets, reflecting their different demographic histories (Table 1). The A. halleri Japan sample set was the least polymorphic of all, with observed heterozygosity Ho=0.00096, nucleotide polymorphism π=0.00128 and a proportion of polymorphic sites equal to 0.0088. At the other extreme, the A. lyrata Plech population was the most polymorphic, with Ho=0.00646, π=0.00758 and a proportion of polymorphic sites equal to 0.0299. These estimations of the background level of nucleotide polymorphism in each sample set were used as internal genomic controls for the study of polymorphism in S-flanking regions, where we were able to interrogate an average of 74,866 sites, containing 3,225 variable biallelic positions (4%).

Detection of the footprints of ancient balancing selection on the S-linked regions

Based on these comprehensive polymorphism data, we combined different approaches to characterise the impact of balancing selection. As a first step, we excluded the potential confounding effect that would arise if mutation rates were higher in the S-flanking regions than in the control regions. Comparison of the mean levels of divergence between A. thaliana and A. lyrata reference genomes showed no evidence for increased divergence in the windows flanking the S-locus, as would be expected if they tended to accumulate more mutations per unit time. Instead, these two regions tended to indicate a slight reduction rather than an increase of divergence, albeit not a significant one (Fig. S1).

Next, we followed a multilocus Hudson-Kreitman-Agade (HKA) approach to compare nucleotide polymorphism within A. lyrata or A. halleri sample sets, taking into account divergence from the outgroup A. thaliana between the 33 S-flanking genes and 67 randomly chosen control genes. The multilocus HKA test showed a highly significant departure from neutral expectation (mean X²=821, P=0, df=33; Table S3), indicating that the two categories of loci differed in their relative patterns of polymorphism. The mean estimate of the selection parameter for the 33 S-flanking genes was above one (mean k=1.46 ; Fig. 2, Table S4), indicating higher polymorphism of the S-flanking genes compared with the control loci. k also tended to increase toward the S-locus, although the magnitude of this pattern varied across samples from different regions and differed between the 5' and 3' flanking regions.

Figure 2: Variation of the mean selection parameter (k) obtained for genes in the S-flanking regions. The thick solid line represents the mean value of k obtained across the six sample sets. The thin solid lines represent the A. lyrata sample sets (square=Plech, circle=Spiterstulen, triangle=North America). The dashed lines represent the A. halleri sample sets (open square=Japan, open circle=Nivelle, open triangle=Mortagne). The threshold value of 1 (no selection) is represented by the horizontal black line.

We then used the new powerful approach of [START_REF] Cheng | Flexible mixture model approaches that accommodate footprint size variability for robust detection of balancing selection[END_REF] that is robust to demographic variations to detect distortions of the site frequency spectrum along the chromosomal fragments and determine the maximum likelihood position of putative targets of balancing selection. We found strong signals of balancing selection in some of the control regions, specifically on chromosomes 3 and 4, but in most cases they were not consistent across all sample sets (Fig. 3). In contrast, the S-locus flanking regions carried consistent signals of balancing selection that were detected across all sample sets, even though they were not always the most extreme (Fig. 3). The exact position of the peak detected in the 25kb windows around the S-locus varied between sample sets (Fig. S2). Overall, these results provide evidence for a strong and consistent footprint of balancing selection on the regions flanking the S-locus.

Figure 3: Manhattan plots for signals of balancing selection (B 2,MAF scores) in the 100 control regions and the S-locus region (red dots and arrows) along the A. lyrata genome. Chromosomes 1-8 of the A. lyrata genome, with the 100 control regions distributed on successive chromosomes represented by an alternation of grey and black dots. The horizontal grey solid line represents the median B 2,MAF scores across SNPs in the genome, and the black horizontal dashed lines represent the top 5, 2.5 and 1% percentiles.

The increased polymorphism of the S-locus flanking regions is mostly caused by an increase of the proportion of polymorphic sites

Then we sought to describe in detail how polymorphism of the S-locus flanking regions compared with the genomic background. To do so, we compared the values of several summary statistics of polymorphism from the S-locus flanking regions to their distribution across the 100 control regions. Specifically, we compared the nucleotide polymorphism (π), the observed heterozygosity (Ho), the mean frequency of the minor allele (MAF) and the proportion of polymorphic sites (number of observed polymorphic sites divided by the total number of sites considered). Significant excess of polymorphism statistics as compared to control regions was found for almost all sample sets in the two 25kb windows immediately flanking the S-locus for Ho (by a factor 1.7-fold in Plech to 6.4-fold in Japan, Fig. 4), π (by a factor 1.6-fold in Plech to 5.8-fold in Japan, Fig. 5), and the proportion of polymorphic sites (by a factor 1.6-fold in Plech to 3.9-fold in Japan, Fig. 6, Table S5). We observed only two exceptions to this pattern. For the +25kb window in Plech and the -25kb window in Spiterstulen, Ho was not significantly higher than in control regions. In stark contrast, the second and third consecutive 25kb windows on either sides of the S-locus generally showed no excess polymorphism as compared to control regions in any sample set, with the exception of the Spiterstulen population, where the -75kb and +50kb windows had a slightly higher proportion of polymorphic sites (Fig. 6).

To verify that the effect we observed was not specific to the particular window size we chose, we used Linear Models to test whether Ho, π and MAF of individual sites of the S-flanking regions (considered as response variables) declined when physical distance increased away from the S-locus. A highly significant negative effect of the distance to the S-locus was observed overall (Table S6), confirming the effect of proximity to the S-locus on polymorphism of sites in the flanking regions. An increase of polymorphism can be explained by 1) an increase of the frequency of mutations at polymorphic sites within each sampling region, 2) an increase of the proportion of polymorphic sites, or 3) a combination of both. Noting that the proportion of polymorphic sites increased in the S-locus flanking region compared with the control regions, we wanted to test whether the allele frequencies at those polymorphic loci were also affected. To do that, we reiterated the analyses above, but considering the polymorphic sites only. We found no difference with respect to control regions when computing the Ho, π and MAF statistics on polymorphic sites only, with a single exception for Ho, which showed a higher value in the +25kb window flanking the S-locus in the North American sample set from A. lyrata (Fig. S3). Hence, the higher polymorphism detected in the S-locus flanking region is essentially due to an increase of the proportion of polymorphic sites rather than to a shift in the allele frequency spectrum. This is confirmed by the absence of deviation of the Tajima's D statistic compared with control regions (Fig. 7). Overall, our results thus show elevated nucleotide polymorphism at the S-locus region, which is mostly caused by a larger number of polymorphic sites rather than by an increased frequency at which the polymorphic sites segregate.

Higher density of putative deleterious mutations in the S-locus flanking region

To determine if the indirect effect of balancing selection described above is associated with the accumulation of a "sheltered" genetic load, we examined the accumulation of 0-fold degenerate sites only, assuming that the majority of amino-acid polymorphisms are deleterious to some extent (Eyre-Walker and [START_REF] Keightley | The distribution of fitness effects of new mutations[END_REF]). Like for total polymorphism above, we observed an increase of polymorphism at 0-fold degenerate sites in the S-locus flanking regions as measured by Ho, π or MAF when compared to control regions (Fig. S4 and S5 for Ho and π respectively), which is mostly due to an increased proportion of polymorphic sites in the first 25kb surrounding the S-locus (Fig. 8, Table S7, Fig. S6). The magnitude of this increase as compared to the genomic background ranged from 1.92 to 3.66-fold across sample sets (Table S7). A GLM restricted to 0-fold sites confirmed the effect of proximity to the S-locus (Table S6). We further compared the ratio of π between 0-fold and 4-fold degenerate sites. If balancing selection decreased the efficacy of the purge of deleterious mutations, we expect an elevation of the π0-fold/π4-fold ratio in the S-locus flanking regions. However, we found no evidence for such an increase in the S-locus flanking regions as compared to the control regions (Fig. 9), with the exception of the -50kb window in the North American sample set of A. lyrata. Finally, to define more precisely the genes that are affected by the sheltered genetic load, we explored the functional annotations of the genes contained in the S-flanking regions (Table S9). The -25kb and +25kb regions, where we found an effect of linkage to the S-locus, contained only eleven annotated genes in the A. lyrata Four of these genes are receptor-like serine/threonine-protein kinases (AT4G21410, AT4G21400, AT4G21390, AT4G21380/ARK3), one is an ubiquitination protein (AT4G21350/Ubox), two are transcription factors (AT4G21340, AT4G21330), one a peptidase (AT4G21323), one a transmembrane protein (AT4G21310), one a tetratricopeptide repeat (TPR)-like superfamily protein (AT4G21300), and a last one is a subunit of Photosystem II (AT4G21280).

Discussion

Elevated polymorphism but no decreased efficacy of purifying selection in the linked region

The segregation of deleterious mutations linked to the S-locus is expected to have important consequences for mating system evolution, affecting the maintenance of SI itself (Porcher & Lande 2005, Gervais et al. 2014) and the conditions for diversification of SI lineages [START_REF] Uyenoyama | Genealogy-dependent variation in viability among self-incompatibility genotypes[END_REF]. Our analysis across several replicate geographic regions in two closely related species is the first comprehensive genomic study to reveal the extent of the S-locus associated sheltered genetic load in a plant genome. In line with theoretical predictions from [START_REF] Schierup | The effect of hitch-hiking on genes linked to a balanced polymorphism in a subdivided population[END_REF], we show that the genomic region directly adjacent to the S-locus presents consistent signals of linked balancing selection and that polymorphism is elevated as compared to the genomic background, by a factor up to 5.8. This observation is important because Contrary to this expectation, however, we found no evidence that the efficacy of purifying selection is decreased for the mutations linked to the S-locus. The linked regions do accumulate more mutations than the control regions, but they segregate at population frequencies that are indistinguishable from those at control genes, and the π 0-fold /π 4-fold ratio is also unchanged. This accumulation of mutations is not due to an increased mutation rate in these genes because divergence from A. thaliana rather tends to be reduced compared with the background. This latter observation is in line with the repeated introgression observed at the S-locus between A. halleri and A. lyrata [START_REF] Castric | Repeated adaptive introgression at a gene under multiallelic balancing selection[END_REF], causing divergence between the two species to be more recent for the S-alleles than for the genomic background. Hence, the main factor causing the elevated polymorphism seems to be the deeper coalescence time among allelic lineages, allowing the accumulation of both neutral and deleterious mutations. One possible explanation for the fact that the distribution of allele frequencies in the linked regions are indistinguishable from those at the control regions can be related to the model proposed by [START_REF] Takahata | A simple genealogical structure of strongly balanced allelic lines and trans-species evolution of polymorphism[END_REF]. This model shows that the genealogical relationships among distinct S-allele lineages under gametophytic SI are expected to be identical to those of neutral genes, except they are expanded by a scaling factor S . If the number of S-alleles in the studied sample sets was extremely large and every chromosome we sampled corresponded to a different S-allele lineage, then the only difference between the genealogies of sequences around the and those in the control regions would lie in the different time scales. Testing this hypothesis would require phasing the flanking sequences with each S-allele to obtain haplotypes. However, based on the published estimates of S-allele frequencies in four of the population samples studied here (the Nivelle population of A. halleri, [START_REF] Llaurens | Does frequency-dependent selection with complex dominance interactions accurately predict allelic frequencies at the self-incompatibility locus in Arabidopsis halleri ?[END_REF]; the Japanese A. halleri sample set, [START_REF] Genete | Genotyping and De Novo Discovery of Allelic Variants at the Brassicaceae Self-Incompatibility Locus from Short-Read Sequencing Data[END_REF]; and the Plech and Spiterstulen populations of A. lyrata, [START_REF] Takou | Maintenance of adaptive dynamics and no detectable load in a range-edge outcrossing plant population[END_REF], it is already clear that randomly sampled chromosomes would be very unlikely to systematically correspond to distinct S-allele lineages. More importantly, the [START_REF] Takahata | A simple genealogical structure of strongly balanced allelic lines and trans-species evolution of polymorphism[END_REF] model is based on a gametophytic SI, while the S-locus of the Brassicaceae functions as a sporophytic SI, such that dominance/recessivity interactions modulate the selective effect between [START_REF] Schierup | Evolutionary dynamics of sporophytic self-incompatibility alleles in plants[END_REF][START_REF] Billiard | A general model to explore complex dominance patterns in plant sporophytic self-incompatibility systems[END_REF]. Whereas in gametophytic SI all S-alleles are expected to segregate at equal population frequencies and hence sharply depart from the neutral site frequency spectrum, in sporophytic SI, the recessive S-alleles are driven to high population frequencies, whereas dominant S-alleles remain relatively rarer. This asymmetry may diminish the contrast between the S-locus and the genomic background. Finally, in sporophytic SI the recessive S-alleles can form homozygous combinations, allowing some purging of deleterious variants, which could help to explain why we detected no apparent decrease of the efficacy of purifying selection. Developing new theoretical models taking into account the structure of the dominance hierarchy between S-alleles will be necessary to fully understand the effect of linkage to a SSI locus. Strikingly, our results are almost the mirror image of the pattern seen at the human MHC by [START_REF] Lenz | Excess of deleterious mutations around HLA genes reveals evolutionary cost of balancing selection[END_REF], where the elevated polymorphism of the genes in the linked genomic region is not due to deleterious mutations being more abundant, but to the fact that each of them tends to segregate at higher population frequency than the genomic background. For the MHC, balancing selection is believed to be driven by pathogen mediated selection (although the exact mechanism remains controversial ; see Spurgin and Richardson 2010), which is very different from the negative frequency dependent selection maintaining diversity at the S-locus. Dominance interactions between the balanced allelic lineages are also not expected at the MHC. It is currently not clear which specific feature of the balancing selection mechanism acting at the S-locus and at the MHC causes these sharply distinct genomic signatures.

A limited extent of the footprint along the chromosome

Because SI is arguably one of the most intense forms of long-term balancing selection, it could a priori be considered a favourable case to detect the footprints of its genomic signature. Yet, a salient feature of our results is the limited extent of the genomic region over which the effect of linkage to the S-locus can be detected, essentially spanning over the immediately 25kb flanking regions only. As compared to other strongly balanced polymorphisms such as sex-determining regions of sex-chromosomes or mating-type loci in fungi, which typically occupy large chromosomal portions, the S-locus itself occupies only 30-110 kb in A. halleri and A. lyrata, and only includes the genes involved directly in the SI recognition machinery [START_REF] Guo | Evolution of the S-locus region in Arabidopsis relatives[END_REF][START_REF] Goubet | Contrasted patterns of molecular evolution in dominant and recessive self-incompatibility haplotypes in Arabidopsis[END_REF]. The large chromosomal regions associated with sex-determining regions of sex-chromosomes are believed to result from the progressive extension of successive inversions that can ultimately capture a large number of genes, eventually expanding across most of the length of a chromosome (Charlesworth et al. 2005 ;[START_REF] Pannell | About PAR: The distinct evolutionary dynamics of the pseudoautosomal region[END_REF]. The classical models for this process entailed sexual antagonism, whereby the inversions selectively fix mutations that are favorable in one sex in the appropriate genetic combination, and the accumulation of deleterious mutations follows from the action of Müller's ratchet once recombination has ceased [START_REF] Rice | The Accumulation of sexually antagonistic genes as a selective agent promoting the evolution of reduced recombination between primitive sex chromosomes[END_REF]. However, recent models have shown that the successive fixation of inversions can still take place even in the absence of sexual antagonism, as a result of the effective masking of recessive deleterious mutations accumulated in the flanking regions of the sex-determining loci (Jay et al. 2021, Lenormand and[START_REF] Lenormand | Y recombination arrest and degeneration in the absence of sexual dimorphism[END_REF]. The reason why the S-locus flanking regions do not undergo this process of inversion, recombination arrest and degeneration in spite of the strong balancing selection it experiences, may be due to the limited size of the region upon which polymorphism is affected and the lack of effective sheltering, preventing the region from extending further efficiently. It would be interesting to investigate genomic patterns of sheltered load in other SI systems, in particular in systems with a collaborative non-self recognition mechanism, such as observed in Solanaceae, as the size of the S-locus is much larger in those systems and apparently includes functional genes unrelated to SI [START_REF] Wu | Sequence analysis of the Petunia inflata S-locus region containing 17 S-Locus F-Box genes and the S-RNase gene involved in self-incompatibility[END_REF].

The peak of polymorphism is robust to sample heterogeneity

Overall, in spite of the different demographic histories, sampling structures and sequencing technologies used, we find qualitatively very similar results across species and sampling regions. This provides strong support for the idea that the contrasts we observed between the and the genomic background are robust to these factors. The difference of demographic histories between populations was expected to modify the levels of polymorphism for control and S-flanking regions (see e.g. [START_REF] Fijarczyk | Detecting balancing selection in genomes: limits and prospects[END_REF]. The A. lyrata Plech population presented the highest level of polymorphism in the genomic background, in line with this population being at the core of the range of European A. lyrata [START_REF] Takou | Maintenance of adaptive dynamics and no detectable load in a range-edge outcrossing plant population[END_REF]). The A. lyrata Spiterstulen population had lower polymorphism in the genomic background, in line with the strong reduction of effective population size it experienced during the colonisation of Norway within the past 100,000 years [START_REF] Takou | Maintenance of adaptive dynamics and no detectable load in a range-edge outcrossing plant population[END_REF]. The Nivelle and Mortagne populations of A. halleri have recently colonised the north of France during the last century from ancestral German populations [START_REF] Pauwels | Multiple origin of metallicolous populations of the pseudometallophyte Arabidopsis halleri (Brassicaceae) in central Europe: the cpDNA testimony[END_REF]. In spite of their composite origin (multiple populations), the North American A. lyrata and the Japanese A. halleri sample sets had the lowest polymorphism, possibly as the results of major demographic bottlenecks they experienced in the course of post glacial colonisations, at least in A. lyrata [START_REF] Clauss | Population genetic structure of Arabidopsis lyrata in Europe[END_REF][START_REF] Kawabe | Patterns of polymorphism and demographic history in natural populations of Arabidopsis lyrata[END_REF]). Although our strategy was not designed to interpret the quantitative differences observed among sample sets, we note that the relative elevation of polymorphism at the S-locus seems to be more pronounced in sample sets with lower baseline levels of diversity across the genome. This is consistent with the observation that the S-locus appears to be sensitive to demographic effects than the genomic background [START_REF] Takou | Maintenance of adaptive dynamics and no detectable load in a range-edge outcrossing plant population[END_REF].

The strongest increase of polymorphism at the S-locus as compared to the genomic background across all samples was found in the two sample sets composed of multiple populations (A. halleri from Japan and A. lyrata from North America). Thus, even though population stratification is to modify the site frequency spectrum, it did not prevent the detection of increased polymorphism in the flanking regions due to balancing selection at the S-locus. We note that [START_REF] Ruggiero | Hitch-hiking to a locus under balancing selection: high sequence diversity and low population subdivision at the S-locus genomic region in Arabidopsis halleri[END_REF] and [START_REF] Roux | Recent and ancient signature of balancing selection around the S-locus in Arabidopsis halleri and A. lyrata[END_REF] similarly used regional samples, and also detected an excess of polymorphism in genes flanking the S-locus, albeit with a much lower level of resolution. Some North American populations have experienced a loss of SI, and have shifted to partial selfing. Selfing is generally expected to reduce the effective rate of recombination, and might thus expand the footprint of balancing selection [START_REF] Ness | Genomic consequences of outcrossing and selfing in Plants[END_REF]). On the other hand, however, selfing may have been expected to relax the intensity of balancing selection on the S-locus. The populations considered here were specifically chosen because they are predominantly outcrossing [START_REF] Foxe | Reconstructing origins of loss of self-incompatibility and selfing in North american Arabidopsis lyrata: a population genetic context[END_REF]), so we expect this effect to be minor.

The nature and number of mutations causing the load

The absence of protein-coding genes within the region itself beyond those directly involved in the SI machinery [START_REF] Goubet | Contrasted patterns of molecular evolution in dominant and recessive self-incompatibility haplotypes in Arabidopsis[END_REF] suggests that mutations causing the sheltered load are likely to lie in the flanking region rather than in the S-locus region itself. While [START_REF] Stone | Sheltered load associated with S-alleles in Solanum carolinense[END_REF], [START_REF] Stift | Inbreeding depression in self-incompatible North-American Arabidopsis lyrata: disentangling genomic and S-locus-specific genetic load[END_REF] and [START_REF] Llaurens | The sheltered genetic load linked to the S locus in plants: new insights from theoretical and empirical approaches in sporophytic self-incompatibility[END_REF] showed that a sheltered load linked to the S-locus could be detected at the phenotypic level, our study provides the first genomic demonstration of an accumulation of potentially deleterious mutations in the S-linked regions. Identifying more precisely the mutations causing the sheltered load would still require fine mapping, but our work suggests that they are most likely to be found in very close proximity to the S-locus (<25kb). The phenotypic traits on which the load was documented varied across the different studies (seed dormancy in Papaver rhoeas, Lane and Lawrence 1995; seed survival in Solanum carolinense, Stone 2004; leaf development and juvenile survival in A. halleri, [START_REF] Llaurens | The sheltered genetic load linked to the S locus in plants: new insights from theoretical and empirical approaches in sporophytic self-incompatibility[END_REF]juvenile survival in Stift et al. 2013; horticultural traits in Rosa, [START_REF] Vieira | The identification of the Rosa S-locus and implications on the evolution of the Rosaceae gametophytic self-incompatibility systems[END_REF]), as would be expected given that the S-locus lies in different genomic environments in distant species, and given that the deleterious mutations are expected to hit the different flanking genes in a random manner. [START_REF] Kawabe | Comparative gene mapping in Arabidopsis lyrata chromosomes 6 and 7 and A. thaliana chromosome IV: evolutionary history, rearrangements and local recombination rates[END_REF] speculated that the low number of genes in the S-genomic region is probably not high enough for a large sheltered load to have an impact on fitness compared to the overall genomic load. Here we show that the genomic interval whose polymorphism is affected by linkage with the S-locus comprises eleven genes. Several of these genes were previously shown to be associated with deleterious phenotypic traits in A. thaliana, or to be ancestral genes shared with the sister family Cleomaceae (Ubox, ARK3 and AT4G21390, [START_REF] Cheng | The Tarenaya hassleriana Genome Provides Insight into Reproductive Trait and Genome Evolution of Crucifers[END_REF]. For instance, mutants of the transcription factor AT4G21330 exhibit abnormal anther morphology at the beginning of stage 4 [START_REF] Zhang | Regulation of Arabidopsis tapetum development and function by DYSFUNCTIONAL TAPETUM1 (DYT1) encoding a putative bHLH transcription factor[END_REF]) and the gene ARK3 has been shown to be implicated in root development [START_REF] Dwyer | A superfamily of S locus-related sequences in Arabidopsis: diverse structures and expression patterns[END_REF]). Hence, it is clear that some of these genes have important functions, and are obvious candidates for the future dissection of the genetic architecture of the mutation load sheltered by the S-locus.

A limitation of our population genetics approach is that it was designed to detect the collective accumulation of mutations rather than individual high-impact mutations. However, it possible that a low number of high-impact mutations, rather than a collection of small-effect mutations, are causing the load. Indeed, the selective dynamics of lethal mutations vs. slightly deleterious mutations can be sharply different (e.g. [START_REF] Clo | Hidden genetic variance contributes to increase the short-term adaptive potential of selfing populations[END_REF], and in the latter case finely dissecting the load at the genetic level will remain challenging. In addition, while our sequence capture approach also includes the intergenic sequences, we quantified the load based on coding sequences only. Previous studies demonstrated that polymorphism on intergenic regions could be under purifying selection [START_REF] Mattila | Impact of demography on linked selection in two outcrossing Brassicaceae species[END_REF] for an example in A. lyrata), so it is also possible that besides the coding sequences, mutations in intergenic sequences contribute to the load, hence making our estimation of the sheltered load an underestimate. Another limitation of our work is the focus on SNPs, while structural variants may also have strong deleterious effects. Long-read sequencing would now be required to achieve a more detailed analysis of these types of polymorphisms. A final limitation of our work is that theoretical models of the effect of SI on the flanking regions have assumed a gametophytic SI system [START_REF] Schierup | The effect of hitch-hiking on genes linked to a balanced polymorphism in a subdivided population[END_REF], while the SI system in Arabidopsis is sporophytic. An exciting next step will be to compare the number and identity of deleterious mutations associated with dominant vs. recessive alleles, a task that will require phasing polymorphisms and is beyond the scope of the present study.

Material and methods

Source plant material

We worked on natural accessions from two closely related species, A. halleri and A. lyrata, each represented by samples from three regions named Japan, Mortagne (France) and Nivelle (France) for A. halleri, and Plech (Germany), Spiterstulen (Norway) and North America for A. lyrata (Table S1). For the Japan, Spiterstulen and Plech samples, we used available whole genome sequencing (WGS) data obtained by [START_REF] Kubota | A genome scan for genes underlying microgeographic-scale local adaptation in a wild arabidopsis species[END_REF] and [START_REF] Takou | Maintenance of adaptive dynamics and no detectable load in a range-edge outcrossing plant population[END_REF]. The Japan sample set was composed of 47 individuals of A. halleri, subsp. gemmifera originating from six different populations (17 individuals from Fujiwara, 17 from Ibuki, 2 from Inotani, 3 from Itamuro, 4 from Minoo and 4 from Okunikkawa; [START_REF] Kubota | A genome scan for genes underlying microgeographic-scale local adaptation in a wild arabidopsis species[END_REF], the Spiterstulen (26 individuals) and Plech (23 individuals) sample sets were from single locations [START_REF] Takou | Maintenance of adaptive dynamics and no detectable load in a range-edge outcrossing plant population[END_REF]. For the three other sample sets, we collected individuals and developed a dedicated targeted enrichment capture approach to sequence the genomic regions of interest. The North american sample set of A. lyrata was composed of 26 individuals from three highly outcrossing populations from the Great Lakes region, named IND (Indiana Dunes National Lakeshore in Michigan, n=8), PIN (Pinery Provincial Park in Ontario, n=10) and TSS (Tobermory Provincial Park in Ontario, n=8) [START_REF] Foxe | Reconstructing origins of loss of self-incompatibility and selfing in North american Arabidopsis lyrata: a population genetic context[END_REF]. We collected 25 individuals from the Nivelle population (50°47'N, 3°47'E, France) and 27 individuals from the closely related Mortagne population (50°47'N, 3°47'E, France). In total, we complemented the 88 individuals with whole genome data with an additional 78 newly sampled individuals that we sequenced with the targeted sequence capture approach described below.

S-locus flanking regions and control regions

To evaluate the effect of balancing selection on the S-locus, we developed an original approach based on the comparison between the patterns of polymorphism of the two flanking regions on either side of the S-locus to those of a set of 100 randomly chosen control regions. The S-locus region can be poorly represented in whole-genome assemblies, so we first sequenced them using two A. halleri BAC clones that we newly obtained following the approach of Goubet et al. (2008) from a BAC library constructed from a mixture of several A. halleri individuals from Italy. These two BAC clones were chosen so as to cover entirely the 5' and 3' regions on either side of the S-locus (BAC clones 37G17 and 21E5, respectively ;10.6084/m9.figshare.16438908). We computed the proportion of coding sequences (CDS) and transposable elements (TEs) on the first 75kb sequences immediately flanking the S-locus on these two BAC clones (but excluding the non-recombining region within the S-locus itself), and we used these two statistics to select a set of matched control regions from across the A. halleri genome [START_REF] Legrand | Differential retention of transposable element-derived sequences in outcrossing Arabidopsis genomes[END_REF]). To do this, we used bedtools [START_REF] Hall | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF] to randomly select 25-kb contiguous genomic intervals. The genomic intervals were retained if their density of CDS and TEs closely matched that of the actual S-locus flanking regions (within 10%). If the proportions of CDS and/or TEs departed from those values, the region was discarded and a new region was picked until a total of 100 genomic intervals was included. The genomic coordinates of the control regions are given in Supplementary Table S9, and their sequences in fasta format are available at 10.6084/m9.figshare.16438908. Because the control regions were defined initially on the A. halleri reference, we used sequence similarity (using YASS, [START_REF] Noé | YASS: enhancing the sensitivity of DNA similarity search[END_REF] to identify orthologous regions along the A. lyrata genome. distribution in control regions with the S-flanking regions, and considered a significant difference when the mean B 2,MAF value was outside the 95% percentile of the distribution in control regions.

To control for a possible difference in mutation rates between genes in the S-locus flanking regions and genes in the control regions, we then compared their pattern of molecular divergence between A. lyrata and A. thaliana (TAIR10 genome) at the sites retained for the polymorphism analysis (i.e. having passed the coverage filter). We identified orthologs as the best hits in the A. thaliana genome using YASS, retaining alignments with a minimum e-value of 0.01 and an identity above 70%. Pairs of sequences were then aligned with clustalOmega (Sievers et al. 2011) and the proportion of divergent sites was determined using a custom Python script (https://github.com/leveveaudrey/analysis-of-polymorphism-S-locus).

We further compared the ratio of within-species polymorphism to between-species divergence [START_REF] Kreitman | A test of neutral molecular evolution based on nucleotide data[END_REF]) using the multilocus maximum likelihood HKA framework developed by [START_REF] Wright | The HKA test revisited: a maximum likelihood ratio test of the standard neutral model[END_REF] and available at https://github.com/rossibarra/MLHKA. The algorithm is currently limited to only one hundred genes, so we tested the 33 S-locus flanking genes and a randomly chosen subset of 67 control genes. Specifically, we compared a model with free mutation at each locus and no selection against a model with free mutation but where each of the 33 S-locus flanking genes are allowed to have their own selection coefficient (k). This parameter corresponds to the relative increase of polymorphism of the S-linked genes compared to genes in the control regions, taking into account differences in divergence between A. lyrata and A. thaliana across loci. We used a log-likelihood ratio test with 33 degrees of freedom to compare the likelihood of these two nested models. Chain length was set to 100,000 and separate analyses were performed for each sample set independently.

Decomposing the signals of balancing selection

We then decomposed the signal of balancing selection across the S-locus flanking regions into a series of elementary statistics. For each site, we estimated the observed heterozygosity (Ho) as the number of observed heterozygous genotypes divided by the number of individuals in the dataset, and the minor allele frequency (MAF). We calculated π at each position using the vcftools -site-pi option [START_REF] Danecek | The variant call format and VCFtools[END_REF]. When a position of the A. lyrata genome was covered but not polymorphic, the Ho, MAF and π statistics were set to 0. For each statistic, we binned SNPs flanking the S-locus into 25kb intervals and compared the distribution of the mean value obtained for sites within non-overlapping windows of 25kb in the flanking regions with the distribution of the mean obtained across the 100 control regions. Finally, we used Linear Models on all the samples cumulated to test for a linear correlation between the exact distance of each SNP to the S-locus along the chromosome and each of the polymorphism statistics listed above with the populations as a random effect in LM. Finally, deviation from neutrality was also tested using Tajima's D for each region of 25kb around the S-locus, for which an excess of intermediate frequency polymorphisms suggests the presence of balancing selection (positive values of D), using the vcftools option -TajimaD [START_REF] Danecek | The variant call format and VCFtools[END_REF].

Quantifying the sheltered load of deleterious mutations

To determine the extent to which the S-locus flanking regions accumulate deleterious mutations, we first reiterated the same analysis with the previous parameters (Ho, MAF, π), but for the 0-fold degenerate sites only (determined using the script NewAnnotateRef.py; [START_REF] Williamson | Evidence for widespread positive and negative selection in coding and conserved noncoding regions of Capsella grandiflora[END_REF]. We assumed that all nonsynonymous changes are deleterious. Because all mutations at 0-fold degenerate sites alter the sequence of the encoded protein, we assumed that these mutations are deleterious (neglecting the rare cases where balancing selection could favour amino acid changes). In contrast, mutations at the 4-fold degenerate sites never alter the encoded amino acid, so we used them as neutral references. For the sake of simplicity, we discarded mutations on the 2-or 3-fold degenerate sites. ,264,458 and 9,451,731 along chromosome 7 of the A. lyrata genome assembly [START_REF] Hu | The Arabidopsis lyrata genome sequence and the basis of rapid genome size change[END_REF]. The S-locus (in grey) contains two protein-coding genes only, SRK and SCR (white rectangles) and is flanked by the ARK3 and Ubox genes (in the 5 ' and 3' directions, respectively). The S-locus region itself was not analysed in the present study. The percentage of CDS in each 25kb window is given on top of the figure. Tables Table 1:

Figure legend

Figure 4 :

 4 Figure 4: Mean Ho around the S-locus and across the control regions from throughout the genome. Each barplot represents the mean value of Ho obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of Ho mean in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 5 :

 5 Figure 5: Mean π around the S-locus and across the control regions from throughout the genome. Each barplot represents the mean value of π obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of π mean in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 6 :

 6 Figure 6: Proportion of polymorphic sites around the S-locus and across the control regions from throughout the genome. Each barplot represents the proportion of polymorphic sites obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of the proportion of polymorphic sites in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 7 :

 7 Figure 7: Tajima's D around the S-locus and across the control regions from throughout the genome. Each barplot represents the mean of Tajima's D obtained in S-flanking windows of 25kb. The distributions (count) in the 100 control regions are represented by a histogram (right). The 97.5% and 2.5% percentiles of the distributions are represented by dashed lines. The median value of the distribution in control regions is represented by black lines.

Figure 8 :

 8 Figure 8: Proportion of polymorphic sites among 0-fold degenerate sites around the S-locus and across the control regions from throughout the genome. Each barplot represents the proportion of polymorphic sites obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of the proportion of polymorphic sites in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 9 :

 9 Figure 9: π 0-fold /π 4-fold ratio around the S-locus and across the control regions from throughout the genome. The bars represent the proportion of polymorphic sites obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of the proportion of polymorphic sites in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 1 :

 1 Figure 1: Schematic representation of the S-locus and its flanking regions.Protein-coding genes in the flanking regions are represented as filled blue rectangles. The genomic regions studied are distributed between positions 9,264,458 and 9,451,731 along chromosome 7 of the A. lyrata genome assembly[START_REF] Hu | The Arabidopsis lyrata genome sequence and the basis of rapid genome size change[END_REF]. The S-locus (in grey) contains two protein-coding genes only, SRK and SCR (white rectangles) and is flanked by the ARK3 and Ubox genes (in the 5' and 3' directions, respectively). The S-locus region itself was not analysed in the present study. The percentage of CDS in each 25kb window is given on top of the figure.

Figure 2 :

 2 Figure 2: Variation of the mean selection parameter (k) obtained for genes in the S-flanking regions. The solid large black line represents the mean value of k obtained across the six sample sets. The black solid lines represent the A. lyrata sample sets (square=Plech, circle=Spiterstulen, triangle=North America). The black dashed lines represent the A. halleri sample sets (open square=Japan, open circle=Nivelle, open triangle=Mortagne). The threshold value of 1 (no selection) is represented by the horizontal black line.

Figure 3 :

 3 Figure 3: Manhattan plots for signals of balancing selection (B 2,MAF scores) in the 100 control regions and the S-locus region (red dots and arrows) along the A. lyrata genome. Chromosomes 1-8 of the A. lyrata genome, with the 100 control regions distributed on successive chromosomes represented by an alternation of grey and black dots. The horizontal grey solid line represents the median B 2,MAF scores across SNPs in the genome, and the black horizontal dashed lines represent the top 5, 2.5 and 1% percentiles.

Figure 4 :

 4 Figure 4: Mean Ho around the S-locus and across the control regions from throughout the genome. Each barplot represents the mean value of Ho obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of Ho mean in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 5 :

 5 Figure 5: Mean π around the S-locus and across the control regions from throughout the genome. Each barplot represents the mean value of π obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of π mean in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 6 :

 6 Figure 6: Proportion of polymorphic sites around the S-locus and across the control regions from throughout the genome. Each barplot represents the proportion of polymorphic sites obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of the proportion of polymorphic sites in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 7 :

 7 Figure 7: Tajima's D around the S-locus and across the control regions from throughout the genome. Each barplot represents the mean of Tajima's D obtained in S-flanking windows of 25kb. The distributions (count) in the 100 control regions are represented by a histogram (right). The 97.5% and 2.5% percentiles of the distributions are represented by dashed lines. The median value of the distribution in control regions is represented by black lines.

Figure 8 :

 8 Figure 8: Proportion of polymorphic sites among 0-fold degenerate sites around the S-locus and across the control regions from throughout the genome. Each barplot represents the proportion of polymorphic sites obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of the proportion of polymorphic sites in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Figure 9 :

 9 Figure 9: π 0-fold /π 4-fold ratio around the S-locus and across the control regions from throughout the genome. The bars represent the proportion of polymorphic sites obtained in non-overlapping regions of 25kb around the S-locus. The distributions (count) of the proportion of polymorphic sites in the 100 control regions are represented by a vertical histogram on the right. The 95% percentile of the distributions is represented by dashed lines. The median value of the distribution in control regions is represented by black lines. *** = observed value above the 99% percentile of control regions, ** = observed value above the 97,5% of control regions, * = observed value above the 95% of control regions.

Variation of the median Ho, π, MAF and proportion of polymorphic sites in control regions in each dataset.

	Species	Sample	Sequencing	Number of	Number of	Ho	π	MAF	Proportion
		names	method a	population	positions				of
				s	considere	x10 -3	x10 -3	x10 -3	polymorphic
					d				sites
									x10 -3
	A.	Japan	WGS	6	953,242	0.96 1.28	0.9	8.8
	halleri
		Nivelle	Capture	1	1,037,607 3.59 4.08	2.79	13.1
		Mortagne	Capture	1	1,059,569 3.71 4.24	2.93	14.6
	A.	Plech	WGS	1	1,190,287 6.46 7.58	5.11	29.9
	lyrata
		Spiterstulen	WGS	1	1,017,504 5.47 5.33	3.74	19.2
		North	Capture	3	503,976	3.15 4.97	3.53	17.9
		America

a The nucleotide sequence polymorphism data obtained by whole genome sequencing (WGS) came from published datasets.

Acknowledgements

This work was supported by a grant from the France-Berkeley Fund (to VC and Rasmus Nielsen); the European Research Council (NOVEL project, grant number 648321); and the Agence Nationale de la Recherche (TE-MoMa project, grant number ANR-18-CE02-0020-01). AL thanks the ERC and the University of Lille for funding her PhD project. The authors thank the UMR 8199 LIGAN-MP Genomics platform (Lille, France), which belongs to the 'Federation de Recherche' 3508 Labex EGID (European Genomics Institute for Diabetes; ANR-10-LABX-46) and was supported by the ANR Equipex 2010 session (ANR-10-EQPX-07-01; 'LIGAN-MP'). The LIGAN-PM Genomics platform (Lille, France) is also supported by the FEDER and the Region des Hauts-de-France. We thank Stephen I. Wright, Rasmus Nielsen, Violaine Llaurens, Sylvain Glémin and Camille Roux for helpful discussions. The authors thank the Région Hauts-de-France, and the Ministère de l'Enseignement Supérieur et de la Recherche (CPER Climibio), and the European Fund for Regional Economic Development for their financial support for the molecular facilities in Lille. This work has been performed using infrastructure and technical support of the Plateforme Serre, cultures et terrains expérimentaux -Université de Lille for the greenhouse/field facilities. Supplementary data include Fasta and Bed files of A. halleri regions and probes used for the sequence capture available online in figshare database at 10.6084/m9.figshare.16438908. All original sequence data are available in the NCBI Short Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra) with accession codes: PRJNA744343. All scripts developed are available in Github (https://github.com/leveveaudrey/analysis-of-polymorphism-S-locus ).

Library preparation, sequence capture and sequencing

For the 78 newly sequenced individuals, we purified DNA from 15 mg of dried leaves of each sample with Chemagic beads (PerkinElmer) following [START_REF] Holtz | Genotyping by sequencing using specific allelic capture to build a high-density genetic map of Durum Wheat[END_REF], using the manufacturer's instructions but with an additional Agencourt AMPure beads (Beckman) purification. DNA was quantified by Qubit and 50 ng of DNA was fragmented mechanically with Bioruptor (Diagenode) to obtain fragments of around 300bp, which we verified using a BioAnalyzer (Agilent) with a DNA HS chip. We prepared indexed genomic libraries using the Nextflex Rapid DNA Seq kit V2.0 (PerkinElmer) using the manufacturer's instructions. Briefly, extremities of sequences were repaired and tailed, ligated with universal adaptors P5/P7 containing multiplexing unique dual index (PerkinElmer), and amplified by five cycles of PCR. We then selected fragments between 150 and 300pb with AMPures beads and pooled all libraries in equimolar proportions.

The pooled libraries then proceeded to a sequence capture protocol using the MyBaits v3 (Ann Arbor, Michigan, USA) approach. Briefly, 120bp RNA probes were designed by MyBaits and synthesised to the complete set of one hundred 25kb control regions as well as the 75kb regions flanking the S-locus on either side, with an average tiling density of 2X (a total of 48,127 probes). In addition to the S-locus flanking regions and the control region, the capture array also contained a set of additional probes that were not used in the frame of the present project but are detailed in Supplementary Information. The indexed genomic libraries were hybridised to the probes overnight at a temperature of 65°C, and were finally sequenced by Illumina MiSeq (300pb, paired-end) by the LIGAN-MP Genomics platform (Lille, France).

Read mapping and variant calling

Raw reads from sequence-capture or WGS datasets (see Supplementary table S1) were mapped onto the complete A. lyrata reference genome (V1.0.23, [START_REF] Hu | The Arabidopsis lyrata genome sequence and the basis of rapid genome size change[END_REF]) using Bowtie2 v2.4.1 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF]. File formats were then converted to BAM using samtools v1.3.1 [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF]) and duplicated reads were removed with the MarkDuplicates program of picard-tools v1.119 (http://broadinstitute.github.io/picard). These steps were performed by the custom Python script sequencing_genome_vcf.py available in https://github.com/leveveaudrey/analysis-of-polymorphism-S-locus. We retained only reads which to the S-locus flanking or control regions. For the sake of consistency, we followed the same procedure for samples sequenced by WGS. Biallelic SNPs in these regions were called using the Genome Analysis Toolkit v. 3.8 (GATK, DePristo et al. 2011) with the option GVCF and a quality score threshold of 60 using vcftool v0.1.15 [START_REF] Danecek | The variant call format and VCFtools[END_REF]. For each sample independently, we computed the distribution of coverage depth across control regions using samtools depth [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF]. We excluded sites with either less than 15 reads aligned or coverage depth above the 97.5 % percentile, as the latter are likely to correspond to repeated sequences (e.g. transposable elements or paralogs). Sites covered by at least 15 reads but containing no SNP were considered as monomorphic. To exclude the possibility of spurious heterozygosity induced by the presence of paralogs, we verified that no SNP was systematically heterozygous across all individuals. The final number of sites in each sample set is summarised in Tables 1 and2.

Footprints of balancing selection

For each sample set, we first evaluated the distribution of the B 2,MAF statistic across all SNPs, which was designed to capture the distortion of the site frequency spectrum along chromosomes caused by linkage to a site under balancing selection [START_REF] Cheng | Flexible mixture model approaches that accommodate footprint size variability for robust detection of balancing selection[END_REF]. We then compared the B 2,MAF


I I Lobanova
email: inna.lobanova@univ-grenoble-alpes.fr

Stéphane Despréaux

Stéphane Labbé

Carola Celada-Casero

Frenk Van Den Berg

Micromagnetic modelling of magnetostriction under uniaxial stress

Introduction

During the process of magnetization ferromagnetic materials undergo changes in their shape or dimensions due to the coupling of magnetic order and strain. In particular, as the magnetization goes up, the magnetostrictive strain increases with magnetic field before reaching a saturating value, λ [START_REF] Jules | On the effects of magnetism upon the dimensions of iron and steel bars[END_REF]. This effect causes energy loss due to frictional heating in susceptible ferromagnetic cores. The effect is also responsible for the low-pitched humming sound that can be heard coming from transformers, where oscillating AC currents produce a changing magnetic field. The reciprocal effect, the change of the magnetic susceptibility in response to an applied field of a material when subjected to a mechanical stress, is called the Villari effect. Thus the mechanical stresses strongly influence the magnetic properties of ferromagnetic steel, affecting especially the hysteresis [START_REF] Singh | Effect of mechanical stress on excess loss of electrical steel sheets[END_REF][START_REF] Aydin | Effect of multi-axial stress on iron losses of electrical steel sheets[END_REF][START_REF] Daniel | Multiscale strategy for the determination of magneto-elastic behaviour: discussion and application to Ni-Zn ferrites[END_REF][START_REF] Ouaddia | Piezomagnetic behavior: experimental observations and multiscale modeling[END_REF]. Therefore the development of appropriate simulation tools is of particular interest for industrials to optimize the material conception.

The micromagnetic theory [6] is a powerful tool to describe the non-linear behaviour of ferromagnetic media. Although the micromagnetic modelling of magnetization process and domain structure are now well mastered, the description of magnetostriction and associated magnetomechanic interactions remains unachieved or badly implemented in the literature [START_REF] Shu | Micromagnetic modeling of magnetostrictive materials under intrinsic stress[END_REF]. Indeed the modelling of the intrinsic stress due to a free strain like magnetostriction is rather complex and therefore the mechanical approach to this problem is relevant. In this paper, we will address the development of a 3D stress dependent micromagnetic model that takes into account the magnetomechanical coupling. When the magnetic materials are subjected to an applied stress, the induced magneto-elastic energy must be taken into account in the Gibbs free energy. The modelling is therefore achieved by introducing elastic and magneto-elastic interactions depending on stress and magnetization in the free energy expression [START_REF] Carbou | Global weak solutions for the Landau-Lifschitz equation with magnetostriction[END_REF][START_REF] Mballa-Mballa | Micromagnetic modeling of magneto-mechanical behavior[END_REF].

The variation of magnetostrictive strain is linked to the domain wall movement and magnetization reversal process, which are usually complex phenomena that depend on the microstructure of magnetic material [START_REF] Singh | Effect of mechanical stress on excess loss of electrical steel sheets[END_REF][START_REF] Aydin | Effect of multi-axial stress on iron losses of electrical steel sheets[END_REF][START_REF] Hubert | Magnetic domains: the analysis of magnetic microstructures[END_REF]. In the present simulations, we have considered the effect of uniaxial stress on the magnetic hysteresis in ferrite with different grain orientations. This ferromagnetic material displays anisotropy of magnetic, magnetostrictive and elastic properties, and therefore the modelling of magnetic properties of different microstructures can help to shed more light on the influence of applied stress and predict its effect on the magnetic properties of the material.

Method

Micromagnetics and domain theory are based on the same variational principle which is derived from thermodynamic principles, as established initially by Landau and Lifshitz in [START_REF] Landau | On the theory of the dispersion of magnetic permeability in ferromagnetic bodies[END_REF] and reviewed in [6] and [START_REF] Döring | Handbuch der Physik: Mikromagnetismus[END_REF]. According to this principle, the vector field of magnetization directions m = M/M s is chosen so that the total free energy reaches an absolute or relative minimum under the constraint |m| = 1, where M is the magnetization vector and M s is the value of magnetization at saturation. In the initial micromagnetic model, the Landau-Lifshitz (LL) equation allows to predict the evolution of magnetization SMM25:374119 that depends on a combination of purely magnetic excitations by external magnetic, crystalline anisotropy, ferromagnetic exchange, and demagnetization fields [6]. In our model we have incorporated the internal elastic and magnetic stresses into the additional magnetostrictive excitation field and its contribution to the free energy [START_REF] Carbou | Global weak solutions for the Landau-Lifschitz equation with magnetostriction[END_REF]. The sum of all energies corresponding to these fields gives:

E tot (M) = E ex + E ani + E Zeeman + E demag + E m-el . (1)
The exchange energy term favours a homogeneous magnetic state such that any gradient of the local magnetization leads to an increase in E ex ,

E ex = A V (∇M) 2 dV, ( 2
)
where A is the exchange stiffness constant.

The energy of a ferromagnet that depends on the direction of the magnetization relative to the structural axes of the material is expressed as

E ani = V Φ(M)dV, (3)
where Φ(M) is the anisotropy energy density whose ex-

pression verifies Φ(M) = K 1 [(a 1 • M) 2 .(a 2 • M) 2 + (a 1 • M) 2 .(a 3 • M) 2 + (a 2 • M) 2 .(a 3 • M) 2 ] + K 2 (a 1 • M) 2 .(a 2 • M) 2 (a 3 •M) 2
in case of a cubic anisotropy. Here K 1 , K 2 are the anisotropy constants, and a 1 , a 2 , a 3 are the anisotropy directions.

The interaction energy of the magnetization vector field with an external field H ext (so called Zeeman energy) is simply:

E Zeeman = -µ 0 V H ext • M dV, (4)
where H ext is the external magnetic field. The stray field magnetic energy tends to demagnetize the material and is expressed as

E d = - µ 0 2 V H d • M dV, (5)
where H d is the demagnetizing magnetic field. And finally, the magneto-elastic energy is the sum of external and intrinsic stress contributions that depend on the magnetization M and the mechanical displacement field u:

E m-el = V λ m : [C : ε(u)-C : (λ m : M⊗M)].MdV, ( 6
)
where C is the elastic constants tensor, ε(u) is the mechanical deformation as a function of the displacement field u, λ m is the magnetostriction constants tensor.

The basic problem consists in minimizing the total free energy (1) and finding the equilibrium magnetization at each value of the external magnetic field, which allows to derive the magnetic hysteretic behaviour [6]. The computation of the equilibrium magnetization state is based on solving the LL equation over time. The discretized time step is chosen in a way that minimizes the magnetic energy. The displacement field is computed by solving the mechanical equilibrium equation ∇(C : u) = ∇(C : λ m : M ⊗ M ) with given material parameters, external magnetic field and equilibrium magnetization [START_REF] Shu | Micromagnetic modeling of magnetostrictive materials under intrinsic stress[END_REF][START_REF] Mballa-Mballa | Micromagnetic modeling of magneto-mechanical behavior[END_REF].

The micromagnetic model described above has been implemented in a self-written software EMicroM by using specific numerical methods [START_REF] Labbe | Numerical simulation of high-frequency behaviour of ferromagnetic materials[END_REF][START_REF] Labbe | Fast solution of quasi-stationary maxwell equations: multi-level toeplitz matrices. application to micromagnerism[END_REF][START_REF] Labbe | Fast computation for large magnetostatic systems adapted for micromagnetism[END_REF]. The problem is discretized using a Finite Volume Approximation which preserves the main properties of the operator given in the continuous model, i.e. positivity and symmetry. Furthermore, the discretization of the demagnetized operator leads to a block-Toeplitz structure for regular 3D cubic grid, which allows to build fast solving method for cubic domains. The representative volume elements (RVE) in this study are cubes 100 × 100 × 100 µm 3 with the resolution restricted to 64 × 64 × 64 (∼ 260 000) cubic cells, owing to the optimization of computing time. We performed computations on artificial 3D microstructures generated using multi-level Voronoi software developed by Tata Steel [START_REF] Kok | A microstructure based model for the mechanical behavior of multiphase steels[END_REF]. The material configuration files contain the information about Euler angles, that allows us to perform the rotation of crystals in reference to the RVE coordinate system. The material parameters used for simulations are given in the Table 1, where M s is the magnetization at saturation, A is the exchange stiffness constant, K 1 , K 2 are the cubic magneto-crystalline anisotropy constants, λ 100 , λ 111 are the magnetostriction constants along <100> and <111> We have computed hysteresis curves for selected microstructures: non-oriented (random) and oriented grains (η-and γ-fibers). There are two random structures: with 8 grains (grain size ∼ 50µm) and 1000 grains (grain size ∼ 13µm). The oriented microstructures both have 1910 grains (grain size ∼ 10µm) but the grains are rotated differently: η-fiber is characterised by <100> || RD and γfiber is characterised by <111> || ND.

Materials

Results and Discussion

The computed hysteresis curves under compressive and tensile uniaxial stresses for randomly oriented grains are displayed in the Fig. For 8 grains configuration the small number of grain orientation leads to a situation, where some of the grains has an easy-axis close to the applied field, and the tensile stress increases the anisotropy to favour the magnetisation vectors alignment with the external magnetic field. Therefore we observe that an increase of tensile stress leads to an increase of coercivity.

When we increase the number of grains to 1000, the dispersion of grain orientations becomes larger. The resulting hysteresis curve is given in the Fig. 2b. Here we do not observe any significant influence of the stress on coercivity, however we see an increase in permeability with the increasing tensile stress (blue and green curves in the Fig. 2b) and a decrease with the compressive uniaxial stress (black curve in Fig. 2b).

The results of hysteresis simulation for grain-oriented microstructures are given in the Fig. 3. We can see that the resulting hysteresis curves show a strong dependence of the properties on the texture. This confirms that our simulation software takes well into account the orientation of grains and the anisotropy of magnetic and elastic properties of the material.

For γ-fiber we observe a strong increase in coercivity with the compressive stress and an increase in permeability for the tensile stress that agrees well with predictions for this texture. For η-fiber increasing tensile stress leads to a decrease of coercive field, however it is still higher than in case of compression. This is due to the fact that for this virtual microstructure the <100> direction is strictly parallel to RD and therefore there is always one easy-axis parallel to RD, and two easy axes in a plain perpendicular to RD. In this situation there are multiple possibilities to minimize the energy, which have an impact on the computed results. This obstacle can be overcome by averaging the calculation result over multiple microstructures.

The presented model allows to calculate the average magnetostrictive strain in the material as an equivalent elastic strain. The results for the longitudinal magnetostrictive strain ε µ || as a function of the magnetisation M are presented in the Fig. 4.

The behaviour of the magnetostriction corresponds well to that observed experimentally [START_REF] Daniel | Multiscale strategy for the determination of magneto-elastic behaviour: discussion and application to Ni-Zn ferrites[END_REF][START_REF] Hubert | Multiscale magneto-elastic modeling of magnetic materials including isotropic second order stress effect[END_REF]. The sign of the strain at the saturation is in coherence with crystallographic orientations of grains in the considered microstructures. The tensile stress increases the longitudinal magnetostrictive strain ε µ || , and, in the case of the non-oriented grains, shifts the maximum to lower magnetization values, whereas the compressive stress decreases ε µ || . It is interesting that magnetostriction behaviour of simulated η-fiber shows similarity to the experimental measurements for DP steel. We also see the similarities in the hysteresis for tensile stress (see Fig. 5). Even though, the reason for this similarity is not evident because of the different microstructures of the 2 steel grades, this shows that the modelling of magnetostriction under uniaxial stress reflects well the hysteresis behaviour and is very sensitive to the microstructure.

Conclusions

The micromagnetic modelling of magnetostriction under uniaxial stress for different grain orientations has been carried out by including the magneto-elastic energy term in the Landau-Lifshitz equation. The computed magnetization curves display a change in the shape of hysteresis loops and the magnetic properties e.g., coercivity and remanence, under the applied stress in a manner which is in reasonable agreement with the reported experimental results of others. The magnetostriction as a function of magnetization follows well the characteristics of considered microstructures and applied uniaxial stress. The obtained results give a solid base for further development of the modelling algorithms as well as computational capacities.

Acknowledgements

The authors are grateful to Prof. O. Hubert from ENSPSay for provided experimental data of hysteresis and magnetostrictive strain in DP steel. The work described in the present paper has been developed within the project "Online Microstructure Analytics" (Ref. OMA, Grant Agreement No. 847296) that has received funding from the Research Fund for Coal and Steel of the European Union, which is gratefully acknowledged. This work was performed using HPC/AI resources from GENCI-IDRIS (Grant 2022-AD010913441) and HPC resources of Gricad. A part of the computations presented in this paper were performed using the GRICAD infrastructure, which is supported by Grenoble research communities.

Figure 1 :

 1 Figure 1: Top (from left to right): microstructure with 1000 grains, each grain has a different colour; crystallographic basic axis rotated in respect to RVE coordinate system for η-fiber and γ-fiber. Bottom: the displacement field distribution in case of tensile stress and compression.

M

  s = 1.7 M A/m A = 18 pJ/m K 1 = 48 kJ/m 3 K 2 = 0 λ 100 = 21 • 10 -6 λ 111 = -21 • 10 -6 C 11 = 228 GP a C 12 = 132 GP a C 44 = 116.5 GP a directions, and C 11 , C 12 , C 44 are the elastic stiffness tensor components in Voigt notation. The external magnetic field and the uniaxial stress are both applied parallel to RD direction. The example of microstructure and stress field distribution u are given in the Fig. 1.

2

 2

Figure 2 :

 2 Figure 2: Hysteresis curves under uniaxial stress for random orientation of grains a) 8 grains; insert: zoom on the region where the curves cross; b) 1000 grains. The applied stress values are given in legend.

Figure 3 :

 3 Figure 3: Hysteresis curves under uniaxial stress for a microstructure with 1910 grains (grain size ∼10 µm) and a) η-fiber ; b) γ-fiber. The applied stress values are given in legend.

Figure 4 :

 4 Figure 4: (a)-(d) Simulated longitudinal magnetostriction strain as a function of magnetization for a) randomly oriented 8 grains ; b) randomly oriented 1000 grains ; c) γ-fiber 1910 grains ; d) η-fiber 1910 grains. Solid lines correspond to the field sweep up, and dashed line -the field sweep down. (e)-(f) Experimentally measured longitudinal magnetostriction strain as a function of magnetization for e) DP steel ; f) ferrite NiZn (ahysteretic) [4]. Stress values are given in legends.

Figure 5 :

 5 Figure5: Comparison between the simulated hysteresis curves for η-fiber and the experimentally measured ones for DP steel[START_REF] Ouaddia | Piezomagnetic behavior: experimental observations and multiscale modeling[END_REF]

Table 1 :

 1 Material parameters