Cosegregation of recombinant chromatids maintains genome-wide heterozygosity in an asexual nematode

In asexual animals, female meiosis is modified to produce diploid oocytes. If meiosis still involves recombination, this is expected to lead to a rapid loss of heterozygosity, with adverse effects on fitness. Many asexuals, however, have a heterozygous genome, the underlying mechanisms being most often unknown. Cytological and population genomic analyses in the nematode Mesorhabditis belari revealed another case of recombining asexual being highly heterozygous genome-wide. We demonstrated that heterozygosity is maintained despite recombination because the recombinant chromatids of each chromosome pair cosegregate during the unique meiotic division. A theoretical model confirmed that this segregation bias is necessary to account for the observed pattern and likely to evolve under a wide range of conditions. Our study uncovers an unexpected type of non-Mendelian genetic inheritance involving cosegregation of recombinant chromatids.

The PDF file includes:

Figs. S1 to S9 Table S1 Supplementary text Legends for data S1 to S3 References
Other Supplementary Material for this manuscript includes the following: Figure S1: Quantification of heterozygosity and estimation of r.Heterozygosity (A) and recombination (B) estimated along a representative contig (contig 110) and distribution of heterozygosity (C) and recombination rates (D) computed on 5000bp windows.Genome-wide average heterozygosity and recombination rate r are represented by a red bar.Other sexual species, taken from the literature, are depicted with a blue bar: a) Homo sapiens (38,39); b) Gasterosteus aculeatus (40); c) Mus musculus (41,42); d) Drosophila melanogaster (20); e) Taenopygia guttata (43); f) Heliconius melpomene (44,45).The assembly comprises 218 contigs, a single representative is shown here.The rho and heterozygosity values for windows of 5000 base pairs for each of the 218 contigs are available in the zenodo archive: https://doi.org/10.5281/zenodo.8112093

Figure S2: Design and expectations for the Edu experiment.
A) Design of EdU pulse-chase method.EdU is pink and DNA is blue.Sketch of M. belari gonad showing a pair of homologous chromosomes at each stage of oogenesis, with each chromatid and the two DNA strands per chromatids (Watson and Crick strands shown in light colors).Before EdU incorporation, each chromatid contains two unlabelled DNA strands (in light blue).Chromatids therefore appear blue.During the next S phase in the presence of EdU, the newly synthesized strains are pink.Hence, each chromatid appears pink.After chromatid segregation and a new S phase without EdU, one chromatid is entirely blue (because it has inherited the first unlabelled DNA strand and a new unlabelled DNA strand has been produced), whereas the other chromatid is pink (because although a new unlabelled DNA strand has been produced, it has inherited the labelled DNA strand from the previous cell cycle).The duration of the EdU exposure and progression through meiosis without EdU has been optimized to obtain such pattern.B) The analysis of chromosome colors at diakinesis strongly suggest that in fact, after the first S phase with EdU, the cells undergo two rounds of S phase in the absence of EdU.The aim of the model is to predict population genomic patterns expected under the Mesorhabditis life cycle.We consider a single neutral locus with an infinite allele model (IAM) of mutation.We consider a subdivided population with K demes, each of size N .The total population is thus N T = KN .The life cycle is as follows (Figure S4 redrawn from the main text).The generations are discrete and non-overlapping.Migration occurs before reproduction according to the island model at rate m.This simple population structure allows modelling the breeding structure of Mesorhabditis with high level of consanguineous mating: the lower the deme size and the migration rates, the higher the level of consanguineous mating.Automictic reproduction that leads to gynogenetic females occurs in proportion σ and sexual reproduction 1 − σ.Under automixis, the rate of loss of heterozygosity (LOH) is β and depends on the underlying mechanism.With standard central fusion β = 1/2 after one crossover and random assortment of chromatids, and β = 0 if there is no recombination.In the proposed model of co-segregation of recombinant chromatids (CRC), β = 1 2 (1 − b), where b represents the segregation bias.If b = 0, there is no bias so that β = 1/2 as for recombinant chromatids under standard central fusion.If b = 1, the segregation is fully biased towards the co-segregation of the two recombinant chromatids.In this model, negative b values could also be considered, with b = −1 corresponding to complete LOH after one meiosis.However, we did not explore this possibility that is not relevant for the Mesorhabditis system.Under sexual reproduction, the proportion of females and males is φ and 1 − φ, respectively.So far, no sexually-produced female has been observed (corresponding to φ=0).However, it is possible that they are produced at low rate in natural populations.If we set σ = 1 and β = 0, this is equivalent to a fully clonal model.If we set σ = 0 and φ = 1/2, this is equivalent to a fully sexual model.

