A C Galucio

J.-F Deü

F Dubois

The G 1 -scheme for Approximation of Fractional Derivatives: Application to the Dynamics of Dissipative Systems

Keywords: Fractional derivatives, viscoelasticity, multi-step scheme, linear dynamics

The Gear scheme is a three-level step algorithm, backward in time and second order accurate, for the approximation of classical time derivatives. In this article, the formal power of this scheme is used to approximate fractional derivative operators, in the context of finite difference methods. Numerical examples are presented and analyzed, in order to show the accuracy of the Gear scheme at the power 1 (G 1 -scheme) when compared to the classical Grünwald-Letnikov approximation. In particular, the combined G 1 -Newmark scheme is shown to be second-order accurate for a fractional damped oscillator problem.

INTRODUCTION

The importance of fractional calculus for modeling viscoelastic materials has been recognized by the mechanical scientific community since the work of [START_REF] Bagley | Fractional calculus-a different approach to the analysis of viscoelastically damped structures[END_REF]. The numerical approximation of such systems has been intensively studied since the late 1980s [START_REF] Padovan | Computational algorithms for FE formulations involving fractional operators[END_REF]. Meanwhile, the numerical community has also becoming interested in the approximation of fractional derivatives. Consider, for example, the pioneering theoretical work of [START_REF] Lubich | Discretized fractional calculus[END_REF], and the state of the art as proposed by [START_REF] Diethelm | Algorithms for the fractional calculus: A selection of numerical methods[END_REF]. Most applications use the discrete convolution formula proposed by Grünwald and Letnikov. Another direction for recent research is autonomous systems, in the context of diffusive representations (Matignon and Montseny, eds., 19981 Trinks and Ruge, 20021 Yu a n and Agrawal, 2002).

In this work, we focus on the application of a numerical method based on the Gear scheme (called a G 1 -scheme herein) to the approximation of fractional derivatives in linear dynamics. Note that preliminary tests of convergence have been performed, and are described elsewhere [START_REF] Galucio | An adaptation of the Gear scheme for fractional derivatives[END_REF]. In what follows, we first describe the G 1 -scheme. Two example applications are then presented and analyzed. The first of these studies a harmonic oscillator with fractional damping, in order both to validate the method and to derive an order of convergence. The use of the G 1 -scheme is then extended to analysis of viscoelastic beams subjected to an applied time-dependent force.

THE G 1 -SCHEME

We now introduce the G 1 -operator, which is based in the Gear scheme, for use in the approximation of fractional derivatives

G 1 2 1 2t 1 1 3 2 2 1 3 I 3 4 3 3 3 4 1 3 43 3 5 2 4 1 (1)
where 2t, which is assumed to be fixed, is the time step.

Let u be a time dependent function, known only by its discretized values u n at each time t n ,wheren is a positive integer. The function u n is approximated by u4t n 5,wheret n 2 n2t. The 1-derivative of u at time t n can be approximated using

4G 1 u5 n 2 1 2t 1 1 3 2 2 1 5 5 j20 g j41 u n3 j (2)
where g is a rational number. The calculation of these G 1 -coefficients is a difficult task, because of the cumulative numerical errors. In order to overcome this difficulty, the method employed here requires analytical calculation of these coefficients (using Matlab symbolic toolbox, in our case). For illustrative purposes, the reader is referred to table 1, where the first ten G 1 -coefficients are presented for three values of 1: 1/3, 1/2, and 3/4.

THE FRACTIONAL DAMPED OSCILLATOR PROBLEM

Consider a fractional single degree of freedom system subjected to a constant step load f for t 6 0 with zero initial conditions. The damping is taken into account by introducing a fractional damping term, or a spring-pot element in the formulation. The corresponding governing equation and initial conditions are

m 6 u 4 c7 1 1 1 u 4 ku 2 f8 t 6 0 u405 2 7 u405 2 0( 3 )
where m and k are mass and stiffness constants, respectively, and c7 1 is a fractional damping constant formed using relaxation time 7 and classical damping constant c. The aim of this section is to solve set of equations (3) with a direct (Newmark) time integration method in conjunction with an approximation of the 1-derivative 1 1 u (G 1 -o r GL-scheme). In order to validate these combinations, the approximate solution produced is compared to an exact solution considered in earlier work [START_REF] Galucio | An adaptation of the Gear scheme for fractional derivatives[END_REF]. Finally, error estimates in L 5 norm are performed. For a fixed time step 2t 2 192 m , this error is computed as

e m 5 2 max89 u4 j2t5 3 u j 9 8 j 2 08 8 2 m (4)
where m is a positive integer.

Algorithm

As mentioned above, the average acceleration algorithm is used to solve equation 3. The displacement history arising from the 1-derivative approximation (damping term) is shifted to the right-hand side of equation 3 [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]. Therefore, using equation 2, the governing equation can be written in its discretized form as

m 6 u n41 4 4k 4 5u n41 2 f n41 4 n41 (5)
where the non-classical terms and arise from the approximation of the 1-derivative:

2 c7 1 2t 1 1 3 2 2 1 (6a) n41 23 c7 1 2t 1 1 3 2 2 1 N 5 k21 g k41 u n413k (6b)
Note that the stiffness term is constant with time, depending only on the size of the time step, which is assumed to be fixed. The modified loading depends on the displacement history.

Results

In all calculations performed below, we assume that m 2 k 2 7 2 f 2 1 in a suitable unit system. The results shown in table 2, and in Figures 1 to 2, assume that c 2 1. Three different values of 1 are tested. The final time is chosen to be T 2 15, for various values of the time step 2t.

It should be remembered that the error estimates in the L 5 norm are obtained using an exact solution, based on formal power series, as mentioned above.

Table 2 shows the error estimates in the L 5 norm. It can be seen that the combined G 1 -Newmark scheme provides second-order accuracy for any value of 1. However, the use of a GL-Newmark algorithm decreases the order of accuracy to 1.

Figures 1 and2 show the evolution of the displacement for two values of 1,a n dt h e associated error estimates for the L 5 norm, for both the G 1 -Newmark and GL-Newmark schemes. In Figures 1(a) and 2(a), the exact solution of equation 3 is presented, along with its numerical approximations (using the GL and G 1 methods), with a time discretization corresponding to 2 6 2 64 time steps. It can be seen that the solution obtained using the G 1 -scheme is very close to the exact solution, while using the GL-method gives in an overestimated result.

Error estimates in the L 5 norms are presented in Figures 1(b) and 2(b). In both cases, the combined G 1 -Newmark scheme gives better accuracy than is achieved using GL-Newmark. The rates of convergence presented in table 2 are computed with 7-9 meshes, otherwise the slopes are wrongly estimated.

It should be emphasized that the order of the fractional derivative does not affect the rate of convergence (see table 2 and Figures 1 and2) when using a Newmark integrator. The influence of 1 is exerted in the mechanical behavior of the fractional damped oscillator by means of a damping factor. The result of this is that, when 1 decreases, the damping (and thus the time required to achieve the quasistatic time solution) increase.

In order to show the influence of added damping, the responses of the oscillator, computed using the G 1 -Newmark scheme, are shown in Figure 3 (a) for three values of the damping constant: c 2 05, c 2 10, and c 2 15 (see equation 3). These results, obtained for a semi-derivative problem, all use 1 2 1921 only the value of c varies. The corresponding error in the L 5 norm, for each result when compared to the exact solution, is shown in Figure 3 (b). The rate of convergence, as given in table 3 (see also table 2) remains the same for all three values of c. As in previous results, the results given by the combined G 1 -Newmark algorithm are of about the second order of accuracy.

EXTENSION TO VISCOELASTIC BEAMS

Viscoelastic Constitutive Equations

The one-dimensional fractional Zener model is adopted to describe the behavior of a viscoelastic material [START_REF] Bagley | Fractional calculus-a different approach to the analysis of viscoelastically damped structures[END_REF] in this section

4t5 4 7 1 1 1 4t5 2 E o 4t5 4 E 5 7 1 1 1 4t5 (7)
where and are the stress and strain, respectively, E o and E 5 are the relaxed and nonrelaxed elastic moduli, and 7 is the relaxation time.

In order to facilitate the numerical implementation of this model, we will introduce an internal variable as an "anelastic" strain function:

1 2 3 9E 5 (8)
This expression, when substituted into equation 7, gives an expression with only one fractional derivative operator.

Algorithm

For the sake of brevity, finite element considerations are not presented in this investigation.

For full details, the reader is referred to the description in a previous work [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]. The governing equation to be solved takes the form

6 7 8 M 6 q n41 4 4K 4 2 25q n41 2 F n41 4 1 n41 q 0 2 7 q 0 2 0 (9)
where M and K are the mass and stiffness matrices, q the degree-of-freedom vector, and F a mechanical load. 2 2 and 1 are the terms arising from the viscoelastic behavior of the beam such that

2 2 c 1 E 5 3 E o E o K (10a) 1 n41 23 c 1 E 5 E o K N 5 k21 g k41 q n413k 1 (10b)
where c 1 2 7 1 947 1 4 2t 1 5. Note that, as in the single degree of freedom problem, the added stiffness matrix 2 2 does not depend on time, while the dissipative force 1 depends on the history of "anelastic" displacements, which are updated at each time step as follows:

