Repeat-Induced Point Mutation and Gene Conversion Coinciding with Heterochromatin Shape the Genome of a Plant-Pathogenic Fungus

ABSTRACT Meiosis is associated with genetic changes in the genome—via recombination, gene conversion, and mutations. The occurrence of gene conversion and mutations during meiosis may further be influenced by the chromatin conformation, similar to the effect of the chromatin conformation on the mitotic mutation rate. To date, however, the exact distribution and type of meiosis-associated changes and the role of the chromatin conformation in this context are largely unexplored. Here, we determine recombination, gene conversion, and de novo mutations using whole-genome sequencing of all meiotic products of 23 individual meioses in Zymoseptoria tritici, an important pathogen of wheat. We confirm a high genome-wide recombination rate of 65 centimorgan (cM)/Mb and see higher recombination rates on the accessory compared to core chromosomes. A substantial fraction of 0.16% of all polymorphic markers was affected by gene conversions, showing a weak GC-bias and occurring at higher frequency in regions of constitutive heterochromatin, indicated by the histone modification H3K9me3. The de novo mutation rate associated with meiosis was approximately three orders of magnitude higher than the corresponding mitotic mutation rate. Importantly, repeat-induced point mutation (RIP), a fungal defense mechanism against duplicated sequences, is active in Z. tritici and responsible for the majority of these de novo meiotic mutations. Our results indicate that the genetic changes associated with meiosis are a major source of variability in the genome of an important plant pathogen and shape its evolutionary trajectory.

some fungi as well as in some plant and algae species (18). Although chromatin configuration and hence the histone modifications are assumed to affect the rate of gene conversions, such an association was so far not identified.
Finally, meiosis is also associated with mutations that occur before or during meiosis (38)(39)(40). Since mutations on average are considered to be deleterious, mutation rates, in general, are low but can also differ greatly between species (41). Meiosis-associated mutation rates, in turn, can differ greatly from the corresponding mitotic mutation rate in the same species, indicating different mechanisms and/or constraints. For example, the germline mutation rate in humans and mice is 1.2 Â 10 28 or 5.7 Â 10 29 per nucleotide per generation, respectively, two orders of magnitude lower than the corresponding mitotic mutation rate (10,42,43). Here, germline mutations might be linked to DSBs and their repair (44), with a higher number of mutations occurring in the vicinity of recombination events (45). Estimates of meiosis-associated mutation rates in different fungal species also vary considerably. In S. cerevisiae, the meiosis-associated mutation rate is 8 Â 10 28 per bp per cell generation (44), much higher than the mitotic mutation rate of 3.3 Â 10 210 per bp per cell division (46). In Neurospora crassa, the meiosis-associated mutation rate is very high at 3.38 Â 10 26 per bp per generation (47), contrasting with a much lower mitotic mutation rate of 6.7 Â 10 210 per bp per cell division (48). Interestingly, this extremely high meiotic mutation rate in N. crassa is caused by repeatinduced point mutation (RIP), a fungal defense mechanism against duplicated sequences ( Fig. S1A) (40,49). RIP is restricted to haploid parental nuclei just prior to karyogamy and meiosis and acts on duplicated sequences of a minimum length of 400 bp in N. crassa (50)(51)(52). Once recognized, both duplicated sequences will be mutated in a C!T manner (50)(51)(52), and RIP can sometimes leak into adjacent non-repetitive regions (53,54). RIP signatures have been detected in the genomic sequences of many fungi; however, active RIP was experimentally confirmed in only a few fungal species (47,55). Hence, there is a growing body of evidence that suggests that the mutational processes prior to and during meiosis differ from those during mitosis. In fungi, the meiotic mutation rate appears to be higher than the mitotic mutation rate, which for some fungi is assumed to be the result of RIP.
Here, we use tetrad analysis to determine genetic changes associated with meiosis in the ascomycete fungus Zymoseptoria tritici, a pathogen of wheat. The haploid genome of Z. tritici comprises 13 core chromosomes and a set of smaller accessory chromosomes (56). These non-essential accessory chromosomes carry a fitness cost (57), are enriched in the facultative heterochromatin mark H3K27me3 (58), have a higher mutation rate during mitosis (59), and exhibit a meiotic drive (60). The availability of complete tetrads for Z. tritici allows us here to address all three major classes of genetic changes associated with sexual reproduction. In particular, the frequency and distribution of mutations associated with sexual reproduction are unknown, although Z. tritici has an asexual and a sexual reproductive cycle, with the latter being the main source of the primary inoculum during the initial stages of the infection (61). The mitotic mutation rate in Z. tritici has been determined experimentally by mutation accumulation experiments at 3.2 Â 10 210 per bp per mitotic cell division (59), which is similar in other fungi (46,48). Although histone modifications affect the mitotic mutation rate of Z. tritici (59,62), it is unknown if the distribution of meiotic recombination events, gene conversion, and meiosis-associated mutations are also influenced by these histone modifications. Finally, although the genome of Z. tritici shows signatures of RIP (56), RIP has so far not been demonstrated experimentally, and the efficacy of this mutational mechanism in this pathogen is not known. Given the fact that 18.6% of the genome of Z. tritici is represented by transposable elements (TEs), it is plausible that RIP is less efficient in Z. tritici than in N. crassa or fails to recognize some duplicated regions (63). Here, we study all major classes of meiosis-associated genetic changes in Z. tritici by analyzing whole-genome sequences of complete tetrads to (i) estimate recombination and gene conversion rates for core and accessory chromosomes, (ii) determine the association between recombination and gene conversions with chromatin modifications, and (iii) estimate meiotic mutation rates in Z. tritici. The use of tetrads allowed us to detect and describe the effects of active RIP and generate a fine-scale map of recombination and gene conversion events and its association with chromatin modifications.