Supplementary Text
To obtain measures of genetic diversity and population structure (F-statistics) we derive recursions on a set of probabilities of identity by descent (IBD), Q i , which leads to heterozygozity and F-statistics measure of the form [51]: and To fully describe the model, we need to follow eight probabilities of IBD, noted Q k 0 when the two genes are sampled in the same individual, Q k 1 when they are sampled in two individuals of the same population, and Q k 2 , when they are sampled in two individuals from different populations.The superscript k stands for the sex of sampled individuals (see Figure S5).
From the parameters describing reproductive pathways (Figure S4) we obtain the proportion of males, sexual and gynogenetic females as: from which we obtain the number of males and females in the populations and the proportions of gynogenetic females among females We also need to introduce the following compound parameters: Υ corresponds to the probability that the two sampled genes have not mutated in one generation, a, respectively b, corresponds to the probabilities that two individuals sampled in a same population, respectively in two different populations, were in the same population before migration.

Recursions
We now write the recursions of the Q i (t + 1) as a function of the Q i (t).At equilibrium Q i (t + 1) = Q i (t) so we remove the time subscript and directly give the equations at equilibrium: (continued on next page) (continued from previous page) As an example, the rationale of the derivation is given for Q f 0 (t + 1).The two copies sampled in a female are IBD first if none has mutated (Υ).Then we must consider that this female is gynogenetic (R g ) or sexual (1 − R g ).If it is gynogenetic, if heterozygozity has been lost (β) the two copies are IBD with probability one, otherwise (1 − β) the probability of IBD is the same as for the mother, so Q f 0 (t).If the female come from sexual reproduction, the two copies are IBD with the same probability as of a random male/female pair at the previous generation, Q mf  1 .The terms in 1/N m and 1/N f that appear in equations for the Q k 1 and Q k 2 correspond to the probability that two different individuals have the same father or mother, respectively.
This system of recursions can be written in the matrix form: where Q is the vector of probabilities of IBD,G is a matrix and C a vector, both depending of the parameters of the model.The solution can be written on the form: where I is the identity matrix.
The F IS statistics measured on genomic data corresponds to F 0,2 in our model as we compare the IBD of two gene copies sampled either within an individual or at random over the whole population.It can be defined either for males, females, or for the whole population by weighting as a function of the sex-ratio:

Simulations
In addition to analytical derivations, an individual-based, multi-locus model was implemented using the SLiM software [50], where each individual's genome is explicitly defined.It allowed to simulate genomic patterns along a chromosome and to introduce deleterious mutations in the model.Each individual is represented by a unique pair of autosome.This chromosome consist of L = 10 5 loci, with a rate of recombination per locus r equal to 1/L = 10 −5 , so that on average one recombination event is observed per gametogenesis.Since M. belari has holocentric chromosomes (there is no specific centromere), the location of the chiasma is randomly drawn according to a uniform distribution along the chromosome.Two types of mutations are introduced, either neutral (no effect on fitness) or deleterious ones, with selection coefficient s, set to 0.01, and dominance coefficient h, set to 0.25, and acting multiplicatively across the genome.The life cycle was the same as described above.Various population sizes and migration rates were explored to modulate the level of inbreeding and we specifically studied the effect of the rate of production of sexual females, φ and the rate of LOH, β.
From the simulations we computed F IS and the linkage-disequilibrium measure r 2 on neutral mutations at the scale of the whole population: with H o the observed heterozygosity (percentage of heterozygous sites) and H e the expected heterozygosity : a sum over all the n mutations present in the population with f i their respective frequencies.
where D = π AB − π A π B with π A and π a the allelic frequency at a first locus A and π B and π b at a second locus B. Only mutations in frequency higher than 5% in the population were used.r 2 was calculated for pairs of loci as the function of their distance on the simulated chromosome.