q n41 1 2 41 3 c 1 5 E 5 3 E o E 5 q n41 3 c 1 N 5 k21 g k41 q n413k
1 It is important to understand that the introduction of q 1 in the formulation does not imply an augmentation of the system. Instead, q 1 can be considered as an intermediate variable in the time scheme (see [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]. As in the previous example, the average acceleration algorithm is used to solve equation 9.

Results

Consider a viscoelastic cantilever beam of length L 2 150 mm, width b 2 25 mm, and thickness h 2 5 mm, discretized with 5 finite elements. The mechanical characteristics of the fictitious viscoelastic material are: Density 2 1000 kg/m 3 , Poisson's ratio 2 05, relaxed elastic modulus E o 2 1 MPa, non-relaxed elastic modulus E 5 2 50 MPa, relaxation time 7 2 1 ms, and order of the fractional derivative 1 2 05. The beam is subjected to a transverse load at its free end such that

F4t5 2 9 F o t9t 1 8 0 t t 1 F o 8 t t 1
where F o 2 001 N, t 1 2 50 ms, and T 2 1 s. The time step used is 2t 2 2ms,andthe whole time history of "anelastic" displacements is used in the calculations. Figure 4, shows the transient responses of the damped viscoelastic beam. The evolution of the tip displacement, and the phase-space diagram are plotted in Figure Figures 4(a) and (b), respectively. As expected, we observe that the oscillations of the viscoelastic beam are damped. It should be pointed out that these preliminary results show the versatility of the G 1 -scheme, since it is implementation is easy, and there are no additional costs relative to the GL-scheme.

CONCLUSIONS

A numerical method based on the Gear scheme for approximation of fractional derivatives is used here to model damping in linear dynamics. This G 1 -scheme is written in terms of a formal power series, the coefficients of which have to be calculated. The numerical evaluation of G 1 -coefficients is delicate, as a result of bad conditioning of the recurrence formula. However, with the help of formal calculus, cumulative numerical errors are avoided.

Two examples are presented and analyzed. In both cases, the average-acceleration algorithm is used to integrate the governing equation. The first example concerns a single degree of freedom oscillator, with fractional damping. In order to validate the presented approach, numerical results for the method proposed are compared to an exact solution for the single degree of freedom problem under a constant load [START_REF] Galucio | An adaptation of the Gear scheme for fractional derivatives[END_REF]. The second example deals with the finite element implementation of a viscoelastic beam subjected to a mechanical load. The combined G 1 -Newmark algorithm seems a promising tool for dynamic problems, as second-order accuracy is obtained in the solutions to both problems.

Figure 1 .

 1 Figure 1. (a) Exact and approximate solutions of equation (3) for 1 2 193 and 2t 2 T 92 6 . (b) Error estimates in the L 5 norm.

Figure 2 .

 2 Figure 2. (a) Exact and approximated solutions of equation (3) for 1 2 394 and 2t 2 T 92 6 . (b) Error estimates in the L 5 norm.

Figure 3 .

 3 Figure 3. (a) Exact and approximate solutions of (3) for 1 2 192 and 2t 2 T 92 9 . (b) Error estimates in the L 5 norm.

Figure 4 .

 4 Figure 4. Damped responses of the viscoelastic beam: (a) Tip displacement versus time1 and (b) phase-space diagram.

Table 1 .

 1 First ten coefficients g j41 of formal power series.

	j						1 2 193	1 2 192	1 2 394
	01	1	1
	1							3	4 9	3	2 3	31
	2							3	7 81	3	1 18	1 12
	3							3	104 2187	3	1 27	3	1 108
	4						3	643 19683	3	17 648	3	1 96
	5					3	4348 177147	3	19 972	3	7 864
	6				3	92809 4782969	3	59 3888	3	193 31104
	7			3	683552 43046721	3	71 5832	3	151 31104
	8			3	5164958 387420489	3	2807 279936	3	5813 1492992
	9	3	358288744 31381059609	3	10627 1259712	3	128713 40310784
	10	3	2805807422 282429536481	3	109159 15116544	3	430313 161243136

Table 2 .

 2 Rates of convergence, computed using the L 5 norm, for three values of 1

		1 2 193	1 2 192	1 2 394
	G 1	1.99	1.96	1.90
	GL	1.00	0.99	0.99

Table 3 .

 3 Rate of convergence computed using the L 5 norm for three values of c.

		c 2 050	c 2 100	c 2 150
	G 1	1.98	1.96	1.93
	GL	0.96	0.99	1.00


Matthieu Simon

Stéphanie Durand

Anthony Ricou

Nathalie Vrielynck

Baptiste Mayjonade

Jérôme Gouzy

Roxane Boyer

Fabrice Roux

Christine Camilleri

Françoise Budar

APOK3 , a pollen killer antidote in Arabidopsis thaliana

Keywords: transmission ratio distortion, pollen killer, poison-antidote, structural variation, Arabidopsis thaliana bioRxiv

à 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

Genetic loci that do not comply with Mendel's laws have been observed since the dawn of genetics. First considered as being genetic curiosities, these loci with transmission ratio distortion (TRD) are now recognized to be common in fungi, plants and animals, with a particularly high incidence in plants [START_REF] Fishman | Standard Deviations: The Biological Bases of Transmission Ratio Distortion[END_REF]. They are of major interest to understand genomic evolution, adaptation and speciation [START_REF] Presgraves | The molecular evolutionary basis of species formation[END_REF][START_REF] Lindholm | The Ecology and Evolutionary Dynamics of Meiotic Drive[END_REF] have been reported only once in A. thaliana [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF], and to our knowledge none has been molecularly deciphered so far in this species.

Meiotic drivers and gamete killers have retained particular attention for their potential in triggering genomic conflicts [START_REF] Lindholm | The Ecology and Evolutionary Dynamics of Meiotic Drive[END_REF][START_REF] Agren | Selfish genetic elements[END_REF]. The most intriguing feature of gamete killers is that killer alleles trigger a defect in the gametes that do not carry them. This can be explained by one of two main genetic models [START_REF] Nunez | Genetic Villains: Killer Meiotic Drivers[END_REF], schematically represented in Figure 1. In both models, all the genes responsible for gamete elimination are tightly linked and the killing factor is produced before (or during) meiosis, and still present in all the developing gametes. In the 'killer-target' model, a partner (the target), necessary for the killing activity, is encoded by the killed allele and expressed after meiosis, thus present only in gametes that carry it. This model applies to the SD system of D. melanogaster and the Sa locus of rice, for example [START_REF] Nunez | Genetic Villains: Killer Meiotic Drivers[END_REF]). In the 'poisonantidote' model, the killer allele also produces an antidote that counteracts the poison, but whereas the poison produced before meiosis subsists in the gametes, the antidote does not, so only the gametes that are able to produce it are rescued. Poison-antidote type segregation distorters have been described in fungi and plants and include the fission yeast wtf genes [START_REF] Nuckolls | wtf genes are prolific dual poisonantidote meiotic drivers[END_REF]) and the rice qHMS7 PK [START_REF] Yu | A selfish genetic element confers non-Mendelian inheritance in rice[END_REF]).

Here we decipher one of the PKs that were uncovered to contribute to a male sterility observed in hybrids between two A. thaliana natural variants, Shahdara (Sha) and Mr-0 [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF]. We show that this PK, located at the bottom of chromosome 3, belongs to the poisonantidote class, and we identified the antidote gene. This gene is chimeric and encodes a mitochondrial protein; it originated and evolved in the species A. thaliana within a highly variable locus. We precisely characterized the locus, showing that it contains at least two killer elements required for its activity, flanking the antidote gene. By exploring the diversity within the species, we found this PK in a number of A. thaliana hybrids. We then showed that in the killer genotypes, the locus strongly differs from those of neutral and killed genotypes by important structural variations, including duplications of the antidote gene. Lastly, we identified in the antidote several polymorphisms that correlate with its protective activity.

. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

The copyright holder for this this version posted September 10, 2021. ; https://doi.org/10.1101/2021.09.08.459474 doi: bioRxiv preprint

Results

Elements leading to a segregation bias at the L3 locus are common in natural variants of

A. thaliana

Amongst the PKs detected in the selfed progeny of Mr-0 x Sha hybrids [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF]), the L3 locus, located at the bottom of chromosome 3, induced a deficit in Sha homozygous progenies. This was linked to the death of pollen grains carrying the Sha allele at this locus. We also observed a strong bias at L3 against the Rak-2 allele in the Mr-0 x Rak-2 F2 population [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF], suggesting that some natural variants other than Sha possess at L3 alleles sensitive to the Mr-0 killer effect. Likewise, other variants than Mr-0 could carry killer alleles.

After crossing 30 accessions belonging to different diversity groups [START_REF] Simon | DNA fingerprinting and new tools for fine-scale discrimination of Arabidopsis thaliana accessions[END_REF] to Mr-0 and/or Sha, we analyzed the segregation at L3 in the progeny of the hybrids to reveal a possible killed or killer behaviour. Similar to Sha, 14 of the 26 accessions tested in crosses with Mr-0 for the killed status showed a bias in the progeny of their hybrid (Table 1). Genotype proportions were consistent with a 1:1 distribution of homozygotes Mr-0 and heterozygotes in most of the biased segregations, as expected if there is a gametophytic defect in these hybrids (Table1 Source Data 1). On the other hand, among 18 accessions analyzed for a killer behaviour, five induced a bias in the progeny of their hybrid with Sha, as did Mr-0 (Table 1), with a distribution of genotypes that was consistent with a gametophytic defect (Table1 Source Data 2). Among the 14 accessions that were tested for both killed and killer status, five accessions, including Col-0, had a neutral behaviour (neither killed nor killer). And, coherently, none of these 14 accessions was found to have both a killed and a killer behaviour (Table 1). natural accessions. (2) Table1-Source Data 2: Segregations at L3 in progenies of crosses of different natural accessions with Sha.

(3) Ita-0 was tested by crossing not with Mr-0 but with two other killer accessions, Ct-1 and Jea. nd: not determined.