RESULTS
Accessory chromosomes show higher recombination rates. To determine the distribution of recombination events during meiosis in Z. tritici, we used previously published tetrads and obtained whole-genome sequences for the tetrad progenies (60). To this end, we included 23 tetrads comprising four ascospore isolates, totaling 92 genomes. An average of 118,772 single nucleotide polymorphisms (SNPs; 0.3% of all analyzed genomic sites) per tetrad was used for the analysis (see Materials and Methods), with the SNP density being similar between core and accessory chromosomes (Fig. S1B). The median distance between SNPs was 61 bp, with few instances of distances between SNPs exceeding 5,000 bp (Fig. S1C). From these data, we identified individual recombination events with the CrossOver tool from the ReCombine package (64) and calculated the recombination rate. A total of 1,138 crossover events were observed, resulting in a genome-wide recombination rate of 65 centimorgan (cM)/Mb, consistent with previously published estimates of the recombination rate in this fungus. Intriguingly, the recombination rate was significantly higher on accessory chromosomes (92.7 cM/Mb) than on the core chromosomes (62.6 cM/Mb) (Fig. 1A, Fig. S1D, Table S1A and B). The recombination rate was negatively correlated with the chromosome length (Pearson's R = -0.76, P = 0.00017) (Fig. 1B), but the absolute number of crossovers per chromosome was positively correlated with the chromosome length (Pearson's R = 0.96, P = 6.2 Â 10 211 ) (Fig. S1E). When two or more crossovers are Distribution of distances between adjacent CO events (CO-CO distance) in the entire genome (pink), on core chromosomes (blue), and on accessory chromosomes (yellow). The red line represents the fitted gamma distribution for the observed CO-CO distances. g values higher than 1 indicate positive CO interference, i.e., a lower than expected number of small CO-CO distances. (D) Gene conversion rates for the genome and core and accessory chromosomes. (E) Violin plot of non-crossover-associated gene conversion (NCO-GCs) and crossover-associated gene conversion (CO-GCs) tract lengths for core and accessory chromosomes. The tract lengths of the gene conversion events spanning TEs are excluded (different letters depict significantly different groups with a P value of ,0.05 in the Wilcoxon test with Bonferroni correction). In panels A and D, P values of paired Wilcoxon tests are shown (*, P , 0.05; **, P , 0.005; ***, P , 0.0005). present in a bivalent, they tend not to occur near each other in many species-a process called crossover interference (65). We accessed the distribution of the distances between COs using the gamma distribution (66)(67)(68) and detected crossover interference (g . 1) for both core and accessory chromosomes (see Fig. 1C). When correlating the crossover frequencies with regions enriched in heterochromatin marks, we found that regions enriched in heterochromatin marks (H3K27me3 or H3K9me3) were associated with higher crossover frequencies, whereas the euchromatin mark H3K4me2 did not show a higher recombination rate than regions lacking all three marks ( Fig. 2A).
To further assess the variation in recombination rate across the genome, we estimated the recombination rates in 20 kb non-overlapping windows. We defined hotspots as those 20 kb windows that showed a significantly (P , 0.001) higher number of recombination events than expected as determined by the Poisson distribution and identified a total of 52 recombination hotspots with more than 4 crossover events per 20 kb window. Recombination rates per window were highly variable, ranging from 0 to 1,196 cM/Mb (Fig. S2). We conclude that recombination varies considerably across the genome, occurs in hotspots, and is higher on the accessory than on the core chromosomes. A more detailed analysis of the immediate vicinity of the crossover events allowed us to identify four motifs that were significantly overrepresented in the 1 kb sequences surrounding the crossover and to correlate recombination rates with short sequence repeats (SSRs) ( Fig. S3C and D). However, a more fine-scale analysis of the nine 1 kb windows with $4 crossover events showed no apparent correlation between the location of TEs, SSRs, genes, and the three histone modifications, H3K4me2, H3K9me3, and H3K27me3 (Fig. S3A). We considered the relevance of high recombination rates with regard to gene evolution and correlated the recombination map with the coordinates of protein-coding genes. In total, we observed 376 genes in the CO hotspots. From the 376 genes observed in the CO hotspots, 6 encode effector candidates and 9 encode CAZymes (Table S1C), which are gene categories with a putative relevance in the pathogenicity of the fungus.
High gene conversion rates are more uniformly distributed within the genome. The genomes of ascospore progenies resulting from one meiotic event provided us with a unique opportunity to characterize the location and distribution of gene conversion events. We therefore identified gene conversion events along the genome and compared the gene conversion rate for different genomic compartments. We identified a total of 890 gene conversion events, 712 associated with non-crossover events (non-crossover gene conversion [NCO-GC], 80%) and 178 associated with crossover events (crossover gene conversion [CO-GC], 20%) ( Table S2A). Each of these gene conversion events was associated with at least three converted markers (SNPs), with the total number of converted markers being 3,529 and 979 for NCO-GC and CO-GC, respectively. We distinguished the gene conversion events on core and accessory chromosomes of Z. tritici and found that 755 gene conversion events were located on core and 135 on accessory chromosomes (Table S2B). We next explored the general patterns of converted SNPs. We observed a weak GC-bias in the gene conversions (binomial test, P = 0.0215; Fig. 2C). Based on genome-wide SNPs identified among the ascospore isolates, we found that the genome-wide gene conversion rate was 1.6 Â 10 23 per SNP (Fig. 1D, Table S2B). The gene conversion rate was significantly higher on the accessory chromosomes than on the core chromosomes (paired Wilcoxon signed rank test, P = 1.0 Â 10 23 ) (Fig. 1D, Fig. 2D). Tract length is considered to be one of the main determinants of gene conversion rates; therefore, we estimated the median conversion tract length for both types of gene conversion events. The median tract length for non-crossover (NCO-GC) was 539 bp, and 432 bp for gene conversions associated with crossovers (CO-GC), respectively (Fig. 1E). Although the tract length on the accessory chromosomes appeared to differ from the corresponding tract length on the core chromosomes, none of these comparisons were significant (e.g., NCO-GC accessory versus core chromosomes P = 0.927) (Fig. 1E).
In the same way as we characterized the distribution of recombination hotspots along the fungal genomes, we used the map of gene conversion events to identify regions with exceptionally high gene conversion rates, here defined as gene conversion hotspots. First, to assess the distribution of gene conversion, we divided the genome into 20 kb non-overlapping windows and calculated the number of gene conversion events in each window. We identified 32 gene conversion hotspots with more than 4 gene conversion events per 20 kb window that were significantly different from the background with P value of , 0.001 defined by the Poisson distribution ( Fig. S4A to C). Of the 145 genes in the gene conversion hotspots, 2 genes encode for effector candidates and seven genes encode for CAZymes (Table S2C). We speculate that the increased rate of gene conversion in these putative pathogenicity-related genes could be a putative mechanism of rapid evolution. However, the majority of the genes associated with gene conversion hotspots (95 genes) encode proteins that have not yet been functionally annotated, preventing any meaningful enrichment analysis. Taken together, the ascospore population allowed us to precisely map and characterize gene conversion events in Z. tritici. We found that rates of gene conversion showed less variation across the genome than recombination rates, but with both being higher on accessory chromosomes than on core chromosomes.
Gene conversion rates are higher in regions enriched in heterochromatin modifications. Previous studies in chicken B-cell lines have identified a correlation of chromatin structure with gene conversion (69), but detailed analysis of the effect of specific histone modifications is mostly missing. Thus, to investigate the potential effect of histone modifications on the rate of gene conversion in Z. tritici, we conducted a correlation analysis of maps of histone modifications and gene conversion rates. We focused on three histone marks which have been previously well characterized in Z. tritici using chromatin immunoprecipitation of antibodies targeting specific histone modifications followed by sequencing (ChIP-seq): the euchromatin mark H3K4me2, the constitutive heterochromatin mark H3K9me3, and the facultative heterochromatin mark H3K27me3 (58) (Fig. 2C, Table S2D). Importantly, we excluded TEs from the full-factorial analysis of gene conversions, and therefore, regions enriched with H3K9me3 and H3K27me3 in this analysis are not associated with TEs. Our analyses showed that higher gene conversion rates associate with regions enriched in heterochromatin marks, particularly in regions solely enriched with H3K9me3, as well as in regions enriched with both H3K9me3 and H3K27me3 (Fig. 2B). Regions enriched solely in H3K27me3 were not associated with higher gene conversion rates, which contrasts with their correlation with the highest crossover frequencies ( Fig. 2A). The median gene conversion rate in the regions enriched with H3K9me3 was 3 Â 10 22 per converted marker per tetrad per meiosis, and the median gene conversion rate in the H3K9me3 and H3K27me3 regions was 1.3 Â 10 22 per converted marker per tetrad per meiosis. Genome regions not enriched with either of the three histone marks, H3K4me2, H3K9me3, and H3K27me3, showed the lowest gene conversion rate of 0.9 Â 10 23 . In summary, the histone maps allowed us to reveal that repressive heterochromatin modifications, especially H3K9me3 but not H3K27me3, are associated with higher gene conversion rates in Z. tritici.
The majority of meiotic mutations are caused by RIP. Recombination and gene conversion as well as meiosis and meiosis-associated processes themselves are considered to be mutagenic (38)(39)(40)45). Hence, we asked to what extent de novo mutations had occurred in the ascospore progenies. These mutations are distinguished by being absent in both of the parental strains but present in the ascospores of the tetrad. Meiotic mutation rates often differ from mitotic mutation rates and are associated with recombination or premeiotic processes (39,(42)(43)(44)(45)(46). We observed a total of 526 de novo mutations that were absent in both of the parental strains for the ascospores of the 23 tetrads. Hence, through the sexual cycle, there were, on average, 22.9 mutations per genome per generation, resulting in a mutation rate of 5.7 Â 10 27 mutations per bp per generation (Fig. 3A, Table S3A). Of these mutations that originated from the sexual cycle, 242 (42%) resided in one particular region that clearly stood out as we mapped meiotic SNPs along the genome. This high number of SNPs were located in a 14,001 bp region on chromosome 3 (Fig. 3A to C). Each of the 242 mutations in this region on chromosome 3 was a CG:TA transition. As we further inspected this region, we noted that the 14 kb region spanning 1,668,370 bp to 1,682,371 bp on chromosome 3 showed an increased sequencing coverage corresponding to a duplication of the region in the parent IPO94269 (Fig. 3B). This means that this duplication on chromosome 3 was present in all 23 meioses and could possibly be identified by RIP, a genome defense mechanism against such duplicated regions (40). A high number of transitions in a duplicated region is the hallmark of such active RIP. Indeed, the high number of mutations in the duplicated region and the fact that all of these mutations are CG:TA transitions indicates active RIP in Z. tritici. Interestingly, the duplicated region on chromosome 3 showed RIP mutations in only 10 of the 23 tetrads (Fig. S5A). Further examination of the sequencing coverage depth revealed that, despite being present in the parental strain IPO94269, the duplicated region was lost in 11 of the tetrads. In these tetrads all 4 ascospores had a normalized sequencing coverage of the duplicated region of approximately 1 (Table S3B), while in 12 tetrads the duplication was still present (normalized sequencing coverage of approximately 2) and therefore followed Mendelian segregation. In 10 of these 12 tetrads that still had the duplication, RIP mutations were detectable in this region, while in the remaining 2 tetrads no RIP mutations occurred in this region (tetrads A03-9 and S1C4_A2) (Fig. S5A). The duplicated region contains six predicted genes, one of which has a functional annotation describing it as similar to transferase hexapeptide domain protein. Therefore, the putative function of the genes affected by the RIP of the duplication cannot yet be deduced.
We further addressed to what extent de novo mutations could be the result of meiosis-associated RIP. Here, we found that 172 of the meiotic de novo mutations (32% of all de novo mutations) were detected in transposable elements (Fig. 3C), and 166 of these 172 mutations (96.5%) were CG:TA transitions, the putative signature of RIP. To examine whether these mutations were also most likely caused by RIP, we determined which class of transposable elements was affected. Of the 172 de novo mutations in TEs, 152 (88.3%) were located in class I transposons, which replicate via an RNA intermediate and are also referred to as copy-and-paste transposons (70)(71)(72). Of these 152 mutations in class I transposons (long interspersed nuclear elements [LINEs] and long terminal repeats [LTRs]), 146 were CG:TA transitions, with 5 TA:GC transitions and 1 G!C transversion, while only 11 CG:TA transitions occurred in class II transposable elements (Fig. 3D). Class II TEs are referred to as cut-and-paste transposons, which are excised and moved to new locations in the genome and hence do not create repeated sequences (71). The mutation rate in class I (copy-and-paste) TEs is 1. Investigations in a few other ascomycete fungi have provided evidence that RIP can affect the vicinity of duplicated regions by "leakage" of mutations into these non-repeated regions (53,54). Based on these previous studies, we here also asked if RIP mutations would be present in regions adjacent to TEs. We subdivided each TE into 40 equally sized windows (average window size, 71 bp) and the two adjacent regions into 20 windows each (window size, 71 bp) and counted the number of de novo mutations in each of these windows. This approach allowed us to show that RIP mutations were not equally distributed along TE sequences. Within TEs the most distal windows showed the highest number of de novo mutations, whereas the more central sequences were less likely to be mutated (Fig. 3E). Further analysis of the effect of the different classes of transposable elements showed that the highest number of RIP mutations in the distal regions of the TEs was mainly due to LTR transposons (Fig. S5B), whereas LINE and class II transposable elements did not show such a variation of the number of mutations along the length of the transposable element (Fig. S5B). In regions adjacent to the TEs a much lower number of mutations occurred than occurred in the TE sequences. Overall, the mutation rate in the regions adjacent to TEs was 7.4 Â 10 28 per bp per generation, which was lower than the genomewide mutation rate of 1.5 Â 10 27 per bp per generation (for regions excluding TEs and excluding the duplication). It therefore appears that leakage has not occurred in the vicinity of the RIP-mutated TEs. Only 21.3% (112 mutations) of all de novo mutations were located in the regions outside the transposable elements and the 14 kb region on chromosome 3. These de novo mutations in non-TE and non-duplicated regions amounted to a significantly higher proportion of TA:CG transitions and transversions (22.3%) (Fig. 3C).
The segregation pattern of the de novo mutations varied between the different compartments, with those in TEs showing a higher proportion of 1:3 segregation (Fig. S5C) than those in the 14 kb duplication or the mutations outside TEs. The segregation pattern of de novo mutations can indicate at what stage a DNA lesion or mismatch may have occurred and at what stage this lesion or mismatch (possibly caused by RIP) was resolved. A 2:2 segregation indicates that the mismatch was resolved prior to the replication cycle of the meiosis, while a 1:3 segregation indicates that the mismatch was resolved only during replication or later. Hence, the higher proportion of 1:3 segregation patterns observed in the TEs could indicate a delayed resolution of the mismatches introduced by RIP in the heterochromatic TEs compared to the resolution in other genomic compartments (Fig. S1A). The 39 neighbor of the targeted cytosine has a strong influence on RIP in N. crassa, with a strong preference for CpA (73,74). We therefore identified the dinucleotide context of the 242 mutations in the 14 kb duplicated regions. Of these, 240 were in a CpA context, and the remaining two were in CpT context. The rate at which a single CpA was mutated during RIP was therefore 0.54% (Fig. S5D). The mutation rate varied along the 14 kb duplication, but this did not appear to correlate with the proportion of dinucleotides in the sequence that were CpA (Fig. S5E). The number of individual mutations that were induced within the duplicated region varied considerably between tetrads (Fig. S5A) but also between the individual ascospores of the tetrads (Fig. S5F). For the 24 ascospores containing the duplication on chromosome 3, the number of mutations within the duplication in individual ascospores varied from 0 to a maximum of 57 (median 11), with most ascospores showing 0 to 10 mutations (Fig. S5F). In summary, RIP has been proposed to be an important player in the genome evolution of Z. tritici. However, this is the first experimental evidence for active RIP in this fungus, and we observed that the large duplication was not always mutated by the RIP mechanism and that the total number of mutations introduced by RIP was low-leaving many cytosines in the duplicated regions unaffected-which together could indicate a low efficiency of RIP in Z. tritici.

DISCUSSION
Here, we used tetrad analysis to estimate recombination rates, gene conversion rates, and de novo mutation rates associated with meiosis from 23 individual meioses in Z. tritici. The ability to dissect genetic events in individual ascospore progenies isolated from tetrads provided us with highly precise maps of meiosis-associated changes along the fungal genome. Our results show (i) higher gene conversion rates and recombination rates on accessory chromosomes compared to core chromosomes, (ii) a correlation of recombination rates and gene conversion rates with histone marks associated with heterochromatin, and (iii) elevated de novo mutation rates during sexual reproduction caused by active RIP .
Our dissection of recombination events during individual meiotic events allows us to confirm the previously reported high recombination rates in Z. tritici-on average, 65 cM/Mb. Interestingly, 35 of the 52 recombination hotspots detected in this study overlapped with previously reported recombination hotspots in Z. tritici, which used distinct isolates from Switzerland, indicating that recombination hotspots may be determined by certain domains or conserved marks along the genome (17). However, these are most likely not small sequence motifs, as there is little similarity between those identified within our study and those from the previous study with Swiss isolates (17). In contrast to previous studies using population genomic data (18), we show that the rate of recombination is, in fact, higher on accessory chromosomes than on core chromosomes (92.7 cM/Mb and 62.6 cM/Mb, respectively) and is negatively correlated with chromosome size. Hence, our results are in line with the general observation that smaller chromosomes tend to have higher recombination rates (75). Earlier population-based recombination rate studies in Z. tritici showed, however, lower recombination rates computed as rho for accessory chromosomes (rho = 0.001) than on core chromosomes (rho = 0.024) (18). Rho is the product of the actual recombination rate and the effective population size Ne, and we speculate that the discrepancy in recombination rate measures may reflect the lower effective population size of the accessory chromosomes compared to the core chromosomes (18). Intriguingly, we could not identify any crossovers on an accessory chromosome in nine instances, despite the presence of the respective homologous chromosomes. A minimum of one crossover per homologous chromosome pair is considered to be required for proper segregation of the homologous chromosomes (4,76). Indeed, in four of these nine instances, we observed such segregation errors. However, in the remaining five instances, the accessory chromosomes properly segregated despite the absence of crossovers. Hence, we conclude that crossovers are not essential for proper chromosome segregation in Z. tritici. spo11 deletion mutants in S. cerevisiae and Sordaria macrospora were previously used to observe the consequences of the absence of DSBs and recombination on the segregation of chromosomes during meiosis. In these fungal species the absence of DSBs and recombination caused widespread chromosome segregation errors, highlighting the importance of recombination for proper segregation (77,78). The relatively high frequency of properly segregated, non-recombined chromosomes in Z. tritici indicates the presence of a non-recombination-dependent segregation system which might also be involved in the meiotic drive system of the accessory chromosomes in Z. tritici (60).
We predict that high levels of gene conversions in Z. tritici can affect allele frequencies to a greater extent than recombination and thereby shape genome composition. The genome-wide gene conversion rate of 1.6 Â 10 23 per SNP identified in our study is approximately 20 times lower than the genome-wide gene conversion rate in S. cerevisiae (3.8 Â 10 22 per SNP) (32,33), but approximately an order of magnitude higher than the genome-wide gene conversion rate in N. crassa (0.7 Â 10 24 to 2.2 Â 10 24 per SNP) (32,33). The variation in the gene conversion rate between species might be influenced by tract length and recombination rate since gene conversion is positively correlated with both characteristics (33,34,36). Our data confirm this since Z. tritici has shorter NCO-GC and CO-GC tract lengths and lower recombination rates (539 bp, 432 bp, and 65 cM/Mb, respectively) than S. cerevisiae (1,681 bp, 1,841 bp, and 375 cM/ Mb) (33). N. crassa in turn shows longer NCO-GC tract lengths but shorter CO-GC tract lengths and a lower recombination rate (950 bp, 284 bp, and 20 cM/Mb, respectively) than Z. tritici (33). Similar to recombination, accessory chromosomes have higher gene conversion rates than the core chromosomes. The smaller size of accessory chromosomes in contrast to core chromosomes could influence the gene conversion rate on accessory chromosomes since smaller chromosomes tend to have higher rates of gene conversion (33).
We found a strong association between chromatin modifications and recombination and gene conversion rates. We speculate that this might be due to two non-exclusive mechanisms: the effect of histone modifications on the location of DSBs and/or their effect on DSB and heteroduplex repair. Although a direct correlation between DSB formation and H3K4me3 has been shown in S. cerevisiae (20), the underlying causality of this correlation and the extent to which it applies to other organisms is unclear (79,80). Local chromatin conformation can influence how DSBs and heteroduplex DNA are repaired, thereby influencing the later stages leading to recombination and gene conversion (81)(82)(83). In contrast to S. cerevisiae, we see an association of heterochromatin as marked by H3K9me3 and H3K27me3 with recombination in Z. tritici. Similarly, and again in contrast to the scenario in yeast, we observe that gene conversion rates are highest in heterochromatin regions in Z. tritici. This suggests that chromatin conformation may affect gene conversion primarily through its effect on DSB and heteroduplex repair and not via its impact on accessibility. DSBs can be resolved by synthesis-dependent strand annealing (SDSA), leading to a homologous recombination-mediated pathway and hence NCO events (84,85). Our results show increased gene conversion rates in regions enriched in H3K9me3. Several studies suggest that H3K9 di-or trimethylated (H3K9me2/3) heterochromatin promotes homologous recombination (HR) (85)(86)(87)(88)(89). DNA double-strand breaks promote chromatin stabilization by H3K9me3 and activate DSB signaling proteins (90). The role of H3K9me3 in promoting homologous recombination could mean that the correlation between H3K9me3 and the gene conversion rate we observed in our study is caused by the effect of H3K9me3 on recombination. However, 80% of the gene conversion events we detected are associated with non-crossover events (NCO-GC). This suggests that the correlation between H3K9me3 and gene conversions is caused not by recombination, but possibly by the effect of histone modifications on the repair of DSBs and heteroduplex DNA. In contrast to DSBs in regions enriched in H3K9me3, DSBs in regions enriched in H3K27me3 were found to be frequently repaired by microhomologymediated end joining (MMEJ), a nonhomologous repair pathway that does not promote homologous recombination and therefore does not lead to gene conversion or crossover (91,92). We find support for a similar effect of chromatin modifications on DSB repair in our study. Indeed, increased repair of DSBs via a nonhomologous pathway in the H3K27me3-enriched regions could explain the reduced gene conversion rate in the regions enriched in both the H3K9me3 and H3K27me3 regions (1.3 Â 10 23 per SNP) compared to regions enriched in H3K9me3 only (2.8 Â 10 23 per SNP) as observed here. Similarly, the absence of an increased gene conversion rate despite the high frequency of crossovers in the regions enriched in H3K27me3 alone could be explained by MMEJ being favored by H3K27me3. However, it is unclear why H3K27me3 is associated with the highest crossover frequencies, since MMEJ should not promote homologous recombination and therefore should result in a lower crossover frequency. In wheat, high crossover rates were associated with H3K27me3, but the mechanism causing this correlation is still unclear (93). We speculate that the observed correlation between H3K27me3 and higher crossover frequencies may be driven by the enrichment of this mark on accessory chromosomes, but functional validation is still missing. In conclusion, we see indications that histone modifications could affect gene conversion rates mainly via their effect on the DSB and heteroduplex repair.
De novo mutations associated with meiosis occurred at a rate of 5.7 Â 10 27 per bp per generation in Z. tritici, which is approximately three orders of magnitude higher than the mitotic mutation rate (3.2 Â 10 210 per site per cell division) which we previously determined in a mutation accumulation experiment (59). De novo mutations were not uniformly distributed across the Z. tritici genome, ranging from the highest mutation rate of 7.5 Â 10 24 per bp per generation within the duplicated region on chromosome 3 to 1.5 Â 10 27 per bp per generation in genomic regions excluding the duplicated region and excluding TEs. Class I and II TEs showed a mutation rate of 1.3 Â 10 26 and 2.8 Â 10 27 per bp per generation, respectively, which is higher than the genome mutation rate but lower than the mutation rate in the duplicated region. This is probably due to the variability of TEs within these two groups, which also contain many non-duplicated sequences. Higher meiotic than mitotic mutation rates have been reported in other fungi such as S. cerevisiae and N. crassa (44,(46)(47)(48). In N. crassa, the difference between mutation rates during mitosis and meiosis is mostly due to RIP, a fungal defense mechanism againt duplicated DNA sequences that occurs in the haploid nuclei just prior to meiosis and that induces CG:TA transitions in duplicated sequences and transposable elements (47,50). In Z. tritici, signatures of past RIP have been found by analyses of genome data (56,94), and here we can confirm that RIP is an active mechanism in Z. tritici. This is evident as 77% of the de novo mutations associated with meiosis were located in class I transposable elements (copy-and-paste elements) and in the 14,001 bp long region on chromosome 3 that was duplicated in the IPO94269 parent. Interestingly, RIP in Z. tritici is not consistently efficient in mutating duplicated sequences, as we find that duplicated sequences in two of the tetrads were not mutated and the number of RIP mutations introduced even in duplicated regions is low. The duplication that was present in the parental strain IPO94269 appears to have been frequently lost during crosses. How and at what stage this loss of the duplication occurred is currently unclear. However, most likely, the loss must have occurred before the RIP process, as no RIP mutations were detected in this region. Until now, RIP was experimentally demonstrated in only a few fungal species (95-100), with lower efficiency in Leptosphaeria maculans and Podospora anserina than in the highly efficient RIP in N.
crassa (97,101). We consider that variation in RIP efficiency may be a more common phenomenon that reflects a trade-off between the evolutionary costs of the mutations introduced by RIP and the evolutionary costs of TEs.
We find that class I (copy-and-paste) TEs show a much higher RIP mutation rate than class II (cut-and-paste). This observation is expected, as the transposition mode of the class I TEs results in duplicated sequences that are recognized by the RIP machinery (70). Interestingly, we find that the RIP mutations are not equally distributed along the TE sequences. The most distal windows composing 5% of the length of the TEs amount to 35% of all de novo meiotic mutations that occurred within the TEs. This high frequency of de novo mutations at the ends of TEs was mostly due to the accumulation of mutations in the distal parts of the LTR transposons. LTR transposons contain terminal repeats of 200 bp to 500 bp in length and therefore already contain duplicated sequences in close vicinity (72). These duplicated terminal sequences in class I transposons seem to be prominent targets for RIP. In various other fungi, RIP leaks into adjacent regions from the duplicated sequences of TEs (53,102,103) and is thought to play a role in the rapid adaptive evolution of effector genes involved in host-pathogen interactions (53,102,103). However, we do not see such leakage from transposable elements in Z. tritici but, rather, see lower mutation rates in the vicinity of transposable elements than in the rest of the genome, which we cannot currently explain. Taking these data together, we show that RIP affects transposable elements in Z. tritici to various degrees and that the mutational environment thus influences the activity of these important drivers of genome evolution.
In conclusion, we show that recombination and gene conversion are correlated with histone modifications in different ways in Z. tritici and that RIP is active, albeit at a lower efficiency than in N. crassa, in this fungus, affecting duplicated sequences as well as TEs. As a result, meiotic mutation rates for Z. tritici are three orders of magnitude higher than the mitotic rates, demonstrating the major impact that genetic changes associated with meiosis have on the genome composition of this important plant pathogen.

MATERIALS AND METHODS
Fungal material. Tetrads used for the sequencing analysis were obtained from the study of Habig et al. (60) and include ascospores isolated from crosses between the Dutch isolates IPO94269 and IPO323 (available from the Westerdijk Institute, Utrecht, The Netherlands, with the accession numbers CBS115943 and CBS115941) and from the crosses between IPO94269 and whole-chromosome isogenic deletion strains (Dchr14 Dchr20) of the reference strain IPO323 generated in the study of Habig et al. (57). All ascospores were cultivated for DNA extraction at 18°C and 200 rpm in YMS (4 g/L yeast extract, 4 g/L malt, 4 g/L sucrose) medium for 5 to 7 days inoculated directly from 280°C glycerol stocks.
Genome sequencing and data analysis. For sequencing, DNA of 84 ascospores was isolated using a phenol-chloroform extraction protocol as described previously (60). Eight ascospores were available per tetrad due to an additional mitosis that follows meiosis II in Z. tritici. A previous study, which sequenced all eight ascospores of two tetrads, confirmed that each tetrad contained four twins as the result of the mitosis (60). Based on the molecular characterization of all ascospores of the 23 tetrads from the previous study (60), we selected four unique ascospores for each tetrad for sequencing. Library preparation and sequencing using an Illumina HiSeq 3000 machine for the 84 ascospores was performed at the Max Planck-Genome-Center, Cologne, Germany. The Illumina read data are available in the Sequence Read Archive under BioProject number PRJNA904559. Please note that two tetrads (A03-4 and A08-1) that were sequenced in an earlier study were also included (60). The Illumina read data for these two tetrads are available at the Sequence Read Archive under BioProject number PRJNA438050. An overview of the included tetrads and ascospores is given in Table S3B. We visually checked the distribution of parental haplotypes within each tetrad. Only 23 tetrads showing a 2:2 segregation of parental haplotypes and crossover events involving the homologous chromosome of two ascospores were included (see example in Fig. S6). One of the original 24 tetrads was excluded from the analysis because it showed a pattern of parental haplotypes on multiple chromosomes that was inconsistent with meiotic recombination and was thought to be not the product of a single meiosis but, rather, a mixture of two or more meiotic events. In addition, the parental strain IPO94269 was sequenced as described above, and its Illumina reads were deposited under the BioProject number PRJNA904744. The Illumina data for the IPO323 parental strains are available under BioProject number PRJNA371572.
Mapping and SNP calling. The reads of 92 ascospores of the 23 tetrads were mapped to the IPO323 reference genome with Bowtie 2 (version 2.3.4.1) (104). This analysis focused only on single nucleotide polymorphisms (SNPs). To obtain a high-quality SNP data set, we performed the SNP calling with two variant callers, GATK (version 4.1.6.0) (105) and SAMtools (version 1.7) (106), and SNPs with a QUAL value of $90 that were called by both variant callers were used for the downstream processing. Variants in regions that contain transposable elements (TEs) were removed from the analysis with BEDTools Intersect (version 2.26.0, option -v) (107) to avoid spurious alignments. From the remaining SNPs, only variants from regions with coverage of .5 in all four ascospores of a tetrad were used for the analysis to avoid false negatives due to the low coverage in one of the spores. Biallelic SNPs with minor allele frequency of .0.9 and with QUAL of .90 that were called by both variant callers were identified by overlapping variant call format (VCF) files from both haplotype callers with BEDTools Intersect and used as a core set of high-fidelity SNPs. VCF files of four spores from the same tetrad were merged with VCFtools (v0.1.15) with the merge option (108) to create a variant file for each tetrad.
High-quality threshold and variant calling accepting only variants identified by two callers will potentially lead to false-negative calls in some of the four spores and hence will affect the segregation ratio obtained for the SNPs. To reduce the risk of false-negative calls, we reintroduced high-fidelity SNPs to any of the other ascospores of a tetrad if there was an indication that it was present but did not satisfy the quality requirements. This means that SNPs which were called in a tetrad but did not meet the quality threshold in some of the ascospores but met the quality threshold in at least one of the ascospores were reintroduced. SNPs on chromosomes that were deleted or absent in the parental strains were removed from the analysis of the respective progeny. Please see Table S3B for an overview of the number of high-fidelity SNPs included for each tetrad.
Identification of recombination events. To detect recombination and gene conversion events in the tetrad progeny, CrossOver.py from the ReCombine package (version 2.1) for tetrad analysis in yeast (64) was modified to fit the genome characteristics of Z. tritici (size and number of chromosomes and the location of the centromeres [58]). Input segregation files for the CrossOver program were generated from merged tetrad VCF files for each tetrad with the custom-made bash script (see Text S1). Each segregation file consisted of seven columns: the first two columns referred to the chromosome and position of the variant, the third column served as a spacer, and the last four columns referred to the presence/ absence of SNPs in four spores. Values of 0 and 1 in the last four columns of a segregation file designated the presence or absence of a variant at a certain position compared to the reference genome. The program initially identifies COs as positions with 2:2 segregation where adjacent markers undergo a reciprocal genotype switch. Gene conversion tracts are then identified as regions of non-2:2 segregation. After the identification of recombination events, all double crossovers separated with a single SNP were filtered out. Gene conversions were filtered for tracts spanning $3 markers. The recombination rate per tetrad (cM/Mb) was calculated with the following formula: The gene conversion rate per tetrad was determined as the proportion of converted markers from the total number of markers identified per tetrad. Furthermore, tract lengths were determined with the midpoint method, i.e., the midpoint between two markers of a different class (e.g., converted versus non-converted) was considered to be the position were the tract started or ended. Tracts spanning TEs were removed for the estimation of tract lengths and recombination rates. Recombination events and gene conversion events detected in this study are listed in the supplemental material (Table S1B and C  and Table S2A to C, respectively). SNPs in TEs were disregarded for the determination of recombination and gene conversion events.
Estimation of meiotic mutation rates. To estimate meiotic mutation rates in Z. tritici, genome-wide SNPs (including SNPs in TEs) satisfying the following criteria were taken into consideration: (i) read depth of .5 in both parental strains and the ascospore progeny and (ii) absent in both parental strains (Table S3A). Before a SNP in the progeny was considered a de novo mutation, both parental sequencing results were manually checked to validate that this SNP was already present but not called in the parental genomes. Only SNPs in the progeny that showed no hints in the parental genomes were included in the subsequent analysis. The per bp mutation rate was calculated as the average number of meiotic mutations per ascus/the reference genome size. To verify in silico-detected meiotic mutations, we performed Sanger sequencing of 20 randomly selected mutations, from which 19 mutations were confirmed (see Text S1).
Detection of duplications. For the detection of duplications in the parental strains, Illumina reads were quality filtered as described above and mapped onto the reference genome using SpeedSeq align followed by structural variation analysis using LUMPY (109) as implemented in the SpeedSeq package (version 0.1.2) (110). The VCF files were filtered using BCFtools (version 1.6) as follows: VCF files were filtered on duplications, genotype (GT, 0/1), quality of .400, and length of ,50,000.
For motif discovery in the vicinity of crossover events, 500 bp upstream and downstream genomic sequence surrounding the individual crossover events were extracted from the genome and used for motif discovery with HOMER (version 4.11) (113). The four most significant de novo-identified motifs with a P value of , 1 Â 10 215 were selected.
Data availability. The Illumina read data are available in the Sequence Read Archive under the BioProjects PRJNA904559, PRJNA438050, PRJNA904744, and PRJNA371572. The Z. tritici IPO323 reference genome is available under accession number GCA_000219625.1.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TEXT S1, PDF file, 0.1 MB. FIG S1, TIF 1-1). E.H.S. is also grateful for support from CIFAR. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.