F IS
The general solution of equation ( 8) can be obtained with the help of Mathematica [52] but it is formidable, so useless for direct biological interpretation.We thus performed numerical explorations in the general case.We also obtained approximations under two limit conditions: φ = 0 (no sexually-produced females) and β = 0 (no LOH).Under these two conditions we also used the standard diffusion limit.First we used the following scaling parameters: Then, we assumed an infinite number of local populations so K, N T → ∞ but that the scaled parameters terms tends towards constant: θ, Φ, B → cte.We also noted M = 4N m, where migration is scaled with the local, N , which can be small.
Assuming no sexually-produced females we obtained: Assuming no LOH we obtained: We give the three expressions for completeness but we can concentrate on the expression for the average F IS In natural populations, there is a strong family structure with most matings occurring between kins.In our model, it corresponds to low migration between demes, so M close to 0. Then, equations (13c) and (14c) become: which both reduce to F IS = (1 − σ)/2 when B = 0 or Φ = 0.As σ ≈ 0.9 in natural populations it corresponds to F IS ≈ 0.05, so close to the observed value (F IS ≈ 0.19, see main text).In an unstructured population (M → ∞, and see N = 1000 on Figure S6), F IS tends to −1 with no LOH (B = 0) (see also [48]).Here, this is compensated by the strong family structure leading to F f IS close to 0 as illustrated on Figure S6 (see also [?]).It is worth noting that just a little bit of sex or LOH (higher than the mutation rate: φ, B > θ) rapidly leads to F f IS close to 1.This is confirmed by simulations as presented (Figure S6 and Figure 4 in the main text).This result is thus very sensitive to the occurrence of sexually produced females and to low level of LOH.For comparison, we also present the case of standard central fusion, clonality and fully sexual reproduction.However, simulations showed that because of deleterious mutations, the range of parameters leading to F f IS close to 0 can be much wider, even when some females are sexually produced and when the CRC is not complete (Figure S7).This is explained by the selection against highly homozygous individuals.This is an important result as a non-zero proportion of sexually produced females is required to explain the low level of linkage disequilibrium observed genome wide as presented below.For comparison, Figures S6 and S7 also shows the results for full sexuality, full clonality and standard central fusion automixis (without CRC).

Linkage disequilibrium
Analytical derivations for linkage disequilibrium would require recursion equations for 40 IBD probabilities.We thus only relied on simulations.
To understand the specific effect of the CRC we first considered a single unstructured population.As expected, LD is low and rapidly decreases with physical distance under sexual reproduction (Figure S8).In contrast, LD is high and independent of physical distance along chromosome under both full clonality because recombination does not occur, and under standard central fusion because recombination is not efficient as individuals are mostly homozygotes.Under complete CRC and without sexually produced females, LD pattern is intermediate: much lower than under clonality or standard central fusion but higher than under full sexuality, and decreasing with physical distance (Figure S8).The reason is that, although new haplotypes are generated by recombination within individuals, they are never associated together through mating, so recombination is not as efficient as under full sexuality.
When we also consider the effect of population structure, LD is globally higher and the difference among reproductive modes are less clear under pure neutrality.However, when deleterious mutations are added the differences become stronger.In particular a very low rate of sexually produced females is sufficient to make CRC similar to full sexuality whereas much higher rates are necessary to erase the signature of clonality or standard central fusion (Figure S9).

Conclusion
Overall, the observed genomic pattern is well predicted by the strong family structure of M. belari populations directed chromatif assortment, sexually produced females at low rate (which can be unobserved in natural conditions), and the occurrence of deleterious mutations

Evolution of the reproductive system
In the second model, we study the evolution of the atypical reproductive system of M. belari from a standard sexual population.Compared to the previous model, the parameters σ (proportion of asexually produced offspring) and β (rate of LOH during automixis) are no longer fixed but under the control of two independent evolving loci.Among sexually produced individuals the proportion of females is set to φ = 1/2.Additional simulations where the sex-ratio can also evolve do not change the results and a sex ratio biased towards males evolved in a second step (not shown).
We started with a burn-in period where the population evolved under sexual reproduction.Then mutations are introduced at one or at the two loci controlling the reproductive mode.Each evolving locus is unlinked with any other locus.Transmission from parents to children of an evolving locus follows the same rules as transmission of loci on the focal chromosome, taking into account the potential LOH through chromosome segregation during automixis.It can take two different values, a resident value, g r , or a mutant value, g m = g r ± δ.We used δ = 0.1 to avoid too long simulations.We assumed additivity between the two alleles so the genotypic value is the sum of the two alleles.At the beginning of simulations, all individuals are homozygous for the resident allele and a mutant allele is randomly introduced in one individual of the population.If this mutant becomes fixed in the population, this mutant allele becomes the resident allele.If the mutant disappears, another mutant is directly reintroduced into the population at random.