Therefore, this segregation distorsion at L3 appears to be widespread in crosses between natural accessions of A. thaliana, and we focused on the Mr-0 x Sha cross to investigate the underlying genetic elements.

The segregation distorter at L3 in Mr-0 x Sha hybrids induces an allele specific impaired pollen development

We previously showed that a plant segregating Sha and Mr-0 alleles only at L3 while fixed Sha in the rest of its nuclear genome (hereafter ShaL3 H ) presented a strong bias against Sha homozygous progenies, which was linked to a deficit in pollen grains carrying the Sha allele [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF]. Here, we tested whether the segregation bias was dependent on the genetic background by comparing the selfed progenies of ShaL3 H with its equivalent in the Mr-0 nuclear background (MrL3 H ). Their progenies showed very similar deficits in Sha homozygotes (Table 2), indicating that the segregation distortion was independent of the fixed parental nuclear background. Accordingly, anthers from both genotypes showed similar proportions of dead pollen (Figure 2). (1) genotyped at marker M1 (described in Figure 5 Source Data 3). (2) Frequency of Sha homozygotes (expected frequency 0.25). ( 3) from [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF]. *** p < 0.001.

We further studied the ShaL3 H genotype, which flowers earlier than MrL3 H . We observed a strong distortion when ShaL3 H was used as male in a cross with Sha, but no distortion when it served as female parent (Table 3). These results confirmed that the segregation distortion was only due to a PK. (1) genotyped at marker M9 (described in Figure 5 Source Data 3). (2) Frequency of Sha homozygotes in the test-cross (expected frequency 0.5). *** p< 0.001; NS, not significant.

We used cytological approaches to determine whether the male dysfunction occurred during meiosis or during pollen development and to specify its timing. The male meiosis of ShaL3 H plants was identical to that of fixed siblings at all stages (Figure 3), excluding a meiotic defect as the source of the segregation bias in ShaL3 H plants. During male gametogenesis, abnormal pollen grains were observed in anthers of ShaL3 H plants from the bicellular stage, after the first 4). The proportion of abnormal pollen increased in the tricellular stage, and about 35% of mature pollen grains were dead in ShaL3 H plants (Figure 4), which reflects an incomplete penetrance since 50% dead pollen would be expected if the PK effect was total.

We concluded that, in these plants, the Sha allele at L3 is poorly transmitted because most of the Sha pollen grains fail to develop properly from the binucleate stage and eventually die.

The PK at L3 contains three genetic elements L3 was previously mapped in a 280 Kb interval at the bottom of chromosome 3 [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF]. We fine-mapped the PK in the genotype ShaL3 H using the presence of a segregation bias in the self-descent of recombinants as a robust phenotypic trait to narrow down the L3 interval: genetic markers fixed for Sha or Mr-0 alleles in recombinants with a significant bias in their progeny were excluded from the candidate interval. This strategy allowed us to map all the genetic elements necessary for the PK activity in an interval hereafter called PK3, corresponding to the region flanked by markers M5 and M13 (68 Kb in Col-0) (Figure 5A).

Out of a total of 4,717 plants genotyped, we found six recombinants between M5 and M13.

