Allorecognition genes drive reproductive isolation in Podospora anserina

Allorecognition, the capacity to discriminate self from conspecific non-self, is a ubiquitous organismal feature typically governed by genes evolving under balancing selection. Here, we show that in the fungus Podospora anserina, allorecognition loci controlling vegetative incompatibility (het genes), define two reproductively isolated groups through pleiotropic effects on sexual compatibility. These two groups emerge from the antagonistic interactions of the unlinked loci het-r (encoding a NOD-like receptor) and het-v (encoding a methyltransferase and an MLKL/HeLo domain protein). Using a combination of genetic and ecological data, supported by simulations, we provide a concrete and molecularly defined example whereby the origin and coexistence of reproductively isolated groups in sympatry is driven by pleiotropic genes under balancing selection.


Individual based simulation of the het-r/het-v interaction in P. anserina
This SLiM v. 3.3.2 1 simulation aims at exploring possible mechanisms through which the interaction of the vegetative incompatibility genes het-v and het-r can result in the formation of two segregated mating groups within a single population of the fungus P. anserina.

Life cycle
The life cycle of Podospora is modelled using two populations. This approach is necessary to take into account critical features of the life cycle of Podospora. Population 1 represents the diploid stage of the life cycle and population 2 the dikaryotic stage. To complete the life cycle, individuals need to cycle through both populations. Diploid individuals of population 1 transition to population 2 through a selfing phase, which is effectively equivalent to the meiosis leading to the dikaryotic stage in the life cycle of Podospora. Individuals of population 2, representing the dikaryotic stage, then produce gametes but are not allowed recombination since in the life cycle of Podospora recombination has already occurred at the transition between the diploid and dikaryotic stages. The gametes are then used to form the new diploid individuals of population 1, completing the life cycle.

Genetic architecture
Each individual carries two haploid genomes, including the het-v and het-r loci on two separate chromosomes. Each locus can be represented by two alleles, V and V1 for het-v and R and r for het-r.

Fitness effects of genotypes
The simulation takes into account the main features of the het-r/v interaction: (i) Lethal genotype: We assumed that the RV combination within a genotype is lethal. This means that after meiosis dikaryons of the genotype RV -XX effectively become monokaryons of genotype XX. This is based on the observation that dikaryons can sometimes disassociate into "sectors" within a culture that have each a different nucleus.
(ii) Balancing selection: Due to vegetative incompatibility, individuals carrying the R allele of het-r have an advantage over individuals carrying the r allele, which is an increasing function of the frequency of V in the population (Supplementary Fig. 16a). The benefit in relative fitness of V over V1 due to the presence of R is identical and also follows the function of panel Supplementary Fig. 16a (with frequency of R on the x-axis and benefit to V on the y-axis). In addition, vegetative incompatibility between V and V1 provides a benefit in relative fitness to the rarer of the two alleles (Supplementary Fig. 16b). See SLiM code for fitness functions in the Zenodo repository or in https://github.com/johannessonlab/HetVPaper.
(iii) Sexual incompatibility: fertilization success of two crosses is reduced as a side-effect of the vegetative incompatibility interactions. The interaction between V and R was set to reduce fertilization success by 40%, which represents a mean of experimentally measured reduction in fertilization success, averaged over parental effects. The interaction of V and V1 also causes a reduction in fertilization success, which we vary from 0% to 40% in the simulations. Due to those two interactions, fertilization success of a RV1 x rV cross is reduced at least 40% due to R/V interaction, and up to 80% depending on the cost of the V/V1 interaction.

Structure of a simulation
Each population is composed of 1 000 individuals (2 000 genomes). The selfing rate of dikaryons is fixed for a given simulation. A population is started with uniform genetic background, either RV1 or rV. After 100 generations, 10 individuals (1% of the population) carrying the alternative genotype are introduced. This type of introduction could represent either the migration of individuals from a neighbouring population, or the second stage of an invasion subsequent to de novo mutations. We chose not to model explicitly mutation rates as we have no information from the natural system on what such rates could be, and introducing the new genotypes by mutations considerably increases simulation runtime. 300 generations after the introduction of the invading genotype, genotype frequencies are sampled every 20 generations over 100 generations and averaged. We study the effect of the following parameters on the outcome of the simulation (genotype frequencies after invasion): the rate of selfing of the population, the intensity of V/V1 balancing selection, the intensity of V/R balancing selection and finally the level of prezygotic isolation due to the V/V1 interaction. For each unique parameter combination, 100 replicated simulations are run. The results present the median and distribution of genotype frequencies across the 100 replicated simulations for each parameter combination. Fig. 6 and Supplementary Fig. 17 show the final genotype frequencies after invasion of the rV genotype in a population of RV1 background. Supplementary Fig. 18 show the reverse invasion scenario. We are primarily interested in the conditions that result in a coexistence of the rV and RV1 genotypes at intermediate frequencies, with the rV genotype at low frequency. The RV genotype is lethal and therefore absent from the population.