Genotype to phenotype map
The value of σ and β are constrained between 0 and 1.The simulations start with a sexual population, with a rate of pseudogamy set to 0 and that can theoretically go up to 1.The rate of LOH β starts at 1/2 and can reach 0 (no LOH, complete cosegregation of recombinant chromatids) or 1 (full LOH, complete segregation of one recombinant and one non-recombinant chromatids).To map allelic values from ] − ∞, +∞[ to [0, 1] we used the following functions: with The parameters a and b define whether the function tends more or less quickly to 1 as x increases.We used a = 0.1 and b = 15 but equilibrium results are not sensitive to chosen values.The rate of pseudogamy and the rate of LOH are controlled by the maternal genotype.

Inbreeding depression
The supposedly strong barrier to the evolution of asexuality (assuming that automixis mechanisms can emerge) is inbreeding depression (ID) due to LOH that leads to the exposure of recessive deleterious alleles.ID is calculated as 1 − w 1 /w 0 with w 0 being the fitness of sexually produced individuals and w 1 the fitness of asexually produced individuals.
A strong ID should hold back the evolution of automixis if the chromosome segregation bias is not strong, preventing LOH.Since we only modelled a single chromosome and central fusion leads to a LOH in only half of offsrpings, the maximum ID is 50%, which is the minimum value preventing the evolution of asexuality.Asexuality should always evolved under such a conditions.To allow a broader range of ID values we assumed an additional cost of LOH and corrected w 1 to w corr 1 = (1 − c) β(n−1) w 1 , where n is the number of chromosomes (fixed to 10) and c is the mean fitness decline caused by LOH per chromosome.

Extinction-recolonisation cycle of subpopulations
In this model, a higher migration rate is set, allowing the evolving loci to be transmitted not too slowly from one population to another.It should also be noted that here the parameters controlling reproduction are no longer fixed rates.It is thus possible that populations may stochastically run out of males, especially if pseudogamy evolves and the percentage of males decreases.If a population has no more males, it becomes extinct.An extinct population is recolonised in the generation following its extinction by a male and at least one female randomly selected from the other population.This pair reproduces and together they form the next generation of the subpopulation.

Figure S3 :
Figure S3: Evidence of Directed Chromatid Assortment.11 images of fixed embryos at one-cell stage.EdU is in pink and DNA is in blue.The count of bicolored chromosomes is summarized in Figure 3.The dotted line shows the chromosome longest axis.Scale bar is 5 µm.Some chromosomes are shown twice, on the transverse and lateral view to better show the chromatid axis.The polar body (with labelled EdU chromatids, as expected) is visible in embryos #1 and #5.The sperm DNA (unlabeled) is visible in embryos # 3, #6 and #7.

Figure S4 :
Figure S4: Mesorhabditis life cycle and parameters of the model.

Figure S5 :
Figure S5: Definition of the eight probabilities of IBD.The grey boxes corresponds to demes.

Figure S6 :
Figure S6: F IS for the different reproductive modes as a function of the rate of sexually produced females (2KN φ) for different levels of population structure (KN = 1000, so N = 1000 corresponds to a single populations and N = 20 to a highly structured one).

Figure S7 :
Figure S7: F IS for automixis with different level of CRC (β = 0 corresponds to full co-segregation and β = 1/2 to random pairing).

Figure S8 :
Figure S8: Pattern of linkage disequilibrium as a function of physical distance (in number of loci) for different rates of sexually produced females (2N Kφ = 0.001, 1, 100) and for different reproductive modes.

Figure S9 :
Figure S9: Pattern of linkage disequilibrium as a function of physical distance (in number of loci) for different rates of sexually produced females (2N Kφ = 0.001, 1, 100) and for different reproductive modes.