Recombination points of these recombinants, 27D6, 25A7, 52D12, 52D7, 8F10BH2 and 23G9 were finely localized. Only three of these plants, 52D7, 8F10BH2 and 23G9, were recombined between M6 and M12, which are very close to M5 and M13, respectively (Figure 5 Source Data 1). None of them presented a bias in its progeny whereas 27D6, 25A7 and 52D12, which are heterozygous between M6 and M12, did. Pollen viability of the recombinants, assessed by Alexander stainings, were consistent with the presence or absence of bias in their progenies (Figure 5B). Further information on the genetic structure of the PK3 was obtained from the three plants recombined between M6 and M12. First, by crossing 52D7 and 8F10BH2 with Sha, we converted their fixed portion of the interval into a heterozygous region (Figure 5A). The offsprings of these new genotypes (named i-52D7 and i-8F10BH2) did not show any segregation bias, whereas their siblings heterozygous along the whole interval did (Figure 5 Source Data 2). This indicated that the parts of the PK3 interval that were heterozygous in 52D7 and in 8F10BH2 both contained elements necessary for the PK activity. We thus segmented the PK3 locus into three genetic intervals, named PK3A, PK3B and PK3C (Figure 5A), PK3A and PK3C both carrying elements necessary for a functional PK. Then, in order to evaluate the PK3B interval, we crossed fixed progenies of 23G9 and 52D7 that inherited the recombination events from their parents, and obtain a plant (23G9#15 x 52D7#7) heterozygous at both PK3A and PK3C and fixed Mr-0 at PK3B (Figure 5A). The absence of segregation bias in the selfed progeny of this plant showed that PK3B also carried an element necessary for the PK activity.

Therefore, each of the three parts of the PK3 interval contains at least one element required for the PK activity. On one hand, Mr-0 alleles were required at PK3A and at PK3C: these two intervals thus contain killer elements. On the other hand, when the Sha allele was absent at PK3B while PK3A and PK3C were heterozygous, the PK was no longer active, indicating that either a target element from Sha was missing, or an antidote from Mr-0 was present in all the pollen grains produced.

The PK3 locus is highly variable

To highlight differences between Sha and Mr-0, we sequenced the entire locus in both accessions. The overall structure of the PK3 locus in Sha was very similar to that of Col-0 (Figure 6A), the main differences being the deletion of the transposable element (TE) AT3G62455, a 1308 bp insertion in the intron of AT3G62460, and an insertion of approximately 1 Kb in the intergenic region between AT3G62540 and AT3G62550. In contrast, the PK3 locus in Mr-0 locus was particularly complex as compared to Col-0 and Sha, showing many structural variations such as large deletions, insertions, duplications and inversions (Figure 6A). Two TEs present in Col-0, AT3G62455 and AT3G62520, were missing in Mr-0. The TEs AT3G62475, AT3G62480 and AT3G62490 that are located in the PK3A region of Sha, were absent from this region in Mr-0, but its PK3B region presented a large insertion of over 20 Kb that contained several TEs including AT3G62475, AT3G62480 and AT3G62490 homologues. Nonetheless, the same protein coding genes are present in the three genotypes, even though Mr-0 has two copies of the gene AT3G62510 and three copies of the genes AT3G62528, AT3G62530 and AT3G62540, with one copy of AT3G62530 and AT3G62540 being inserted into the second intron of AT3G62610 (Figure 6A).

Because the PK3 locus is highly rearranged between these three accessions, we looked at other variants of known status for the PK phenotype. The entire genomes of 10 such variants, four killers, three killed and three neutral, were de novo sequenced. PK3 sequence alignments revealed structural variations relative to Col-0 in all the killers (Figure 6B). On the contrary, PK3 loci of killed and neutral natural variants, excepted Bur-0, are mostly colinear with Col-0 (Figure 6B). When comparing the synteny of protein coding genes at the locus between A. thaliana accessions and A. lyrata (Figure 7), we observed that the locus structure in Ita-0 is similar to that of A. lyrata. This locus structure is also found in other Brassicaceae, such Boechera stricta (Figure 7). In contrast, Col-0, the two other neutral accessions Blh-1 and Oy-0, and Sha have a duplication of the A. lyrata gene AL5G45290, that encodes a pentatricopeptide repeat protein (PPR), resulting in two nearly identical genes (AT3G62470 and AT3G62540).

Five additional protein coding genes (AT3G62499 to AT3G62530) compared to the A. lyrata and Ita-0 sequences were found between the two PPR paralogues in these accessions. We also observed a great variability in the number of copies of all these genes according to the accessions. Compared to Col-0, Blh-1, Oy-0 and Sha, where they were present once, some of them were absent in the other killed accessions Cvi-0 and Are-10. On the contrary, in the killers, they have undergone a variable number of duplications, with one to four copies depending on both the genes and the accessions (Figure 7).

Given the complexity and diversity of the locus, no obvious candidate gene appeared for killer or target/antidote elements. For the killer elements, an additional difficulty results from the fact that two elements are necessary for the killing activity, one located in the PK3A interval and the other in the PK3C interval, thereby preventing to draw conclusions from the comparison of killer and non-killer alleles in a single interval. Moreover, the intervals PK3A, PK3B and PK3C were delimited by the mapping recombinants between Sha and Mr-0, and, due to the structural differences between the accessions, it is not possible to infer their frontiers in the other genotypes, at least for the limit between PK3A and PK3B. However, it can be noted that the locus of the neutral accession Bur-0 is nearly identical to that of the killer Jea: only the copy of AT3G62530 inserted into AT3G62610 in the PK3C interval in all the killers is missing in Bur-0 (Figure 7). This makes this copy, which is absent in all the non-killer alleles, a gene to be tested for killer activity, although it cannot be excluded that Bur-0 lost its killer activity due to another polymorphism in PK3C or in its unknown PK3A killer element.

We thus focused on the PK3B interval, which contains either a target element in killed alleles or an antidote in killer and neutral alleles. However, comparison of the PK3B sequences did not reveal a gene specific for the killed alleles that could encode a target, nor a gene specific for the killer and neutral alleles that could encode an antidote. We therefore exploited mutants in each of the genes present in this region.

AT3G62530, expressed in young developing pollen, encodes the antidote.

Col-0 had a neutral behaviour regarding the PK (Table 1). That means that, in a poison-antidote model, Col-0 would carry an antidote element whose inactivation should let Col-0 pollen unprotected against Mr-0 killer. A hybrid between Col-0 with an inactivated version of the antidote and Mr-0 would then present dead pollen, and produce an F2 with a segregation bias against the Col-0 allele at L3. In order to test the poison-antidote model, we thus used T-DNA insertion mutants available in this genetic background for the eight PK3B genes. None of the homozygous mutants (Col mut ) was distinguishable at the phenotypic level from Col-0 in our greenhouse conditions. We crossed each Col mut with an early-flowering Mr-0 genotype carrying a KO mutation in the FRIGIDA gene, hereafter named Mrfri (see Materials and Methods). All the Mrfri x Col mut F1 plants exhibited only viable pollen, except the hybrid with a T-DNA insertion in AT3G62530 (Col mut530 ) that presented aborted pollen (Figure 8). Then, we genotyped the F2 families at L3, which all followed the expected Mendelian proportions, with the exception of the Mrfri x Col mut530 and Mrfri x Col mut499 F2s. The latter presented a slightly biased segregation against the Col mut499 allele, but when the Mrfri x Col mut499 F1 was crossed as female or male parent using Mrfri as a tester, no mutant allele transmission bias was detected (Figure 8). In contrast, the Mrfri x Col mut530 F2 population showed a strong bias against Col mut530 homozygous genotypes, as those observed in hybrids with an active PK. We analysed the mutant allele transmission in the presence of a killer allele by crossing the Mrfri x Col mut530 F1 as female or male with Mrfri as a tester, and we observed a significant bias against the Col mut530 allele only when it came from pollen (Figure 8). At this step, we characterized the T-DNA insertion in AT3G62530, we checked that it really inactivated the gene and did not induce pollen abortion nor a bias in the transmission of the mutation to the progeny of the heterozygous Col mut530 mutant (Table 4). These results strongly suggested that AT3G62530 encoded an antidote. If this is the case, the mutation in AT3G62530 should have no effect on pollen viability nor on allele segregation in the absence of a killer allele. Indeed, no dead pollen was observed in the anthers of homozygous nor heterozygous Col mut530 (Figure 8 Supplement 2), and no segregation bias against the Col mut530 allele was detected in the self-progeny of the heterozygous Col mut530 nor in the progeny of the Sha x Col mut530 hybrid (Table 4). We concluded that the presence of a killer allele was necessary to trigger the death of Col mut530 pollen. This was confirmed by crossing Col mut530 with another killer accession, Ct-1: while no bias was found in the progeny of the Ct-1 x Col-0 hybrid, a bias against the Col mut530 allele was observed in the progeny of Ct-1 x Col mut530 (Table 4), which was very similar to the bias against the Sha allele observed in the progeny of the Ct-1 x Sha cross. (1) genotyped with marker M1. (2) genotyped with PCR primers used to characterise the mutant (Figure 8 Supplement 2). (3) The progeny of the F1 (Ct-1 x Col-0) presents no TRD (Table 4 Supplement 1). ( 4) frequency of Col mut530 homozygotes (expected frequency 0.25). *** p< 0.001; NS, not significant.

In the frame of the poison-antidote model, the antidote must be expressed in cells that need to be protected from the poison elements, in particular in developing pollen. Indeed, RT-PCR assays showed that AT3G62530, in addition to being expressed in leaves, is expressed in microspores (young developing pollen, before the first pollen mitosis, Figure 4) from plants either Col-0, Sha or Mr-0 at the locus (Figure 8). In addition, AT3G62530 seems one of the most expressed gene in microspores amongst those of the PK3B interval, in the three genotypes.

All together, the above results fit perfectly with a poison-antidote system where AT3G62530 codes the antidote, Col-0 and Mr-0 carrying functional forms of the antidote while Sha has a non-functional antidote allele. We named the gene APOK3, for ANTIDOTE OF POLLEN KILLER ON CHROMOSOME 3.

APOK3 encodes a chimeric protein addressed to mitochondria

Because APOK3 is expressed at similar levels in microspores carrying the Sha allele, not protected by an antidote, and microspores carrying Col-0 or Mr-0 alleles, which must harbor an active antidote (Figure 8), we hypothesized that the allelic differences for antidote activity are due to differences at the protein level.

In Col-0, APOK3 encodes a protein of 221 amino acids, belonging to the ARM-repeat superfamily (https://www.arabidopsis.org/), defined by the presence of tandem repeats generally forming alpha-helices (https://supfam.org/). Further analysis of structural domains identified three HEAT repeat domains (Figure 9A). Different parts of APOK3 are very similar to parts of proteins encoded by the genes AT3G62460 (98% identity on residues 1 to 44), AT3G43260 (66% identity and 73% similarity on residues 43 to 142) and AT3G58180 (57% identity and 78% similarity on residues 81 to 221) (Figure 9B). These three genes are located on the chromosome 3, and it is interesting to note that AT3G62460 is in the PK3A interval. AT3G43260 and AT3G58180 are both annotated as related to deoxyhypusine hydroxylases of other organisms, but a close examination showed that the genuine Col-0 deoxyhypusine hydroxylase is encoded by AT3G58180 (Figure 9 Supplement 1). The part of the AT3G62460 protein shared by APOK3 includes a mitochondria-targeting peptide (Figure 9A), suggesting that APOK3 is addressed to mitochondria. Indeed, APOK3 has been repeatedly found in A. thaliana mitochondrial proteomes [START_REF] Heazlewood | Experimental analysis of the Arabidopsis mitochondrial proteome highlights signaling and regulatory components, provides assessment of targeting prediction programs, and indicates plant-specific mitochondrial proteins[END_REF][START_REF] Klodmann | Defining the protein complex proteome of plant mitochondria[END_REF][START_REF] Taylor | The Arabidopsis thaliana 2-D gel mitochondrial proteome: Refining the value of reference maps for assessing protein abundance, contaminants and post-translational modifications[END_REF][START_REF] Konig | The mitochondrial lysine acetylome of Arabidopsis[END_REF][START_REF] Senkler | The mitochondrial complexome of Arabidopsis thaliana[END_REF]. We therefore hypothesized that, if the antidote acts in the mitochondria, the strength of the bias due to the PK could be influenced by the mother plant cytoplasmic background. As a first insight in this direction, we compared the segregation distortion in the progenies of plants heterozygous Sha/Mr-0 at PK3 differing only by their cytoplasmic backgrounds, and we observed that the bias was stronger in the Sha than in the Mr-0 cytoplasmic background (Figure 9C). We concluded that the antidote function of APOK3 is likely to be sensitive to variation in the mitochondrial genome.

APOK3 has undergone several duplication events within killer PK3 loci

Mr-0 has two strictly identical copies of APOK3 in the PK3B interval, which have the same structure as Col-0 and Sha genes. Mr-0 also has a third copy inserted with other sequences in the intron of AT3G62610 in the PK3C interval (Figure 10A), but the N-terminal part of this copy differs from that of the other two and is not homologous to AT3G62460, which suggests that this copy is functionally different from the others; it was thus named APOK3-like.

Our analysis of the natural variation at the PK3 locus revealed that all the killers analysed have multiple copies of APOK3 (in yellow on Figure 7), with Mr-0 and Ct-1 having two copies, and Cant-1, Jea and Shigu-2 having three copies each. In addition, these five accessions have an APOK3-like copy inserted in AT3G62610. Bur-0, which has a neutral behavior, also possesses three copies of APOK3, but no APOK3-like, even if it has an insertion of AT3G62540 in the intron of AT3G62610. The other neutral accessions and the killed ones have only one copy of APOK3, except Ita-0, in which no copy exists at the locus (Figure 7). No other copy of APOK3 was found elsewhere in the Ita-0 genome, neither by searching in the genomic sequence nor by PCR amplification in Ita-0 genomic DNA. Similarly, we did not find any gene encoding a protein closer to APOK3 than to AT3G58180 neither in A. lyrata nor in other Brassicaceae sequences available in the databases. Altogether, these results suggest that APOK3 is specific to A. thaliana and has evolved within the species.

Antidote and non-antidote forms of APOK3 differ by three amino acids

In order to further explore the sequence variation of APOK3 in relation with its antidote activity, we amplified and sequenced APOK3 in all the accessions whose behavior for the PK3 had been determined (Table 1). The APOK3 copies found in each killer accession whose genome was de novo sequenced were identical. In Etna-0, the only killer accession for which no genomic sequence was available, we obtained a unique Sanger sequence for APOK3, indicating that, whatever the number of copies it has, they are identical. Phylogenetic analysis clustered the sequences into three clades, one grouping all copies found in accessions possessing the antidote, i.e. killer and neutral. A second cluster groups most of copies from the killed accessions, with the exception of Are-1, Are-10 and Kas-2 which branch together as a third clade (Figure 10B).

We identified 11 different APOK3 haplotypes, five from non-killed accessions and six from killed accessions (Table 5). One polymorphism located 36 pb upstream of the ATG start codon and three non-synonymous SNPs in the coding sequence distinguished non-killed from killed alleles. Consequently, the APOK3 proteins with an antidote activity differ from those with no antidote activity by the three amino acid changes C85S, V101D and C105R (Table 5).

Interestingly, the last two of these amino acids are located in the first HEAT-repeat domain of the protein, suggesting these changes could modify the protein interactions and could make it functional or non-functional as an antidote.

F L S V A N S T D Y R G V non-antidote C S Y D S T C S V D C V L (1)
The names of the accessions whose entire PK3 sequences are known are in bold. The Col-0 sequence (TAIR 10) is used as reference. (2) Base 1 of the nucleotide sequence is the A of the start codon of the protein. The nucleotides in the exons are in upper case and those in the untranslated regions in lower case. Nucleotides identical to the reference sequence are represented as dots. Multiple sequence alignment was done with ClustalW (https://www.genome.jp/tools-bin/clustalw). SNPs and amino acids specific of antidote vs non-antidote forms are in bold.

Discussion

As pointed out by [START_REF] Burga | Toxin-Antidote Elements Across the Tree of Life[END_REF], the discovery of poison-antidote systems in eukaryotes is most often fortuitous, especially because nearly all these systems are species specific. In A.

thaliana too, we found PKs adventitiously when looking at a hybrid male sterility, of which PKs turned out to be unexpected components, in a cross between Sha and Mr-0 natural variants [START_REF] Gobron | A cryptic cytoplasmic male sterility unveils a possible gynodioecious past for Arabidopsis thaliana[END_REF][START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF]. In this study, we focused on one of these PKs, responsible for a deficit in Sha homozygotes at the bottom of chromosome 3 in Mr-0/Sha hybrid progenies, due to the death of pollen grains carrying this allele. A TRD at the same locus was found in the progenies of 19 out of 44 hybrids involving either Sha or Mr-0 as one of the parents, indicating that killed or killer alleles are not rare in A. thaliana (Table 1). This is in agreement with the report by [START_REF] Seymour | Transmission ratio distortion is frequent in Arabidopsis thaliana controlled crosses[END_REF] of frequent TRD in F2 populations from 80 A. thaliana founders, although these authors captured all types of segregation distorters, using a genome wide detection approach. Indeed, F2 distortion at one locus can result from allelic interactions either at the same locus or between different loci [START_REF] Fishman | Standard Deviations: The Biological Bases of Transmission Ratio Distortion[END_REF]. We cannot formally exclude that some of the biases observed in our F2 families are partially caused by other allelic interactions than the PK3 studied here. However, our subsequent results on sequence variation in the antidote gene APOK3, which grouped antidote and non-antidote forms of the genes in distinct clades (Figure 10B), support our conclusion that we indeed detected the effect of PK3 in these crosses. In addition, in biased F2 populations, most killer:Hz ratios complied with the 1:1 ratio that is expected for a bias due to a gametophytic defect controlled by a single locus. The exceptions concerned the accessions Kas-2, Kly-1, N13 and Cvi-0 in crosses with Mr-0, for which F2 genotypes presented an excess in heterozygotes (>55%). These cases could be explained by the presence, at loci partially linked to the PK3, of other segregation distorters that affect the transmission of the Mr-0 allele.

Our cytological analyses showed that, in ShaL3 H plants, a defect of pollen development is visible as soon as the bicellular stage (Figure 4). The proportion of affected pollen then increases during development but reaches only approximately 35% of dead pollen in mature anthers instead of 50% expected if the killer effect was total. This incomplete penetrance of the killer effect is supported by the fact that plants homozygous for the killed allele are found in the progeny of heterozygotes. Indeed, if considering that all the 35% dead pollen carry a Sha allele, we theoretically expect 11,5% of Sha homozygotes in the progeny, which is very close to the 10% observed (Table 2). Incomplete penetrance was also observed in other reported gamete killers and meiotic drivers, for example in tomato [START_REF] Rick | Abortion of male and female gametes in the tomato determined by allelic interaction[END_REF] or rice [START_REF] Matsubara | Relationship between transmission ratio distortion and genetic divergence in intraspecific rice crosses[END_REF]. In some cases, as in distorters from Drosophila (Larracuente and Presgraves 2012), rice [START_REF] Koide | Complex genetic nature of sex-independent transmission ratio distortion in Asian rice species: the involvement of unlinked modifiers and sex-specific mechanisms[END_REF]) and yellow monkeyflower [START_REF] Finseth | Selfish chromosomal drive shapes recent centromeric histone evolution in monkeyflowers[END_REF], the incomplete penetrance of the phenotype was due to the presence of unlinked modifier loci that modulate the strength of the distorter. In this regard, it is interesting to note that the strength of the biases that we observed in the F2s vary between crosses, either with Sha or with Mr-0 (Table 1 Source Data 1 & 2). It is conceivable that this variation be due to the presence, in some hybrids, of unlinked epistatic loci involved in the sterility caused by PK3, even if in the case of Mr-0/Sha we observed no difference in the biases between plants with Sha and Mr-0 fixed nuclear backgrounds (Table 2). In contrast, changing the Mr-0 cytoplasmic background to the Sha one had a strong effect on the elimination of the Sha allele in pollen (Figure 9C), making it almost absolute (only one Sha homozygous plant among 512 progenies). Considering that APOK3 is a protein targeted to mitochondria, it is tempting to hypothesize that variations in the mitochondrial genome can modulate pollen sensitivity to the killer. In this case, the mitochondrial genome could be considered as a 'modifier locus' of the PK. Interestingly, the antidote of the rice qHMS7 locus is also addressed to mitochondria [START_REF] Yu | A selfish genetic element confers non-Mendelian inheritance in rice[END_REF]) but the effect of the cytoplasmic background was not explored. To date, however, there is no indication that these poison-antidote systems from rice and Arabidopsis have functional similarities.

Our diversity survey revealed that some natural variants are neutral for the PK activity. Neutral alleles have been reported not only for gamete killers in rice [START_REF] Koide | Lineage-specific gene acquisition or loss is involved in interspecific hybrid sterility in rice[END_REF][START_REF] Yu | Hybrid Sterility in Rice (Oryza sativa L.) Involves the Tetratricopeptide Repeat Domain Containing Protein[END_REF][START_REF] Wang | Identification and confirmation of three neutral alleles conferring wide compatibility in inter-subspecific hybrids of rice (Oryza sativa L.) using near-isogenic lines[END_REF][START_REF] Liu | Identification of neutral genes at pollen sterility loci Sd and Se of cultivated rice (Oryza sativa) with wild rice (O. rufipogon) origin[END_REF]) and tomato [START_REF] Rick | The tomato ge locus: linkage relations and geographic distribution of alleles[END_REF]), but also for spore killers in fungi such as Neurospora [START_REF] Turner | Geographic distribution of neurospora spore killer strains and strains resistant to killing[END_REF]) and for the meiotic driver wtf in fission yeast [START_REF] Nuckolls | wtf genes are prolific dual poisonantidote meiotic drivers[END_REF]. We took advantage of the neutral behavior of Col-0 to exploit available mutants affected in each of the protein coding genes present in the interval that contains the sensitivity factor of the PK. This allowed us to establish that the PK3 functions as a poison-antidote system and to identify the gene encoding the antidote, APOK3, as AT3G62530. Inactivation of APOK3 in the neutral Col-0 allele turned it sensitive to a killer allele, as exemplified in hybrids with Mr-0 (Figure 8) and Ct-1 (Table 4), whereas hybrids with Col-0 or Sha had fully viable pollen and did not show any segregation distortion at the locus in their progenies (Table 4). These results demonstrate that the Col-0 allele of APOK3 protects pollen from the effect of killer alleles in hybrids, which is the definition of an antidote. In addition, three residues in the protein sequence and one SNP in the 5'-UTR are strictly associated with the antidote forms of the gene (Table 5), suggesting their involvement in its protective activity.