Additional remarks
Looking at Supplementary Figs. 17 and 18, we can see that the coexistence of the rV and RV1 genotypes is favoured by a high selfing rate. In fact, if selfing rate is low, the invading genotype is unable to enter the population. This is quite intuitive: because the RV combination is lethal, an invader carrying the V allele in an outcrossing population where all individuals carry R has no chance to invade, and vice versa with R invading V. If selfing rate is high enough however, selective advantage to the rare allele due to vegetative incompatibility kicks in and helps the invasion of the new genotype up to intermediate frequencies.
The selfing rate, strength of balancing selection (V/R and V/V1) and the strength of prezygotic isolation then all interact to produce the final genotype frequencies after invasion.
We can see that both the V/R and V/V1 balancing selections favour the invasion of the new genotype and the maintenance of rV and RV1 at intermediate frequencies.
The effect of V/V1 balancing selection is particularly visible in Fig. 6 in the main text, selfing rate 75%. In some cases, invasion is successful but the two dominant genotypes are rV and rV1. This may occur if the strength of the V/V1 balancing selection is strong enough to overpower the V/R balancing selection. An increased selfing rate, or increased strength of prezygotic isolation shifts the situation to an rV and RV1 coexistence (see for example Supplementary Fig. 17a, selfing 75% and compare the different outcomes with increasing strength of prezygotic isolation).
In conclusion, coexistence of the rV and RV1 genotypes is made possible by high selfing rates and R/V balancing selection. In addition, V/V1 balancing selection and prezygotic isolation favor that scenario, with prezygotic isolation and selfing contributing to lowered frequencies of the intermediate rV1 genotype.

Supplementary Figure 1.
Pleiotropic effects on the sexual function of genetically identified het genes in Podospora anserina. The diagram gives the two allelic categories (centre) of each of the nine identified het loci of P. anserina (top). The het loci can be classified as allelic or non-allelic systems (het-v is involved in two systems). Vegetative and sexual (in)compatibility reactions are given by a black (non-pleiotropic) or grey (pleiotropic) half arrow, whose directions go from the male parent to the female parent. For instance, the diagram specifies that a male Z2 x female Z1 cross is sterile while the reciprocal cross is fertile. The skull symbol indicates that a male S x female s cross can lead to meiotic drive of s (killing of S spores). For all systems marked as non-allelic, crosses between incompatible genotypes also lead to hybrid lethality, in the form of self-incompatible progeny. Note that het-c, het-d and het-e loci are multi-allelic and that a simplified nomenclature is used here where incompatible interaction are denoted with capital letters, while small letters are nonreactive (neutral) alleles. Based on refs. 2,3 .

Allelic systems
Non-allelic systems Vegetative incompatibility and sexual compatibility Vegetative incompatibility and sexual incompatibility Spore-killing