APOK3 molecular function still remains to be elucidated. It is annotated belonging to the ARM-repeat superfamily, which mediate numerous cellular processes including signal transduction, cytoskeletal regulation, nuclear import, transcriptional regulation and ubiquitination [START_REF] Samuel | Multifunctional arm repeat domains in plants[END_REF]. The remarkable features of this protein are its interacting domains, its mitochondrial location and its chimeric structure. APOK3 has three HEAT repeat domains (Figure 9A), predicted to form super-helixes [START_REF] Andrade | HEAT repeats in the Huntington's disease protein[END_REF] and to mediate protein-protein interactions [START_REF] Andrade | Comparison of ARM and HEAT protein repeats[END_REF]). More specific investigations will be needed to determine if these domains are involved in interactions with other mitochondrial proteins. In this context, it is interesting to note that APOK3 was reported to associate with ISOCITRATE DEHYDROGENASE 1, a regulatory subunit of NAD+-dependent enzyme of the tricarboxylic acid cycle (TCA), and SUCCINATE DEHYDROGENASE 4, a subunit of the mitochondrial respiratory complex II, in a study exploring the protein-protein interaction network of the plant TCA cycle [START_REF] Zhang | The Extra-Pathway Interactome of the TCA Cycle: Expected and Unexpected Metabolic Interactions[END_REF]). To date, the relevance to APOK3 function of its ability to bind zinc bivalent ions, reported by Tan et al. [START_REF] Tan | Divalent metal ions in plant mitochondria and their role in interactions with proteins and oxidative stress-induced damage to respiratory function[END_REF]) is difficult to assess. However, these indications make excellent entry points to functional studies aiming at discovering how APOK3 fulfills its antidote activity. Indeed, this question has been elucidated in only very few systems, among which the wtf spore killer of fission yeast in which the antidote was shown to coassemble with the corresponding poison and address the toxic aggregates to sequestering vacuoles [START_REF] Nuckolls | The wtf4 meiotic driver utilizes controlled protein aggregation to generate selective cell death[END_REF]. Interestingly, these authors also found that genes involved in the mitochondrial functioning counteract the poison activity, suggesting a role, yet to be uncovered, of mitochondria in this antidote mechanism. From our results, it seems that APOK3 has no other biological role than being the antidote for PK3. In Col-0 the apok3 KO mutation did not reveal any obvious phenotype alterations, at least in our laboratory standard conditions. Moreover, APOK3 is missing in the Moroccan accession Ita-0. This supports that APOK3 is not essential to A. thaliana. In addition to be absent in Ita-0, which belongs to an ancient lineage of A. thaliana (Arabidopsis Genome Consortium 2016; [START_REF] Durvasula | African genomes illuminate the early history and transition to selfing in Arabidopsis thaliana[END_REF]), APOK3 has no ortholog in the A. thaliana closest relative A. lyrata (Figure 7) and was not found in the available Brassicaceae sequences. This argues for this gene being specific to the A. thaliana species. It is thus likely that APOK3 originated after A. thaliana divergence from its common ancestor with A. lyrata. The gene structure suggests that APOK3 originated from a duplication of a member of the AT3G43260 and AT3G58180 gene family, followed by the recruitment of the AT3G62460 mitochondrial targeting sequence. In Oryza sativa ssp japonica, a comparative genome analysis detected 28 new O. sativa specific genes on chromosome 3, and 14 of which were found to be chimerical [START_REF] Zhang | High occurrence of functional new chimeric genes in survey of rice chromosome 3 short arm genome sequences[END_REF]. However, as far as we are aware, their biological functions remain unknown. Because of the absence of APOK3 in any related species, we could not conclusively infer whether non-antidote forms preceded antidote ones or rather derived from them. However, because A. thaliana Madeiran strains such as Are-1 and Are-10 have been described as archaic [START_REF] Fulgione | Madeiran Arabidopsis thaliana Reveals Ancient Long-Range Colonization and Clarifies Demography in Eurasia[END_REF], our phylogenetic analysis of APOK3 (Figure 10B) suggests that antidote forms evolved more recently. This would be an interesting question to address in order to understand the evolutionary trajectory of the gene in relation to its antidote function. Indeed, poison-antidote systems are often considered as selfish genetic elements, but their selfish nature is not always easy to properly establish, requiring thorough population and evolutionary genetic studies not amenable in all cases [START_REF] Sweigart | Making a Murderer: The Evolutionary Framing of Hybrid Gamete-Killers[END_REF]). In the present case, the question is still open, and deserves to be treated. It would involve the exploration of APOK3 sequence diversity in a wider sampling of natural variants, in particular in polymorphic A. thaliana natural populations. The examination of the organization and gene content of 11 genotypes of different PK3 status in addition to Col-0, Sha and Mr-0 provide interesting information, though, revealing the strikingly complex and variable structure of the PK3 locus in particular in killer alleles. In eukaryotes, poison-antidote elements are generally found within regions of high divergence and structural variation [START_REF] Burga | Toxin-Antidote Elements Across the Tree of Life[END_REF]). These are regions with low recombination (Larracuente and Presgraves 2012), which could favor the poison-antidote system. Indeed, the antidote must be inherited together with the killer. In our case, the structural differences between Sha and Mr-0 alleles easily explain the scarcity of recombinants found inside the PK3 interval during the fine-mapping of the locus. The use of these rare recombinants and of derived crosses allowed us to establish that the PK3 contains at least three genetic elements necessary for the PK activity (Figure 5), the antidote being flanked by two intervals each carrying elements necessary for the killer activity. Poison-antidote systems whose genetic factors have been identified so far in eukaryotes other than fungi most commonly involve two components [START_REF] Burga | Toxin-Antidote Elements Across the Tree of Life[END_REF]. However, in plants, three components poison-antidote systems have been previously described in rice distorter loci. In rice S1 [START_REF] Xie | An asymmetric allelic interaction drives allele transmission bias in interspecific rice hybrids[END_REF], three linked genes are necessary for toxicity; in rice S5 [START_REF] Yang | A killer-protector system regulates both hybrid sterility and segregation distortion in rice[END_REF], two sporophytic factors are necessary for the killer activity but they are not carried by the same allele. From our mapping results, we cannot completely exclude that PK3A or PK3C also carry rescue activities not redundant with PK3B, nor that PK3B also carries killer activity not redundant with PK3A and PK3C. However, it is noticeable that this structure prevents the killer activity to be isolated from the antidote by a recombination event, since such an event would also separate the two mandatory killer elements.

The PK3 locus is particularly dynamic and prone to structural variations. Besides the Ita-0 allele that is structurally similar to that of A. lyrata, with orthologous genes in the same order and orientation (Figure 7), the other killed and the neutral alleles (excepted Bur-0) present the same organization of the locus (Figure 6B), even if some protein coding genes are missing in Are-10 and Cvi-0 (Figure 7). Among the neutral alleles, Bur-0 is an interesting exception: it resembles a killer both in structure (Figure 6B) and gene content (Figure 7), and we hypothesize it evolved from a killer allele that lost its killing capacity. Comparison of the Bur-0 sequence with the very similar allele of Jea will probably help identifying killer elements of PK3. The most complex alleles are found in killers (Figure 6B), which are highly variable, with different groups of genes duplicated in one or the other orientation and in different relative positions (Figure 7). Having more than one copy of APOK3 is one remarkable common feature of the killer alleles. It is likely that having multiple copies is important for killer alleles not to be the victims of their own activity. Gene dosage has also been reported to be critical in a poison-antidote system in C. elegans [START_REF] Mani | A mechanistic basis for the coordinated regulation of pharyngeal morphogenesis in Caenorhabditis elegans by LIN-35/Rb and UBC-18-ARI-1[END_REF]. In addition, the different copies of APOK3 found in each killer are identical, but may be slightly different between accessions, indicating that recent duplications occurred independently in the lineages of the different variants we examined.

These results suggest that these alleles have experienced a more intense structural evolution that neutral and killed ones, and raises the questions of the mechanisms and evolutive forces leading to these structures. The presence of several transposable elements at the locus, and their variation in occurrence between Sha, Mr-0 and Col-0 (Figure 6A) could be relevant to this question, since TE mobility is a main source of genetic variation in A. thaliana [START_REF] Baduel | Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana[END_REF]). It has also been suggested that the proximity of transposons has facilitated duplication of the fission yeast wtf genes [START_REF] Eickbush | Killer Meiotic Drive and Dynamic Evolution of the wtf Gene Family[END_REF].

Even if gamete killers likely exist in all plant species, none had been investigated in A. thaliana until now, to our knowlege. The PK we dissected here has some general features of eukaryotic poison-antidote systems, including its species-specific nature and its presence within a hyper variable locus. The layout of the killer alleles is particular, with at least three mandatory elements and diverse duplications of sequence blocks that contain antidote genes trapped between killer elements. Continuing to exploit the natural variation of the species should help identifying the killer elements, and provide clues towards the underlying mechanisms responsible for the PK activity, the role of the mitochondria, and eventually the forces driving the evolution of the locus.

Materials and Methods

Plant material and growth conditions