Supplementary Figure 2. (above)
Sliding window analysis (10kb-long with steps of 1kb) of representative chromosomes 2, 4, 5, 6 and 7 with values of genetic diversity (measured as either the pairwise nucleotide diversity π or as Watterson's theta θW) and the Tajima's D statistic. Relevant loci are marked, namely het genes, the centromere, and meiotic drivers of the Spok family. Note that Spok3 and/or Spok4 can be found at different locations in the genome depending on the strain.

Supplementary Figure 3.
Linkage disequilibrium decay in the Wageningen collection of P. anserina. The r 2 statistic was estimated from within 30 windows (50 kb long) randomly sampled along each chromosome. Hexagonal bins are coloured according to the log10 count of measurements for a given distance vs r 2 combination. The light blue curve corresponds to a nonlinear regression model following Remington et al. 4 , while the red line is a generalized additive mode smoothing. Red points correspond to mean r 2 values of 1 kb distance bins.

Supplementary Figure 4.
Changes of allele frequencies through time for het-q (data shown from years with more than five samples). Het−q (chromosome 7) Samples with mating data are coloured green or lilac based on their PCoA grouping of mating success (Fig. 2a). Gray points have no mating success data.

Supplementary Figure 6.
Genetic differentiation between the two mating groups. The Fst statistic was calculated in sliding windows (10kb-long with steps of 1kb) and plotted in dark blue. The light red areas represent the maximum Fst values reached by 1000 permutations of group membership. Marked loci include het genes, the centromere, the mating type locus (MAT), and meiotic drivers of the Spok family. Note that Spok3 and/or Spok4 can be found at different locations in the genome depending on the strain.

Supplementary Figure 7.
The pairwise nucleotide diversity between mating groups (Dxy) in sliding windows of 10 kb (1 kb steps). Relevant loci are marked, namely het genes, the centromere, and meiotic drivers of the Spok family. Note that Spok3 and/or Spok4 can be found at different locations in the genome depending on the strain.   Figure 10.
The region encompassing Pa_5_12710 (het-Vb) and Pa_5_12720 (het-Va) is responsible for low fertility of rV x RV1 crosses. Crosses were set up between the given genotypes by fertilization (also known as spermatization). In the top row, the haploid RV1 strain is used as male (m.) parent to fertilize a haploid female (f.) parent. In the bottom row, the RV1 strain is used as female parent. The white arrows point to rare fruiting bodies. In both directions, the Δ12810-12690 deletion restores full fertility. Reintroduction of the region encompassing the Pa_5_12710 and Pa_5_12720 genes in the deletion strain reduces fertility, indicating that this region is responsible for the sexual incompatibility. The reduction on fertility is stronger when the female parent is of RV1 genotype.  Patterns of segregation of het-r and het-v. In P. anserina the products of meiosis are duplicated by a mitosis event, and then rearranged in paired non-sister nuclei within two pairs of identical spores per ascus (only two shown for simplicity). A given locus might undergo first or second division segregation (FDS and SDS) during meiosis with a frequency dependent on the distance to the centromere. As het-r and het-v are on different chromosomes, their segregation is independent, which would result in equal proportions of the resulting spore genotypes if there was no incompatibility (0.25 each). As the R and V alleles are incompatible, a fraction of the spores dies upon germination, but the exact proportion is unknown. If a dikaryotic spore contains one RV nucleus, the incompatibility reaction might also kill the other nucleus during germination (red lines), in which case only 47.5% of all the nuclei produced by the cross survive. Alternatively, the other nucleus might escape by sectoring into a haploid mycelium while the RV counterpart dies (75% survival). As the V product is diffusible, it might still kill the sectoring RV1 nuclei (purple lines). That would result in the survival of 65% of nuclei per cross. Expression of het-Vb alone has the same phenotype as the absence of both genes (same phenotype as V1 recipient). In contrast, expression of het-Va alone abolishes incompatibility to V. Note that the two mutants of het-Va (in grey and light green) have different phenotypes. While het-Va Y233A behaves as null mutant (with respect to allelic V/V1 incompatibility), het-Va Y233F has a wild-type behaviour, in the sense that this allele abolishes incompatibility to V, as does wild-type het-Va. Yet, in contrast to wild-type, joined expression of het-Va Y233F with het-Vb does not confer incompatibility to V1. The phenotype conferred by het-Va Y233F is thus different both from wild-type and from het-Va Y233A. The activity of the HET-Va Y233F allele could be explained by a potential dominant negative effect on the methyltransferase activity of wild-type het-Va or by a residual activity of this more conservative change. (b) A possible mechanistic model for V/V1 incompatibility based on the experimental results presented in (a) and on the predicted functions of the HET-Va and HET-Vb protein domains is given. In this model, HET-Va SET domain protein deposits a methylation mark on an unknown cellular target. This target is thus marked as "self" in V strains and not recognized by the HET-Vb reader and cell death inducing protein.

Supplementary
In V1 strains, both proteins are absent and the target is not marked. Upon fusion of V and V1 strains, unmarked targets are detected as "non-self" by the HET-Vb reader which induces cell death. The fact that expression of HET-Vb alone in V1 does not lead to incompatibility rules out a simple scenario where HET-Vb directly recognizes the unmarked target. Rather, presence of HET-Va is also required to bring about cell death. This could occur for instance if HET-Vb recognizes a methylation state formed transiently when naive (V1-type target) is exposed to HET-Va, as depicted in the graphic model. In the absence of population data of all species, it remains unclear if the presence of het-v is ancestral to the species complex. However, the RI groups in P. anserina cluster together for all genes, confirming that the divergence beyond the het-v locus happened after P. anserina diverged from other species. Branch support (bootstraps) is indicated above the tree branches. Branch lengths are proportional to the scale bar (nucleotide substitution per site). Rooting was done arbitrarily with CBS112942.