A. thaliana natural accessions were provided by the Versailles Arabidopsis Stock Center (http://publiclines.versailles.inrae.fr/) and T-DNA lines were provided by the NASC (http://arabidopsis.info/).

In the analysis of the Sha/Mr bias at L3, all the genotypes used had the Mr-0 cytoplasm, unless specified. The ShaL3 H genotype was previously designated [Mr]ShaL3 H [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF]. It has the Mr-0 cytoplasmic genomes and a Sha nuclear background except at L3, i.e. between markers M1 and M15 (Figure 5 Source Data 3), where it is heterozygous. Here we constructed the MrL3 H genotype by backcrossing the (Mr-0 x Sha) F1 by Mr-0 three times. We selected a third backcross progeny heterozygous at L3, and at the loci L1 and L5 previously described in [START_REF] Simon | Genomic Conflicts that Cause Pollen Mortality and Raise Reproductive Barriers in Arabidopsis thaliana[END_REF] and we chose in its self-descent a plant homozygous Mr-0 at L1 and L5 while heterozygous from the marker M0 (Figure 5 Source Data 3) to the bottom of chromosome 3, which we called MrL3 H .

In order to save time, we constructed an early-flowering version of Mr-0 by inactivating FRIGIDA (AT4G00650) by CrispR-Cas9 [START_REF] Doudna | Genome editing. The new frontier of genome engineering with CRISPR-Cas9[END_REF]. We cloned two guide-RNAs (275rev and 981forw, Supplemental File 1) in a pDe-Cas9-DsRed binary vector [START_REF] Morineau | Dual Fatty Acid Elongase Complex Interactions in Arabidopsis[END_REF] and introduced this construct by floral dipping [START_REF] Clough | Floral dip: a simplified method for Agrobacterium-mediated transformation of Arabidopsis thaliana[END_REF] in Mr-0 plants. We isolated several T1 transformants, some of them flowering earlier than Mr-0. After segregating out the T-DNA, we selected one T2 plant homozygous for a stop mutation in the FRIGIDA first exon. This early version of Mr-0, called Mrfri, flowers 39 days after sowing in our greenhouse conditions. Mrfri produces viable pollen, whereas its cross with Sha induces pollen lethality, as does Mr-0 (Figure 8 Supplement 1). As expected, it results in a segregation distortion against Sha alleles at L3 in the progeny of Mrfri x Sha plants.

Plants were grown in the greenhouse under long-day conditions (16h day, 8h night) with additional artificial light (105 µE/m2/sec) when necessary.

Cytological analyses

DAPI staining of spread male meiotic chromosomes was performed according to [START_REF] Ross | A light microscopic atlas of meiosis in Arabidopsis thaliana[END_REF]. All stages were observed in 3 independent ShaL3 H plants and in parental controls.

Observations were made using a Zeiss Axio Imager2 microscope and photographs were taken using an AxioCam MRm (Zeiss) camera. Propidium iodide and Alexander staining [START_REF] Alexander | Differential staining of aborted and nonaborted pollen[END_REF] of pollen were performed as described in [START_REF] Durand | A restorer-offertility-like pentatricopeptide repeat protein promotes cytoplasmic male sterility in Arabidopsis thaliana[END_REF].

Fine-mapping and genotyping

DNA extractions were conducted on leaves from seedlings as described by [START_REF] Loudet | Bay-0 x Shahdara recombinant inbred line population: a powerful tool for the genetic dissection of complex traits in Arabidopsis[END_REF].

Markers for the fine-mapping are described in Figure 5 Source Data 1. For CAPS markers M3 and M12 we used Cac8I (NEB) and Bsp1407I (ThermoFisher) restriction enzymes, respectively. Other SNPs were genotyped by sequencing. For the fine-mapping of PK3, we genotyped a total of 4,717 plants and identified 42 recombinants between markers M1 and M15.

We selected 23 informative recombinants that were tested for segregation distortion at the locus by genotyping their self-descent progenies with appropriate markers (Figure 5 Source Data 1).

All the other primers used in this work for mutant characterization, gene expression, PCR amplification and DNA sequencing are listed in Supplemental File 1.

Purification of microspores

This step was carried out at the Imagerie-Gif facility (https://www.i2bc.paris-saclay.fr).

Microspores were isolated from flower buds by chopping with a razor blade in 0.1M Mannitol.

The crude suspension was filtered through a 50-µm nylon mesh (Sysmex-Partec) and collected in polypropylene tubes at 4°C. Microspores were sorted by flow cytometry using a MoFlo Astrios_EQ© cytometer (Beckman Coulter, Roissy, France) in PuraFlow sheath fluid (Beckman Coulter) at 25 psi (pounds per square inch), with a 100-micron nozzle. We performed sorting with ~43 kHz drop drive frequency, plates voltage of 4000-4500 V and an amplitude of 30-50 V. Sorting was performed in purity mode. The 488-SSC (Side Scatter) parameter was set as threshold. Microspores gate was first set on FITC-A (526/50 band-pass for autofluorescence) versus SSC plot. Then, only low Forward Scatter (FSC) events were kept. The singlet gate was established using autofluorescence parameters. Accuracy of gating was determined post-sorting using microscopy with transmitted light and 40X dry. The flow cytometer-sorted microspores (300k) were collected in 1.5-ml tubes containing 300 µL Trireagent® and conserved at -80°C until subsequent RNA extraction.

RNA extractions and RT-PCR

We mixed on ice 500 µl of microspore suspension in TriReagent with 300 µl glass beads (Sigma, G8772-100G glass beads acid washed) and two sterilized metals beads (3mm diameter). The tubes were shaked 2 x 2,5 min at a 1/30 frequency in a Retsch MM400 mixer mill and the solution without beads was recovered in a new RNase-free tube by making a small hole at the bottom of the tube and centrifuging 2 min at 3,900 rpm at 4°C. The beads were rinsed with 500 µl of TRIreagent and centrifuged again. The eluate was centrifuged 2min at 11,000 rpm, 4°C to remove cellular debris. RNA was extracted using the extraction RNA kit (Zymo research). Leaf RNA was prepared as above except that leaves were grinded in Trizol with one metal bead.

Reverse transcription was proceeded with 200 ng of purified RNA using the Maxima reverse transcriptase (Thermo Scientific) primed with oligodT in 50 µL. After heat inactivation of the RT, 2 µL of cDNA were used for each PCR. (http://www.codoncode.com/aligner/). For Mr-0, because of the complexity of the locus, a fosmid library was built using the "Copycontrol fosmid library production kit" with pCC2FOS TM vector and T1 EPI300 TM E. coli strain (Illumina technologies). The library titled 7.600 cfu/mL with an average size of inserts of 35 Kb. Purification of fosmid DNA was performed using the FosmidMax DNA purification kitTM (Illumina technologies).

Sanger sequencing and annotation of

The transfer of structural and functional annotation from the Col-0 reference genome to Sha and Mr-0 was done with the EGN-EP transfer pipeline [START_REF] Sallet | EuGene: An Automated Integrative Gene Finder for Eukaryotes and Prokaryotes[END_REF]). These annotations were manually corrected, especially in Mr-0, to take into account copy number variations. TE annotation was performed with the TEannot tool [START_REF] Quesneville | Combined evidence annotation of transposable elements in genome sequences[END_REF], http://urgi.versailles.inra.fr/) and hits of more than 1 Kb were retained.

De novo genome sequencing, assembly and annotation

High molecular weight DNA was extracted from 3 week-old plantlets using a protocol modified from [START_REF] Mayjonade | Extraction of high-molecular-weight genomic DNA for long-read sequencing of single molecules[END_REF] which is described on Protocols.io [START_REF] Russo | HMW DNA extraction from diverse plants species for PacBio and Nanopore sequencing[END_REF]). The following library preparation and sequencing were performed at the GeT-PlaGe core facility, INRAE Toulouse.

Nanopore sequencing: ONT libraries were prepared using the EXP-NBD103 and SQK-LSK109 kits according to the manufacturer's instructions and using 4 µg of 40Kb sheared DNA (Megaruptor, Diagenode) as input. Pools of six samples were sequenced on one R9.4.1 flowcell.

Between 0.014 and 0.020 pmol of library were loaded on each flowcell and sequenced on a PromethION instrument for 72 hours.

Illumina sequencing: Illumina libraries were prepared using the Illumina TruSeq Nano DNA HT Library Prep Kit according to the manufacter's instructions. Libraries were then sequenced with 2x150bp paired-end reads on an Hiseq3000 instrument (Illumina).

Nanopore sequence datasets ranging from 4Gb to 14Gb with a minimum N50 of 22Kb were assembled with CANU [START_REF] Koren | Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation[END_REF]) software release 1.9 (genomeSize=125M parameter). Mitochondrion and chloroplast genomes were assembled with CANU after selection of previously corrected CANU reads longer than 40Kb mapped with minimap2 (Li 2018) (-x asm5) on a databank built with Col-0 chloroplast and Col-0 and Landsberg mitochondrion genomes (NCBI accessions: NC_000932, NC_037304, JF729202). Spurious contigs were identified with minimap2 (-x asm5) and removed from the raw CANU assemblies.

A spurious contig was defined as either a contig mapped on mitochondrion or chloroplast genomes or a contig mapped on a larger contig with a hit spanning at least 80% of its length.

Nuclear contigs were scaffolded with AllMaps [START_REF] Tang | ALLMAPS: robust scaffold ordering based on multiple maps[END_REF] using Col-0 genome as reference. Then, two rounds of consensus polishing using Illumina paired-end data were performed. For each round, the paired-end reads were first mapped with bwa (Li and Durbin 2009) (0.7.17-r1188 debian) with a minimal score of 50 (-T) then pilon software [START_REF] Walker | Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement[END_REF]) (version 1.23) was used to generate the polished consensus sequences (--fix snps,indels --flank 20 -min-depth 20). Unanchored contigs shorter than 40kb were removed (considered as likely contaminants as very few Illumina reads were mapped on these short contigs).

The genomes were annotated using the EuGène software version 4.2a via the integrative pipeline egnep [START_REF] Sallet | EuGene: An Automated Integrative Gene Finder for Eukaryotes and Prokaryotes[END_REF]) version 1.6 (http://eugene.toulouse.inra.fr/Downloads/egnep-Linux-x86_64.1.6.tar.gz). In addition to standard tools for repeat masking and lncRNA prediction included in egnep, the most recent and comprehensive annotation of A. thaliana available was used both as training dataset and as source of evidences. The peptide database of Araport (Krishnakumar et al. 2015) version 11 (201606) was used for the similarity searches performed with NCBI-BLASTX software version 2.2.31+ (hits spanning more than 80% of the protein length were retained). The corresponding cDNA database (201606) was used as transcriptional evidences. The cDNAs were mapped with gmap [START_REF] Wu | GMAP: a genomic mapping and alignment program for mRNA and EST sequences[END_REF] version 2017-09-05 (hits spanning more than 30% of the transcript length at a minimum identity percentage of 94 were retained).

In accessions with several copies of APOK3, occasional sequence errors subsist in some copies in the automatic assemblies. To correct these errors, APOK3 was Sanger sequenced after gene amplification with the primers AT3G62530F1 and AT3G62530R1 (Supplemental File 1).

Sequence availability

Whole genome sequences and PK3 locus sequences are available at the doi listed in Figure 6 Source Data 1. Start and end of the PK3 locus were determined by homology with the Col-0 reference sequence. The locus is represented as a horizontal line with the killer and killed alleles as black and white boxes, respectively. For these general models, there is no hypothesis on the number of genes present in each box. In both models, the black allele is expressed before meiosis, and the killer (K) or poison (P) persists in all the gametes. In the killer-target model, only gametes with the white allele express the target (T), which interacts with K to trigger gamete death (represented as dashed outline). In the poison-antidote model, the black allele also produces a short-lived antidote (A) that counteracts P, so only the gametes that inherit the black allele are efficiently protected.

  Sha and Mr-0 PK3 locus PCR and sequencing primers are listed in Supplemental File 1. Amplification products were sequenced by Beckman Coulter Genomics (http://www.beckmangenomics.com)

Figure Supplements : Figure 8

 :8 Figure Supplements: Figure 8 Supplement 1: Mrfri has the same killer behaviour as Mr-0.

Figure 8

 8 Figure 8 Supplement 2: Genetic and phenotypic characterization of the Col mut530 mutant.

Figure 4

 4 Figure 4 Source Data 1: countings of abortive and normal grains in developing and mature pollen of plants with PK3 and parental controls.

Figure 5

 5 Figure 5 Source Data 1: Segregation analyses in the self descent of 23 recombinants identified in PK3 fine mapping.

Figure 5

 5 Figure 5 Source Data 2: Segregation at L3 in self progenies of hybrid plants derived from recombinants 52D7, 8F10BH2 and 23G9.

Figure 5

 5 Figure 5 Source Data 3: Genetic markers.

Figure 6

 6 Figure 6 Source Data 1: PK3 locus DNA sequence resources.

Figure 8

 8 Figure 8 Source Data 1: annotated genes in the PK3B interval.

Figure 8

 8 Figure 8 Source Data 2: Segregation analyses of Mrfri x Col mut F2 families.

Figure 8

 8 Figure 8 Source Data 3: Transmission of the Col mut499 and Col mut530 alleles from female and male sides.

Figure 8

 8 Figure 8 Source Data 4: raw gels of RT-PCR experiments on PK3B genes.

Figure 9

 9 Figure 9 Source Data 1: Segregation at L3 in reciprocal self progenies of hybrid plants derived from recombinants 52D7, 8F10BH2 and 18D1BC7.

FigureFigure 1 :

 1 Figure 10 Source Data 1: multifasta of APOK3 gene sequences from 27 accessions of knwon status at PK3.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Plants heterozygous at L3 in either Sha or Mr-0 nuclear backgrounds show dead pollen. Typical pollen phenotypes (Alexander staining of anthers) of ShaL3 H and MrL3 H plants compared to the parental accessions Sha and Mr-0. Viable pollen grains are stained in red, aborted pollen grains appear in blue. Scale bars: 100µm. Graphical representations of ShaL3 H and MrL3 H genotypes are displayed on the right.

Figure 5 :Figure 8 :Figure 8 Figure 8 Figure 9 :Figure 9 Figure 10 :

 58889910 Figure 5: Fine mapping of the PK3 locus. A: The upper panel shows the mapping of PK3 in the L3 locus delimited by the genetic markers M1 and M15 (Figure 5-Source Data 1). Markers M5 and M13 delimit the PK3 interval, where all the PK elements required to cause the segregation bias are present. Positions of markers on the TAIR10 Col-0 genomic sequence are given. The most relevant recombinants are shown, with the frequencies of Sha homozygotes (f Sha ) observed in their selfed progenies indicated on the right. TRD, transmission ratio distortion (chi square test for 1:2:1 segregation) *** p < 0.001, NS: not significant. The lower panel shows the further dissection of the PK3 locus into PK3A, PK3B and PK3C (Figure 5-Source Data 2). Because some markers are very close to each other (Figure 5-Source Data 1), this diagram is not to scale. Genotypes used to dissect the PK3 are presented: i-52D7 and i-8F10BH2 result from crosses of 52D7 and 8F10BH2, respectively, with Sha; the 23G9#15 x 52D7#7 genotype, with two heterozygous regions flanking a fixed central region of the PK3, was produced by crossing appropriate selfed progenies of 23G9 and 52D7 recombinants. All markers are described in Figure 5-Source Data 3. B: Pollen viability in recombinants used to delineate the PK3 interval. Representative images of Alexander staining of anthers from recombinants presented in A. Aborted pollen grains are blue whereas living pollen grains are red. Scale bar: 50µm.

Table 1 :

 1 . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in Killed or killer behaviors of natural accessions

	bioRxiv preprint	doi:	; https://doi.org/10.1101/2021.09.08.459474	this version posted September 10, 2021.	The copyright holder for this

Table 2 : Segregation analyses at L3 of selfed progenies in two different nuclear backgrounds

 2

	Genotype	Sha	Number of plants (1) Hz Mr-0	Total	p χ 2 (1:2:1)	f Sha (2)
	MrL3 H	14	83	73	170	1.2 10 -9 ***	0.08
	ShaL3 H (3)	18	87	73	178	4.0 10 -8 ***	0.10

Table 3 : Transmission of the Sha allele in a heterozygous context from the male and female sides

 3

	Mother plant	Number of plants (1) Sha Hz Total	p χ 2 (1:1)	f Sha (2)
	Sha x ShaL3 H	40	132	172	2.3 10 -12***	0.23
	ShaL3 H x Sha	78	90	168	0.4 NS	0.46

  International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

	bioRxiv preprint	doi:	; https://doi.org/10.1101/2021.09.08.459474	. CC-BY 4.0 The copyright holder for this this version posted September 10, 2021.
	pollen division (Figure

Table 4 : Segregation at L3 in F2s where the Col mut530 allele is confronted to a killed (Sha), neutral (Col-0) or killer (Ct-1) allele

 4

			Number of plants
	Genotype	not				p χ 2 (1:2:1)	f Col mut530 (4)
		Col mut530 Hz Col mut530 Total
	F2 (Sha x Col mut530 ) (1)	46	101	37	184	0.3 NS	0.20
	heterozygote Col mut530 /Col-0 (2)	49	85	46	180	0.5 NS	0.26
	F2 (Ct-1 x Col mut530 ) (3)	74	89	21	184	2 10 -7 ***	0.11

Table 5 : APOK3 haplotype diversity among 27 accessions of known status for PK3 SNP position (2)

 5

		Accessions (1)	-36 6	17 27 56 68 152 241 248 253 332 341 363 373 378 384 415 583 765 840
		Col-0, Shigu-2, Oy-0, Blh-1, N16, Koch-1	c	T	T	T	T	C	T	G	A	A	a	t	A	A	T	C	G	G	t	T
	antidotes	Mr-0, Etna-2 Ct-1 Bur-0 APOK3-1 & 2	. . .	. . .	. . .	. . .	. . .	. . .	. . .	T . .	. . .	. . .	. . .	. . .	. T .	. . .	. . .	. . .	. . .	. . C	. g g	. . .
		Jea, Cant-1, Bur-0 APOK3-3, Lov-5	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	.	g	.
		Cvi-0, N13, Kly-1, Nov-01, Zal-3, Stepn-1	t	.	.	.	.	.	.	.	.	T	.	.	.	T	.	T	.	.	g	.
	non-antidotes	Kidr-1 Sorbo, Kz-1, Rak-2 Sha Kas-2	t t t t	. . . .	. . . .	. . . .	. . . C	. . A .	. A A .	. . . .	. . . C	T T T T	. . . .	c . . .	. . . .	T T T T	. . . G	T T T T	. T T .	. . . .	g g g g	. . . G
		Are-10, Are-1	g	C	G	C	.	.	.	.	C	T	g	.	.	T	G	T	.	.	g	.
		Amino acid change position			6		19 23 51 81 83 85			98 101 103 105 115 141
		antidote

Table 4

 4 Supplement 1: Control of segregation at L3 in the progenies of F1 (Ct-1 x Col-0).

	Figure 9 Supplement 1: Tree of protein sequences for deoxyhypusine hydroxylases from
	different taxa and their homologues in A. thaliana.
	Source Data Files:
	Table 1 Source Data 1: Segregations at L3 in progenies of crosses of Mr-0 with different natural
	accessions.
	Table 2 Source Data 2: Segregations at L3 in progenies of crosses of different natural accessions
	with Sha.

Acknowledgments

We gratefully acknowledge Adeline Simon for her help in bioinformatics, Fabien Nogué for his advices in the CrispR Cas9 mutagenesis and Filipe Borges, Nicolas Valentin and Mickaël Bourge for the microspore FACS sorting. Sebastien Santini (CNRS/AMU IGS UMR7256) and the PACA Bioinfo platform (supported by IBISA) are acknowledged for the availability and management of the phylogeny.fr website. The present work has benefited from Imagerie-Gif core facility supported by l'Agence Nationale de la Recherche (ANR-11-EQPX-0029/Morphoscope, ANR-10-INBS-04/FranceBioImaging; ANR-11-IDEX-0003-02/ Saclay Plant Sciences) and from the support of IJPB Plant Observatory technological platforms. The IJPB benefits from the support of Saclay Plant Sciences-SPS (ANR-17-EUR-0007). This work was funded in part by INRAE Biology and Plant Breeding department (EVOLOX and POLLEN grants) and by the Région Midi-Pyrénées (CLIMARES project) and the LABEX TULIP (ANR-10-LABX-41).

Competing interest statement

We declare no competing interests. A: PK3 locus in Col-0, Sha and Mr-0. Beige and blue plain arrows represent protein coding genes and TEs, respectively, with their orientations. Gene labels on the Col-0 locus correspond to the last three digits of AT3G62xxx gene identifiers in the reference sequence (TAIR10), the same numbers on Sha and Mr-0 loci indicate genes homologous to the Col-0 genes. Black, red and blue dotted lines respectively represent insertions/deletions, inversions and duplications. The limits of the PK3A, PK3B and PK3C intervals in Sha and Mr-0 are indicated (dashed broken lines). B: Dot plots of PK3 sequences of 12 accessions (x axes) against the reference Col-0 (y axes), generated at NCBI using the nucleotide Blast tool with the option 'align two or more sequences' and the following settings: max target = 10, Expect threshold = 0.001, word size = 256 (Altschul et al. 1990). The names of killer, neutral and killed accessions are in blue, green and red, respectively. The sequences can be accessed by DOI numbers indicated in Figure 6 Source Data 1, with positions of the PK3 sequences. Alignment the PK3 loci from thirteen A. thaliana accessions and two related species, Arabidopsis lyrata and Boechera stricta, drawn to highlight synteny between protein coding genes. The scheme is not to scale and TEs are not represented. For A. lyrata and B. stricta the structural annotation was obtained from Phytozome (https://phytozome.jgi.doe.gov/). The names of killer, neutral and killed A. thaliana accessions are written in blue, green and red, respectively. Plain coloured arrows represent coding genes with their orientations, each colour representing orthologues and paralogues of a same gene. Gene labels correspond to the last three digits of AT3G62xxx, AL5G45xxx and Bostr.13158s0xxx gene identifiers in A. thaliana (Col-0 reference sequence TAIR10), A. lyrata (V2.1) and B. stricta (V1.2), respectively. Black arrows above the Mr-0 and Sha loci delimit the PK3A, PK3B and PK3C intervals from left to right.