Fanny E Hartmann
email: fanny.hartmann@u-psud.fr

Ricardo C Rodríguez De La Vega

Jean-Tristan Brandenburg

Fantin Carpentier

Tatiana Giraud

Gene presence-absence polymorphism in castrating anther-smut fungi: recent gene gains and phylogeographic structure

Keywords: number variation, pathogen, population structure, adaptive variation

Gene presence-absence polymorphisms segregating within species are a significant source of genetic variation but have been little investigated to date in natural populations. In plant pathogens, the gain or loss of genes encoding proteins interacting directly with the host, such as secreted proteins, probably plays an important role in coevolution and local adaptation. We investigated gene presence-absence polymorphism in populations of two closely related species of castrating anther-smut fungi, Microbotryum lychnidis-dioicae (MvSl) and M. silenes-dioicae (MvSd), from across Europe, on the basis of Illumina genome sequencing data and high-quality genome references. We observed presence-absence polymorphism for 186 autosomal genes (2% of all genes) in MvSl, and only 51 autosomal genes in MvSd. Distinct genes displayed presenceabsence polymorphism in the two species. Genes displaying presence-absence polymorphism were frequently located in subtelomeric and centromeric regions and close to repetitive elements, and comparison with outgroups indicated that most were present in a single species, being recently acquired through duplications in multiple-gene families. Gene presence-absence polymorphism in MvSl showed a phylogeographic structure corresponding to clusters detected based on SNPs. In addition, gene absence alleles were rare within species and skewed toward low-frequency variants. These findings are consistent with a deleterious or neutral effect for most gene presence-absence polymorphism. Some of the observed gene loss and gain events may however be adaptive, as suggested by the putative functions of the corresponding encoded proteins (e.g., secreted proteins) or their localization within previously identified selective sweeps. The adaptive roles in plant and anther-smut fungi interactions of candidate genes however need to be experimentally tested in future studies.

Introduction

Gene presence-absence polymorphisms segregating within species are probably a widespread and important source of genetic variation [START_REF] Conrad | A high-resolution survey of deletion polymorphism in the human genome[END_REF][START_REF] Conrad | Origins and functional impact of copy number variation in the human genome[END_REF]. Several gene copy number changes have been shown to be adaptive, due to their effects on gene expression and phenotypic diversity [START_REF] Henrichsen | Copy number variants, diseases and gene expression[END_REF][START_REF] Orozco | Copy number variation influences gene expression and metabolic traits in mice[END_REF]). However, the proportion and type of genes affected by presence-absence polymorphism in natural populations have been little explored, and such studies are only now being facilitated by new sequencing technologies [START_REF] Schrider | Gene copy-number polymorphism in nature[END_REF][START_REF] Albalat | Evolution by gene loss[END_REF].

Several fungal plant pathogen species have small, compact genomes [START_REF] Möller | Evolution and genome architecture in fungal plant pathogens[END_REF], rendering them highly suitable for the fine mapping of intraspecific differences in gene content, as low repeat content and high-quality genome assemblies facilitate structural variant discovery [START_REF] Guan | Structural variation detection using next-generation sequencing data: A comparative technical review[END_REF]. Furthermore, some of the best-documented examples of adaptive evolution through gene loss or gene gain have been identified in fungal pathogen-host plant interactions [START_REF] Fouché | The birth and death of effectors in rapidly evolving filamentous pathogen genomes[END_REF]. During host infection, fungal pathogens secrete proteins (i.e. effectors) that may be recognized by the host in a gene-for-gene relationship and trigger plant defense responses [START_REF] Jones | The plant immune system[END_REF]. Losses or gains of genes encoding proteins interacting with the host defense system have a strong impact on fitness [START_REF] Presti | Fungal effectors and plant susceptibility[END_REF]. In agricultural ecosystems, the loss of entire effector genes, leading to the emergence of virulence, has been documented in several crop pathogens [START_REF] Orbach | A telomeric avirulence gene determines efficacy for the rice blast resistance gene Pi-ta[END_REF][START_REF] Schürch | Molecular population genetic analysis differentiates two virulence mechanisms of the fungal avirulence gene NIP1[END_REF][START_REF] Gout | Lost in the middle of nowhere: the AvrLm1 avirulence gene of the Dothideomycete Leptosphaeria maculans[END_REF][START_REF] Gout | Genome structure impacts molecular evolution at the AvrLm1 avirulence locus of the plant pathogen Leptosphaeria maculans[END_REF]Hartmann et al. 2017). The acquisition of new effector genes through horizontal gene transfer or duplication events is also associated with the emergence of virulence [START_REF] Friesen | Emergence of a new disease as a result of interspecific virulence gene transfer[END_REF][START_REF] Khang | Genome organization and evolution of the AVR-Pita avirulence gene family in the Magnaporthe grisea species complex[END_REF][START_REF] Jonge R De | Tomato immune receptor Ve1 recognizes effector of multiple fungal pathogens uncovered by genome and RNA sequencing[END_REF]. The widespread use of new resistant host genotypes carrying specific resistance genes probably imposes a strong selective pressure on fungal pathogen populations in agricultural crops, promoting loss of the cognate effector gene or the retention of recently gained genes with potentially beneficial effects in terms of virulence, according to the 'arms race' model of evolution [START_REF] Brown | Plant-parasite coevolution: bridging the gap between genetics and ecology[END_REF]. Comparisons of multiple genomes within agricultural pathogen species have shown that gene content variation can affect hundreds of genes, including a high proportion of effector genes and genes located in rapidly evolving genomic regions, in particular regions rich in repetitive elements [START_REF] Gout | Lost in the middle of nowhere: the AvrLm1 avirulence gene of the Dothideomycete Leptosphaeria maculans[END_REF][START_REF] Yoshida | Association genetics reveals three novel avirulence genes from the rice blast fungal pathogen Magnaporthe oryzae[END_REF][START_REF] Yoshida | Host specialization of the blast fungus Magnaporthe oryzae is associated with dynamic gain and loss of genes linked to transposable elements[END_REF][START_REF] Plissonneau | The evolution of orphan regions in genomes of a fungal pathogen of wheat[END_REF]Hartmann and Croll 2017;[START_REF] Plissonneau | Pangenome analyses of the wheat pathogen Zymoseptoria tritici reveal the structural basis of a highly plastic eukaryotic genome[END_REF].

Despite its importance for co-evolution, we know very little about the extent of gene presenceabsence polymorphism in fungal plant pathogens in natural ecosystems. Several factors may drive specific dynamics of gene content change in natural populations. For example, wild plant parasites undergo frequent population extinction-recolonization events, with major consequences for infection success and population genetic structure [START_REF] Tack | Ecological and evolutionary implications of spatial heterogeneity during the off-season for a wild plant pathogen[END_REF]. High environmental heterogeneity and seed or spore banks also affect pathogen population structure and the process of co-evolution in natural conditions [START_REF] Laine | Heterogeneous selection promotes maintenance of polymorphism in host-parasite interactions[END_REF][START_REF] Koopmann | Fisher-Wright model with deterministic seed bank and selection[END_REF]. Host ecotype and environmental factors also have a strong impact on local pathogen populations, driving local adaptation [START_REF] Laine | Host ecotype generates evolutionary and epidemiological divergence across a pathogen metapopulation[END_REF][START_REF] Stam | The wild tomato species Solanum chilense shows variation in pathogen resistance between geographically distinct populations[END_REF]. In natural conditions, the co-evolutionary dynamics of plant-parasite interactions appear to follow essentially a 'trench warfare' model, favoring balanced polymorphism at loci involved in host-parasite co-evolution [START_REF] Brown | Plant-parasite coevolution: bridging the gap between genetics and ecology[END_REF]. In addition, the virulence of pathogen populations in wild pathosystems is usually a quantitative trait [START_REF] Alexander | Spread of anther-smut disease (Ustilago violacea) and character correlations in a genetically variable experimental population of Silene alba[END_REF][START_REF] Laine | Spatial variation in disease resistance: from molecules to metapopulations[END_REF][START_REF] Thrall | Variation in resistance and virulence among demes of a plant host-pathogen metapopulation[END_REF][START_REF] Koskela | Resistance and tolerance in a host plantholoparasitic plant interaction: genetic variation and costs[END_REF], suggesting that presence-absence variations of effector 6 gene may be a less frequent mechanism of co-evolution than in crops, in which virulence is more often a binary trait, with resistant varieties completely impeding infection [START_REF] Möller | Evolution and genome architecture in fungal plant pathogens[END_REF].

Microbotryum lychnidis-dioicae (MvSl) and M. silenes-dioicae (MvSd) are two closely related species of castrating anther-smut fungi. These highly specialized plant pathogens parasitize the white campion Silene latifolia and the red campion S. dioica, respectively (Le [START_REF] Gac | Phylogenetic evidence of host-specific cryptic species in the anther smut fungus[END_REF][START_REF] Refrégier | Cophylogeny of the anther smut fungi and their caryophyllaceous hosts: prevalence of host shifts and importance of delimiting parasite species for inferring cospeciation[END_REF]. Virulence is mostly a quantitative trait, and no gene-for-gene relationship has been identified in anther-smut fungi [START_REF] Alexander | Genotypic variation in plant disease resistance-physiological resistance in relation to field disease transmission[END_REF][START_REF] Alexander | Spread of anther-smut disease (Ustilago violacea) and character correlations in a genetically variable experimental population of Silene alba[END_REF][START_REF] Biere | Sex-specific costs of resistance to the fungal pathogen Ustilago violacea (Microbotryum violaceum) in Silene Alba[END_REF][START_REF] Chung | Variation in resistance to multiple pathogen species: anther smuts of Silene uniflora[END_REF]). Populations of MvSl are genetically more diverse than those of MvSd, and form distinct genetic clusters, corresponding to the footprints of southern glacial refugia. By contrast, MvSd displays little population differentiation across Europe [START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]. The strong co-structure between the anthersmut pathogen MvSl and its host, together with cross-inoculation experiments, are consistent with plant local adaptation [START_REF] Kaltz | Local maladaptation in the anther-smut fungus Microbotryum violaceum to its host plant Silene latifolia: evidence from a cross-inoculation experiment[END_REF][START_REF] Feurtey | Strong phylogeographic co-structure between the anther-smut fungus and its white campion host[END_REF]. MvSl has undergone numerous selective sweeps, consistent with rapid co-evolution, affecting almost 17% of the genome, whereas only 2% of the genome is affected in MvSd [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]). These two fungal species have total haploid genome size of about 30 Mb; hyphae in plants are dikaryotic and teliospores are diploid. Both species have dimorphic mating-type chromosomes containing large regions without recombination, but with extensive rearrangements and gene losses due to permanent heterozygosity [START_REF] Fontanillas | Degeneration of the non-recombining regions in the mating-type chromosomes of the anther-smut fungi[END_REF][START_REF] Badouin | Chaos of rearrangements in the mating-type chromosomes of the anthersmut fungus Microbotryum lychnidis-dioicae[END_REF][START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]. By contrast, the autosomes display low levels of heterozygosity, due to high rates of selfing and the rarity of outcrossing events [START_REF] Giraud | Mating system of the anther smut fungus Microbotryum violaceum: selfing under heterothallism[END_REF][START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF].

We investigated homozygous gene presence-absence polymorphisms on autosomes of the two sister anther-smut species MvSl and MvSd, using available genome sequence data and highquality reference genome assemblies. The genes displaying presence-absence polymorphism were mostly recently acquired, in a single species, through duplications in multiple-gene families and their absence alleles segregated at low frequencies. The phylogeographic structure of presence-absence polymorphism in MvSl corresponded to the previously identified genetic clusters based on SNPs. Altogether, these findings suggest that most gene presence-absence polymorphism is neutral. Nevertheless, the putative functions of some genes affected by presence-absence polymorphism (e.g., secreted proteins) or their localization within previously identified selective sweeps suggest that some gene loss or gain events may be adaptive. However, functional validation in future studies will be needed to test adaptive roles in plant and anthersmut fungi interactions.

Materials and Methods

Genome data used for gene presence-absence polymorphism calling

We analyzed the genome sequences of 39 MvSl isolates and 19 MvSd isolates collected from across Western Europe, previously obtained with Illumina paired-end sequencing technology, with a mean coverage of 100x [START_REF] Whittle | Recent and massive expansion of the mating-type-specific region in the smut fungus Microbotryum[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]. In total, 14 strains were sequenced in these previous studies as haploids of the a 1 or a 2 mating type, whereas the remaining 47 strains were sequenced as diploids. We downloaded raw data publicly available from the NCBI Short Read Archive (SRA) under the BioProject IDs PRJNA295022 and PRJNA269361 (See summary in Supplementary Table S1). We studied gene presence-absence polymorphism in MvSl and MvSd, using the high-quality reference genomes of the MvSl-1064 strain and the MvSd-1303 strain, previously obtained and annotated for gene models (see summary in Supplementary Table S2, [START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]). We focused on autosomes, which are highly homozygous in these species [START_REF] Badouin | Chaos of rearrangements in the mating-type chromosomes of the anthersmut fungus Microbotryum lychnidis-dioicae[END_REF]. We therefore used only one of the two available haploid genomes for each strain for further analyses (the haploid genomes of the a 1 mating type). We used the MvSl-1064-A1-R4 assembly from the European Nucleotide Archive (ENA), accession number ERS1013679, and the MvSd-1303-D assembly, accession number ERS1436592 [START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF].

Read mapping

We mapped Illumina reads against the reference genomes of each species. We trimmed Illumina raw reads of the MvSl and MvSd genomes for sequence quality and removed adapter sequences with the Cutadapt v1.12 software [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. We used the options: -q 10,10; --minimumlength 50; -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC; -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATT.

We aligned MvSl trimmed reads against the MvSl-1064-A1-R4 genome and MvSd trimmed reads against the MvSd-1303-D a 1 genome. We used the short-read aligner bowtie2 v2.1.0 [START_REF] Langmead | Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[END_REF] for the mapping of trimmed reads with the following software options: - reference genomes, respectively. As differences in coverage between samples may lead to a bias in the specificity and sensitivity of copy number variation detection [START_REF] Guan | Structural variation detection using next-generation sequencing data: A comparative technical review[END_REF], we normalized the coverage of all isolates to a single value (30X), by sampling reads at random from alignment files with the samtools v1.3.1 software [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF]).

Gene presence-absence polymorphism calling

We used a combination of read depth-based and split read-based methods to call homozygous missing fragments in the genomes analyzed, relative to the reference genomes of the MvSl-1064

and MvSd-1303 strains. Mating-type chromosomes contain an extensive region without recombination that is permanently sheltered in a heterozygous state, which led to genomic decay and in particular the loss of hundreds of genes [START_REF] Fontanillas | Degeneration of the non-recombining regions in the mating-type chromosomes of the anther-smut fungi[END_REF][START_REF] Badouin | Chaos of rearrangements in the mating-type chromosomes of the anthersmut fungus Microbotryum lychnidis-dioicae[END_REF]. We therefore focused here exclusively on autosomes. We first used the algorithm implemented in CNVnator v0.3, which uses all mapped reads and performs a statistical analysis of read coverage in bins along the genome sequence [START_REF] Abyzov | CNVnator: an approach to discover, genotype, and characterize typical and atypical CNVs from family and population genome sequencing[END_REF]. We set the bin size at 100 bp and retained only deletion calls, i.e missing fragments, fulfilling the following criteria: p-value < 0.05, length > 500 bp, q0 < 0.4 and normalized read coverage < 0.3. As a second method, we used the Pindel v0.5.7 software, which uses information about paired reads for which only one end can be mapped. Based on the anchor point of the mapped read, the insert size and the direction of the unmapped read, Pindel predicts missing fragments by breaking the unmapped read into fragments and mapping them separately [START_REF] Ye | Pindel: a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads[END_REF]). We converted .bam alignment files into Pindel input format with the sam2pindel tool of Pindel. We ran Pindel with default options and used the Pindel pindel2vcf tool to convert output files into variant calling format. We retained only homozygous deletion calls of more than 500 bp, supported by at least 15 reads (allele depth >15). The calling of missing fragments is more challenging in short contigs, due to frequent poor-quality assembly and high repeat content. We therefore restricted our analyses to the largest contigs (minimum size of 147 kb), which summed covered 90% of the length of the autosomal genome of the MvSl-1064 and MvSd-1303 strains. To detect genes present in the reference genomes that were absent from the Illumina genomes, we retained the missing fragments called with each method that covered at least one gene, over >90% of its length. We assessed the overlaps between missing fragments and gene models masked for repeats, with the bedtools "intersect" command [START_REF] Quinlan | BEDTools: a flexible suite of utilities for comparing genomic features[END_REF]. We considered gene absence events identified by both the read depth-based method and the split read-based method as single gene absence events. We determined the total length of the genomic region affected by missing fragments for each strain by summing the lengths of missing fragments detected by the two methods. To avoid counting twice the lengths of overlapping missing fragments detected by the two methods, we took into account the reunion of overlapping missing fragments.

Quality control of detected gene presence-absence polymorphism

To assess the overall quality of gene presence-absence polymorphism calls, we performed several quality control steps. First, we checked the gene absence events called with our pipeline using Illumina sequencing of the very same strain (Illumina-resequenced MvSl-1064 strain, [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF])) against the reference high-quality genome of the MvSl-1064 strain (Branco et al. 11 2017). Second, we used a comparison of high-quality genome assemblies of two MvSl strains, MvSl-1064 (used as a reference for gene presence-absence calling) and a newly sequenced strain, MvSl-1318 (see section "Analyses of strain-specific genes using de novo genome assemblies" below). We compared the number of gene absence events called with our pipeline by using the Illumina-resequenced MvSl-1064 strain [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]) against the MvSl-1318 reference genome with the number of genes found to be specific to the MvSl-1318 reference genome (i.e., lacking from the MvSl-1064 reference genome) in the comparison of the two high-quality genome assemblies. We studied global synteny between genome assemblies using the nucmer command from MUMmer package v3.1 [START_REF] Kurtz | Versatile and open software for comparing large genomes[END_REF]). To identify gene sequences of the MvSl-1318 strain that were absent in the MvSl-1064 strain, we performed blastn analyses in the genome of the MvSl-1064 strain using the sequences of predicted genes of the MvSl-1318 strain as a query. Genes were considered as absent if no blast hit was found on the orthologous contig with a minimum identity of 90%, a bit score value of 100 and a length of at least 90% of the gene length. We used the Integrative Genomics Viewer (IGV) browser to visually inspect read mapping in genomic regions with detected gene absence events [START_REF] Thorvaldsdóttir | Integrative Genomics Viewer (IGV): highperformance genomics data visualization and exploration[END_REF][START_REF] Robinson | Integrative genomics viewer[END_REF]. Finally, we validated the presence-absence polymorphism inferred in silico for two genes using polymerase chain reactions (PCR) on DNAs from the same strains whose genomes were studied. Primers were designed in the gene coding sequences using Primer 3.0 [START_REF] Rozen | Primer3 on the WWW for general users and for biologist programmers[END_REF]. For the gene MvSl-1064-A1-R4_MC01-1g01717 in the MvSl-1064 reference genome, we performed PCR in nine MvSl strains using 5' GATGTCGATGCGCTCTTTGT 3' and 3' CGTCATCAGTGTGCCCTTTT 5' as forward and reverse primers, respectively.

For the gene

MvSdioicae_1303_FR02_D_N206_PbcR_C014g07422 in the MvSd-1303 reference genome, we performed PCR in nine MvSd strains using 5' GACATCAGGCACCACTCACA 3' and 3' ATCCACCCGTCAAATTCGCA 5' as forward and reverse primers, respectively. PCR reactions were conducted in a 30 µl volume containing 5-10 ng genomic DNA, 0.4mM each of forward and backward primers, 0.25mM dNTP, 0.75U Taq polymerase (DreamTaq, Thermo Fisher, Inc.), PCR buffer. PCR products were amplified for 30 cycles. The resulting amplicons were examined on 1.5% agarose gels. We also tested for significant differences in gene absence calls between haploid and diploid strains within MvSl and MvSd populations.

Gene ontology analysis and functional annotation

We assigned genes to functional categories using InterproScan v5.24-63.0 (Zdobnov and Apweiler 2001) that provided information for protein family (Pfam), gene ontology (GO) and pathways. We used SignalP v4.1 [START_REF] Petersen | SignalP 4.0: discriminating signal peptides from transmembrane regions[END_REF]) and TMHMM v2.0 [START_REF] Krogh | Predicting transmembrane protein topology with a hidden markov model: application to complete genomes[END_REF]) to predict putative secreted proteins. We defined secreted proteins as proteins with a predicted signal peptide but no predicted transmembrane helices. Putative small secreted proteins were defined as secreted proteins with a size inferior to 300 amino acids and unknown functional domain. We filtered out from the set of genes affected by presence-absence polymorphism all genes with a transposable element activity related domain, i.e. reverse transcriptase (PF00078, PF07727), integrase (PF00665, PF13976), Helitron helicase-like domain (PF14214) and GAGpolypeptide of long terminal repeat (LTR) copia-type (PF14223). GO enrichment analysis of genes affected by presence-absence polymorphism was performed with the hypergeometric test implemented in the R package GOstats. GO enrichment p-values were adjusted for multiple testing with the Bonferroni correction and [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] correction using the p.adjust() R function.

Annotation of repetitive elements

Annotation of repetitive elements was performed using a set of repeat sequences (available upon request) detected de novo in a large set of 41 genomes of the Microbotryum genus. Briefly, repeat motifs were detected using a combination of two tools, LTRharvest [START_REF] Ellinghaus | LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons[END_REF]) that performs high-quality annotation of LTR retrotransposons and RepeatModeler [START_REF] Smit | RepeatModeler Open-1.0[END_REF]) that uses results from three other programs, RECON, RepeatScout, and Tandem Repeats Finder. Detected repeat motifs were clustered based on sequence similarities (min. 90% of similarity) and the best representative sequences per cluster, i.e. the centroid sequences, were kept. A set of 3599 repeat sequences was built after removing redundant sequences among centroid sequences, and keeping only transposable element sequences present in at least seven genomes. We annotated genome assemblies using the set of 3599 repeat sequences as library in RepeatMasker v4.0.3 [START_REF] Smit | RepeatMasker Open-4.0[END_REF]. We masked gene models for transposable elements in all the genomes. Telomeric and centromeric repeats were identified by performing a blastn search of telomeric and centromeric motifs previously identified in the MvSl-1064 genome [START_REF] Badouin | Chaos of rearrangements in the mating-type chromosomes of the anthersmut fungus Microbotryum lychnidis-dioicae[END_REF]. Subtelomeres were defined as regions within 100 kb of the contig extremities containing telomeric repeats. Centromeres were defined as regions within 100 kb of the contig extremities containing centromeric repeats. To assess expression levels of genes affected by presence-absence polymorphisms, we used RNAseq data collected using the MvSl-1064 strain (i.e. the same strain as the MvSl reference genome) both using in vitro cultures and in vivo late infection [START_REF] Perlin | Sex and parasites: genomic and transcriptomic analysis of Microbotryum lychnidis-dioicae, the biotrophic and plant-castrating anther smut fungus[END_REF]. We trimmed raw reads for using cutadapt v1.12 [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]) and mapped trimmed reads to the MvSl-1064 genome using TopHat v2.1.1 [START_REF] Trapnell | TopHat: discovering splice junctions with RNA-Seq[END_REF]). We used the featureCounts program from the Subread v1.5.3 package [START_REF] Liao | featureCounts: an efficient general purpose program for assigning sequence reads to genomic features[END_REF] to calculate absolute transcript abundance (i.e. number of mapped reads for each transcript). For each tested condition, we computed reads per kilobase of transcript per million mapped reads (RPKM) values using a customized R script.

Expression data

Comparative genomics analyses

For comparative genomics analyses, we used three additional genomes from previously sequenced Microbotryum species [START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]): M. violaceum s. str. collected on S. nutans (MvSn ERS1013671), M. lagerheimii collected on S. vulgaris (MvSv ERS1013677), and M.

intermedium collected on Salvia pratensis (Mint ERS1324257). To investigate gene orthology relationships between species, we compared the predicted protein sequences of the genome assemblies of MvSl-1064, MvSd-1303, MvSn, MvSv and Mint, with blastp+ v2.30. The output was processed in orthAgogue [START_REF] Ekseth | orthAgogue: an agile tool for the rapid prediction of orthology relations[END_REF]) that uses Markov clustering to obtain orthologous groups. To study the synteny between the MvSl-1064 and MvSd-1303 genomes, we used the nucmer command from MUMmer package v3.1 [START_REF] Kurtz | Versatile and open software for comparing large genomes[END_REF]) and performed blastn+ v2.30 analysis on repeat masked genomic sequences. Blast hits were stringently filtered for a minimum alignment length of 500 bp, a minimum identity of 80%, and a bit score value of 80.

Population structure and statistical analyses

Principal component analyses were computed in the R package {ADEGENET} (dudi.pca function ;[START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF][START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF]. Neighbor-joining trees were performed with the R package {Ape} [START_REF] Paradis | APE: Analyses of Phylogenetics and Evolution in R language[END_REF]). These analyses were performed based on unique events of gene presence-absence polymorphism, i.e. using only a single gene per missing fragments, to be conservative regarding the number of gene gain or loss events. To compare the population structure based on gene presence-absence polymorphism with the population structure based on SNPs, we called genome-wide SNPs for MvSl strains against the MvSl-1064-A1-R4

genome and for MvSd strains against the MvSd-1303-D a 1 genome. For SNP calling, we used GATK version 3.7 [START_REF] Mckenna | The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data[END_REF]. We ran HaplotypeCaller on each strain individually using a diploid mode. We performed joint variant calls using GenotypeGVCFs on a merged gvcf variant file. SNP call filtering for quality was performed using VariantFiltration and following the GATK Good Practices for hard-filtering of variants (QUAL < 250; QD < 2; MQ < 30.0; -10.5 > MQRankSum > 10.5; -5 > ReadPosRankSum > 5; FS > 60;SOR > 3). We kept only bi-allelic SNPs with a high genotyping rate (> 90%) on autosomes. In total, we detected 216,878 bi-allelic SNPs in MvSl and 33,694 bi-allelic SNPs in MvSd. The variant call format files are available upon request. Pairwise F ST were calculated for SNPs and gene presence-absence polymorphisms between genetic clusters using the R package {hierfstats} [START_REF] Goudet | hierfstat, a package for r to compute and test hierarchical F-statistics[END_REF] implementing Yang's algorithm [START_REF] Yang | Estimating hierarchical F-Statistics[END_REF]. To test the significance of the F ST value obtained between the MvSd populations that we found differentiated based on gene presence-absence polymorphism, we performed a permutation test with 1000 bootstraps, assigning randomly cluster identity to each strain, keeping the cluster sizes constant, and computing mean F ST values for each permutation. We performed Kruskal-Wallis rank sum test and Pearson's chi-squared test in the R software v3.3.3. Plots were performed using the software Circos v0.67-7 [START_REF] Krzywinski | Circos: an information aesthetic for comparative genomics[END_REF]) and the R package {ggplot2} [START_REF] Wickham | ggplot2: elegant graphics for data analysis[END_REF].

Analyses of strain-specific genes using de novo genome assemblies

For identifying genes absent from the reference genome but present in the genomes of other MvSl strains, we performed de novo genome assemblies. As MvSl showed a population structure with four genetic clusters [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]), we performed de novo assemblies for three MvSl strains belonging to the genetic clusters different from that of the reference MvSl-1064 genome.

We ). We only used the a 1 haploid genome for our analyses. Gene models were masked for repeats.

To detect strain-specific genes in the assembled genomes of the MvSl-1318, MvSl-140-01 and MvSl-1005 strains, we applied the same procedure as the one used to detect genes present in the reference genomes and absent from the Illumina genomes, for consistency. We mapped Illumina sequencing data of the MvSl-1064 strain [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]) against the three genome assemblies and called gene presence-absence polymorphism with our gene presence-absence detection pipeline. We restricted our analyses to the largest contigs, which covered 90% of the autosomal genome length in the three genomes (minimum size of 317'619 bp, 11'561 bp and 11'985 bp in the MvSl-1318 strain, MvSl-140-01 strain, and MvSl-1005 strain respectively). We identified autosomal contigs in the three genome assemblies by global synteny analyses with the genome assemblies of the MvSl-1064 strain using the nucmer command from MUMmer package [START_REF] Kurtz | Versatile and open software for comparing large genomes[END_REF]). As the MvSl-140-01 and MvSl-1005 strains were sequnced as diploid, we excluded any contigs with large synteny with the mating-type chromosomes a 1 and a 2 . The MvSl-1064-A2-R4 assembly (the haploid genomes of a 2 mating type) was accessed from ENA at accession number ERS459551. To analyze population structure of MvSl-1318-specific genes, we called gene presence-absence polymorphism with our pipeline in all MvSl strains against the MvSl-1318 assembly and selected missing fragments affecting the MvSl-1318-specific genes. Mean genome-wide mapping coverage of MvSl strains ranged from 53X to 345X on the MvSl-1318-D reference genome.

Results

Gene presence-absence event calling and quality control

We called homozygous gene presence-absence polymorphisms by comparing the available high- S2C). In MvSl, 473 missing fragments affected a single gene and 38

affected groups of at least five genes. In MvSd, 48 missing fragments affected a single gene and one affected a group of five genes. In total, we detected 953 gene absence events in MvSl, and 161 in MvSd. The number of missing genes detected per strain ranged from two to 50, with a median value of 24 genes in MvSl, and one to 15 with a median value of eight genes in MvSd (Supplementary figure S2D).

Several quality control checks indicated that the inferences for the missing fragments were reliable. First, no missing fragments were detected following genome resequencing of the MvSl-1064 strain [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF], the source of the reference genome sequence. Second, we found a large overlap between the missing fragments detected with our pipelines and those detected by comparing high-quality genome assemblies for two MvSl strains. Out of 15 genes detected to be missing in one strain with our pipeline, four genes were detected present using blast analyses on de novo assemblies (Supplementary Text S1; Supplementary figure S3); however, blast analyses can match paralogs. The visual inspection of read mapping in genomic regions with detected gene absence events using a genome browser suggested a rate of false positive of about 15% for all analyzed genomes. Finally, PCR tests validated in silico inferences for two genes affected by presence-absence polymorphism. The gene MvSl-1064-A1-R4_MC01-1g01717 couldn't not be amplified in five MvSl strains that were predicted to lack the gene in silico. The gene MvSdioicae_1303_FR02_D_N206_PbcR_C014g07422 couldn't not be amplified in three MvSd strains that were predicted to lack the gene in silico.

Contrasting gene presence-absence polymorphism patterns between species

We found that different genes and different numbers of genes displayed presence-absence polymorphism in MvSl and MvSd. We identified 186 autosomal genes displaying presenceabsence polymorphism in the set of 38 MvSl strains studied, corresponding to 2% of the total gene content of MvSl-1064 autosomes (Figure 1, right panel; Supplementary Table S3). We identified 51 autosomal genes displaying presence-absence polymorphism in the set of 19 MvSd strains, corresponding to 0.6% of the total gene content of MvSd-1303 autosomes (Figure 1, left panel; Supplementary Table S3). Fewer genes displayed presence-absence polymorphism in MvSd than in MvSl, consistent with the lower SNP and microsatellite diversity previously reported in MvSd [START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]). In agreement with this diversity difference, we identified 216,878 bi-allelic SNPs in MvSl and 33,694 bi-allelic SNPs in MvSd by mapping MvSl and MvSd Illumina genomes on their respective reference genome. Gene absence alleles segregated at low frequency in both species, with a median frequency of 5% in MvSl and 10% in MvSd (Figure 1E; Supplementary Table S3).

We searched for orthology relationships between the genes identified as displaying presenceabsence polymorphism on MvSl-1064 and MvSd-1303 autosomes. The two closely related species shared 7,869 groups of orthologous autosomal genes, containing 83.3% and 78.2% of the total gene content of MvSl-1064 and MvSd-1303 autosomes, respectively. In total, 59% of the genes displaying presence-absence polymorphism in MvSl and 49% of those in MvSd belonged to orthologous groups present in both species (Figure 2B). However, only four orthologous groups were affected by presence-absence polymorphism in both species (Figure 2A). Gene presence-absence polymorphism affected different genomic regions in the two species (Figure 1G).

Functional characteristics, evolutionary history and genomic environment of the genes displaying presence-absence polymorphism

Most of the genes displaying presence-absence polymorphism encoded proteins of unknown function, to a greater extent than genome average. Indeed, a total of 86% S4). We found no enrichment in any specific cellular compartment. In particular, genes encoding secreted proteins, which frequently play a major role in the interactions of fungal pathogens with their hosts [START_REF] Fouché | The birth and death of effectors in rapidly evolving filamentous pathogen genomes[END_REF], were not overrepresented among the genes displaying presence-absence polymorphism. Three genes encoding small secreted proteins were nevertheless absent from some MvSl and MvSd strains. In MvSl, genes encoding a secreted peptidase and a histidine phosphatase, both potentially involved in pathogenicity [START_REF] Monod | Secreted proteases from pathogenic fungi[END_REF][START_REF] Vincent | Editorial: How can secretomics help unravel the secrets of plant-microbe interactions? Front[END_REF][START_REF] Albataineh | Regulatory roles of phosphorylation in model and pathogenic fungi[END_REF], also displayed presence-absence polymorphism (Supplementary Table S3). Five genes encoding proteins involved in RNA interference activity (Piwi domain, Argonaute domain) were absent in one MvSd strain (Supplementary Table S3). Some plant pathogens use the RNA interference machinery to prevent the expression of host immunity genes [START_REF] Weiberg | Small RNAs-the secret agents in the plant-pathogen interactions[END_REF].

We investigated the evolutionary history of genes displaying presence-absence polymorphism by searching for orthologs in three Microbotryum species closely related to the sister species MvSl and MvSd, using the available genomes of M. violaceum s. str. (MvSn), M. lagerheimii (MvSv), and M. intermedium (Mint) [START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]. Most of the genes displaying presence-absence polymorphism (46.7% in MvSl and 54.8% in MvSd, Figure 2B) were recently acquired i.e., were present only in one or both of the sister species MvSl and MvSd and not in the other Microbotryum species. Only 5.9% of the genes displaying presence-absence polymorphism in MvSl and 9.8% of those displaying presence-absence polymorphism in MvSd were present in all five Microbotryum species. This finding indicates that presence-absence polymorphism is rare in ancient genes, that were present in the common ancestor of the five Microbotryum species. Most of the genes displaying presence-absence polymorphism in MvSl and MvSd instead corresponded to species-specific genes, i.e. genes recently gained in a single species, after their divergence from other species (Figure 2B).

The genes displaying presence-absence polymorphism in MvSl were found in subtelomeric regions (defined as regions within 100 kb of the contig extremities and containing telomeric repeats) more frequently than would be expected by chance (Figure 3A). A total of 10.7% of genes in subtelomeric regions were affected by presence-absence polymorphism compared to 2% in the whole genome. We checked by visual inspection that the gene absence events called in subtelomeric regions were not artifacts due to low mapping quality. The small number of subtelomeric repeats in the MvSd-1303 genome precluded testing for an enrichment in gene presence-absence polymorphism in subtelomeric regions in MvSd, but gene absence events were frequently found close to the ends of the contigs (i.e. within 100 kb of the contig extremities).

The genes displaying presence-absence polymorphism in MvSl were also found closer to centromeric regions (defined as regions within 100 kb of the contig extremities and containing centromeric repeats) than would be expected by chance (Figure 3B). The genes displaying presence-absence polymorphism were also significantly closer to transposable elements than other genes. In MvSl, 42% of the genes displaying presence-absence polymorphism were located within 5kb of the nearest transposable element, against only 19% of other genes (Figure 3C). In MvSd, 69% of genes displaying presence-absence polymorphism were located within 5kb of the nearest transposable element, against only 23% of other genes (Figure 3D). We found significant differences in the abundances of the transposable element families closest to the genes displaying presence-absence polymorphism and those closest to other genes in MvSl, but not in MvSd (Figure 3E-F). LTR transposons, the most frequent family of repeats in the MvSl-1064 reference genome [START_REF] Hood | Repetitive DNA in the automictic fungus Microbotryum violaceum[END_REF][START_REF] Perlin | Sex and parasites: genomic and transcriptomic analysis of Microbotryum lychnidis-dioicae, the biotrophic and plant-castrating anther smut fungus[END_REF], were significantly more frequently found closest to genes displaying presence-absence polymorphism than to other genes (one-tailed Fisher's exact test p-value = 0.0007345 in MvSl; not significant in MvSd). Conversely, members of the Helitron family were less frequently found closest to genes displaying presence-absence polymorphism than to other genes in MvSl (one-tailed Fisher's exact test p-value = 0.0001393; not significant in MvSd).

Population structure of gene presence-absence events

We investigated the population structure of gene presence-absence polymorphism in each species. For MvSl, we identified four genetic clusters based on the three first components of a principal component analysis and a dendrogram (Figure 4A-D). These clusters were consistent with the population subdivision previously described based on genome-wide SNPs (Figure 4C)

and corresponding to the phylogeographic footprints of glacial refugia [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]). The congruence of the genetic subdivisions based on SNPs and gene presence-absence polymorphism in MvSl resulted from the large proportion of cluster-specific presence or absence among the genes affected by presence-absence polymorphisms (64%; Figure 5A). Within MvSl, we found higher mean pairwise F ST for SNPs (mean value = 0.32 for pairwise comparisons) than for gene presence-absence polymorphisms (mean value = 0.15). Within MvSd, we observed low levels of population structure (Figure 4E-H), consistent with previous findings for SNPs (Figure 4G) and microsatellites [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF][START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF]). However, a principal component analysis of gene presence-absence polymorphism suggested here the existence of genetic differentiation between strains across one east-west longitudinal gradient (Figure 4E-F). When calculating mean pairwise F ST between the eight MvSd strains sampled in Eastern Europe and the eleven MvSd strains sampled in Western Europe (Figure 4E), we found higher differentiation for gene presence-absence (mean pairwise F ST value = 0.12) than for SNPs (mean pairwise F ST value = 0.08). The observed mean F ST value calculated between these two groups for gene presenceabsence polymorphisms was higher than expected by chance, lying within the 1% extreme values in 1000 permutations.

We investigated the possible presence of gene presence-absence polymorphisms in regions previously identified to have been affected by a recent selective sweep in the MvSl Northwestern genetic cluster [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF] MvSl Northwestern strain were located within a selective sweep region, including three genes encoding secreted proteins (Figure 5B-C; Supplementary Table S3). In selective sweep regions, seven gene absence were rare (frequency < 5%) and three gene absence events were almost fixed (frequency > 70%), as expected for selective sweeps (Supplementary figure S4). Some of the gene presence-absence events, particularly those concerning genes encoding secreted proteins, may thus be adaptive, although functional studies are needed to confirm this.

For identifying genes absent from the MvSl-1064 genome reference but present in the genomes of other MvSl strains, we generated a high-quality genome assembly for the MvSl-1318 strain from the MvSl Eastern cluster and we assembled de novo the Illumina genomes of the MvSl-140-01 strain from the MvSl Northwestern cluster, and of the MvSl-1005 strain from the MvSl Southwestern cluster (Figure 4A; Supplementary Text S2). Using our gene presence-absence detection pipeline, we showed that the MvSl-1318 genome contained 11 genes absent from the MvSl-1064 genome (Figure 6A). We identified nine genes in the MvSl-140-01 genome and 15 in the MvSl-1005 genome that were absent from the MvSl-1064 genome. The number of strainspecific genes identified in the three de novo assembled genomes (i.e. absent from the MvSl-1064 genome) was consistent with the number of genes present in the MvSl-1064 genome and absent from other MvSl strains (Supplementary figure S2D). The population structure of the presenceabsence polymorphism of strain-specific genes in MvSl confirmed the MvSl phylogeographic structure corresponding to the plant local adaptation (Figure 6C-D).

Discussion

Gene content variation affects recently gained genes in anther-smut fungi

Using a robust gene presence-absence calling procedure based on two different read mapping methods, we characterized the degree of gene content variation on autosomes of the fungal anther-smut species MvSl and MvSd. Autosomal chromosomes are highly homozygous and have a total assembly size of about 30 Mb in both species [START_REF] Badouin | Chaos of rearrangements in the mating-type chromosomes of the anthersmut fungus Microbotryum lychnidis-dioicae[END_REF][START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]).

The read depth-based and split read-based methods were complementary, as they were based on different read mapping approaches and detected different gene presence-absence polymorphism events. We detected two to 50 missing genes per strain relative to the reference genome for the species concerned. We performed de novo genome assemblies for three additional MvSl strains to validate the degree of intraspecific gene content variation between two strains, as we detected about a dozen genes that were absent from the MvSl reference genome but present in other strains. Our results supported the view that one single reference genome underestimates the gene space of fungal species [START_REF] Plissonneau | Pangenome analyses of the wheat pathogen Zymoseptoria tritici reveal the structural basis of a highly plastic eukaryotic genome[END_REF][START_REF] Peter | Genome evolution across 1,011 Saccharomyces cerevisiae isolates[END_REF]. The level of gene presenceabsence polymorphism in anther-smut pathogens and the sizes of the missing DNA fragments were consistent with findings for yeasts. In Saccharomyces cerevisiae, gene content variation concerns a median of 31 genes per strain [START_REF] Bergström | A high-definition view of functional genetic variation from natural yeast genomes[END_REF]. We found that gene presenceabsence polymorphism contributed to the genetic variation of anther-smut fungi, consistent with the view that gene gains and losses are important sources of fungal genetic variation [START_REF] Taylor | Sources of fungal genetic variation and associating It with phenotypic diversity[END_REF][START_REF] Fouché | The birth and death of effectors in rapidly evolving filamentous pathogen genomes[END_REF].

The functional characteristics of the genes displaying presence-absence polymorphism were similar in MvSl and MvSd: mostly genes of unknown function, with low expression levels, small sizes and enriched in genes from multiple-copy gene families. These functional characteristics suggest that the genes displaying presence-absence polymorphism do not have housekeeping functions and have been gained recently. Our finding that a large proportion of the genes displaying presence-absence polymorphism belong to multiple-gene families suggests that these genes originate from recent duplication events in multiple-gene families. Following a duplication event, a gene is gained and one of the copies can diverge extensively and acquires a new function [START_REF] Tautz | The evolutionary origin of orphan genes[END_REF]. Short length and low level of expression are also common features of genes that have recently arisen de novo, as reported for Zymoseptoria tritici, a fungal pathogen of wheat [START_REF] Plissonneau | The evolution of orphan regions in genomes of a fungal pathogen of wheat[END_REF]Hartmann and Croll 2017). A search for orthologs in outgroup species also suggested that most of the presence-absence polymorphisms observed probably resulted from recent gene gains in MvSl or MvSd. New genes may arise from noncoding DNA through the spontaneous evolution of an open reading frame and the gain of cisregulatory elements [START_REF] Carvunis | Proto-genes and de novo gene birth[END_REF][START_REF] Mclysaght | New genes from non-coding sequence: the role of de novo protein-coding genes in eukaryotic evolutionary innovation[END_REF], or through horizontal gene transfer (Marcet-Houben and Gabaldón 2010; [START_REF] Slot | Horizontal transfer of a large and highly toxic secondary metabolic gene cluster between fungi[END_REF][START_REF] Ropars | Adaptive horizontal gene transfers between multiple cheese-associated fungi[END_REF], We found that gene presence-absence polymorphism was particularly prevalent in subtelomeric regions and close to repeats, as previously reported for yeasts and for the plant pathogens M. oryzae and Z. tritici [START_REF] Steenwyk | Copy number variation contributes to cryptic genetic variation in outbreak lineages of Cryptococcus gattii from the North American Pacific Northwest[END_REF][START_REF] Yoshida | Host specialization of the blast fungus Magnaporthe oryzae is associated with dynamic gain and loss of genes linked to transposable elements[END_REF][START_REF] Plissonneau | The evolution of orphan regions in genomes of a fungal pathogen of wheat[END_REF][START_REF] Gallone | Domestication and divergence of Saccharomyces cerevisiae beer yeasts[END_REF], consistent with the view that duplication and excision events may be driven by repeated elements. The high content of repeats found in subtelomeres and centromeres (Figure 1B-1D) might at least partly explain the enrichment of gene presence-absence polymorphism events in these genomic regions.

Gene presence-absence polymorphism and host-pathogen coevolution in natural environments

Most of the gene presence-absence polymorphism detected is likely neutral. Indeed, previously identified selective sweep regions were not enriched in gene presence-absence events. In addition, gene absence alleles were skewed toward low frequency variants in the two anther-smut fungal species. Neutral evolution for gene presence-absence polymorphism is also consistent with previous findings in other plant fungal pathogens. Whole-genome sequence comparisons between

Magnaporthe oryzae strains specific to rice, millet, wheat or oat revealed for instance that most of the genes displaying presence-absence polymorphism were genes whose functional domains were present multiple times in the genome, i.e. genes likely displaying high levels of functional redundancy [START_REF] Yoshida | Host specialization of the blast fungus Magnaporthe oryzae is associated with dynamic gain and loss of genes linked to transposable elements[END_REF]). In the wheat pathogen Z. tritici, mainly weak negative selection and neutrality have been shown to act on presence-absence polymorphism of recently gained genes, but also divergent selection in some cases (Hartmann and Croll 2017). In the anther-smut pathogens MvSl and MvSd, we also found that a few of the recently gained genes may be adaptive, with predicted functions potentially involved in host-pathogen interactions, such as secreted proteins, and/or being located within regions subject to recent selective sweeps.

The enrichment in functions linked to cellular stress responses of genes displaying presenceabsence polymorphism may also reflect adaptation to specific environments. The loss of genes encoding proteins involved in RNA interference activity might also be adaptive, as some plant pathogens use the RNA interference machinery to prevent the expression of host immunity genes [START_REF] Weiberg | Small RNAs-the secret agents in the plant-pathogen interactions[END_REF], but this has not been investigated in anther-smut fungi. We found that different genomic regions displayed gene presence-absence polymorphism in MvSl and MvSd specialized on different host plants, with differences in the affected secreted protein-encoding genes, in particular. The genes involved in the co-evolution of anther-smut fungi with their host plants are unknown [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF], and the genes displaying presence-absence polymorphism within MvSl and MvSd populations and corresponding to secreted proteins and/or located within selective sweeps represent potentially interesting candidates for future exploration.

Functional validation assays using gene transformation tools recently established for anther-smut fungi [START_REF] Toh | Reliable transformation system for Microbotryum lychnidis dioicae informed by genome and transcriptome project[END_REF]) should be performed in future studies to validate an adaptive role of the gene presence-absence events with relevant functions of genomic regions.

We identified only five genes encoding secreted proteins displaying presence-absence polymorphism in MvSl, and three in MvSd. These numbers correspond to lower percentages of total gene content than reported for several fungal crop pathogens, such as the rice blast pathogen M. oryzae, the wheat pathogens Z. tritici and Stagonospora nodorum, and the generalist crop pathogens Verticillium dahliae and Fusarium oxysporum [START_REF] Yoshida | Association genetics reveals three novel avirulence genes from the rice blast fungal pathogen Magnaporthe oryzae[END_REF][START_REF] Ma | Comparative genomics reveals mobile pathogenicity chromosomes in Fusarium[END_REF][START_REF] Syme | Resequencing and comparative genomics of Stagonospora nodorum: Sectional gene absence and effector discovery[END_REF][START_REF] Jonge R De | Extensive chromosomal reshuffling drives evolution of virulence in an asexual pathogen[END_REF][START_REF] Plissonneau | The evolution of orphan regions in genomes of a fungal pathogen of wheat[END_REF]. Many gene-for-gene relationships have been documented in crop-pathogen systems, several of which corresponded to recent losses or gains of genes in the fungus, to prevent host recognition by the plant defense system [START_REF] Orbach | A telomeric avirulence gene determines efficacy for the rice blast resistance gene Pi-ta[END_REF][START_REF] Gout | Lost in the middle of nowhere: the AvrLm1 avirulence gene of the Dothideomycete Leptosphaeria maculans[END_REF][START_REF] Gout | Genome structure impacts molecular evolution at the AvrLm1 avirulence locus of the plant pathogen Leptosphaeria maculans[END_REF][START_REF] Jonge R De | Tomato immune receptor Ve1 recognizes effector of multiple fungal pathogens uncovered by genome and RNA sequencing[END_REF]Hartmann et al. 2017;[START_REF] Fouché | The birth and death of effectors in rapidly evolving filamentous pathogen genomes[END_REF].

By contrast, no gene-for-gene relationship has been reported for anther-smut fungi, despite extensive studies on genetic variation in host resistance and fungal pathogenicity. Instead, the probability of infection is a quantitative trait [START_REF] Alexander | Genotypic variation in plant disease resistance-physiological resistance in relation to field disease transmission[END_REF][START_REF] Alexander | Spread of anther-smut disease (Ustilago violacea) and character correlations in a genetically variable experimental population of Silene alba[END_REF][START_REF] Biere | Sex-specific costs of resistance to the fungal pathogen Ustilago violacea (Microbotryum violaceum) in Silene Alba[END_REF][START_REF] Chung | Variation in resistance to multiple pathogen species: anther smuts of Silene uniflora[END_REF], which likely corresponds to a control of host recognition by other genetic changes than gene loss or gain. Differences in co-evolutionary dynamics between anthropized and natural pathosystems, such as the 'arms race' and 'trench warfare' evolution models [START_REF] Brown | Plant-parasite coevolution: bridging the gap between genetics and ecology[END_REF] may also result in different mechanisms of adaptation to the host. More studies of this type are required, on non-crop pathogens in particular, to test hypotheses concerning the evolutionary causes of structural variation frequencies affecting genes involved in co-evolutionary processes.

Diversity and phylogeographic structure

Gene presence-absence polymorphisms were more frequent in MvSl than in MvSd (2% vs. 0.5% of all autosomal genes). The degree of gene presence-absence polymorphism in MvSd may be slightly underestimated, due to the lower reference genome quality for MvSd than for MvSl and the lower number of sequenced strains (Supplementary Table S2). However, this is unlikely to account for the magnitude of the difference observed, which was consistent with those reported for SNPs and microsatellites [START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]).

The analysis of gene presence-absence polymorphism in MvSl revealed a population structure highly similar to that reported on the basis of SNPs and microsatellite data, with the same four genetic clusters corresponding to footprints of glacial refugia [START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]). The addition of MvSl strains [START_REF] Whittle | Recent and massive expansion of the mating-type-specific region in the smut fungus Microbotryum[END_REF] compared to previous phylogeographic studies revealed very high levels of genetic diversity among the strains of the Southwestern cluster, as expected for strains derived from southern glacial refugia, which have been little explored to date [START_REF] Vercken | Glacial refugia in pathogens: European genetic structure of anther smut pathogens on Silene latifolia and Silene dioica[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]. The correspondence between the structures revealed by analyses of neutral markers and those revealed by analyses of gene presence-absence polymorphism further reinforce the view that most of such polymorphism is neutral. In particular, gene gain or loss within gene families may be neutral due to functional redundancy [START_REF] Albalat | Evolution by gene loss[END_REF]. However, the host plant S. latifolia also displays the same genetic subdivision in Europe and local adaptation between clusters: plants from a given cluster are more resistant to MvSl strains from the same cluster [START_REF] Feurtey | Strong phylogeographic co-structure between the anther-smut fungus and its white campion host[END_REF]. Some of the gene presence-absence polymorphism events that are congruent with these clusters may, therefore, also correspond to adaptive polymorphism resulting from co-evolution with the host plant. Interestingly, this analysis of gene presence-absence polymorphism suggested the existence of a geographical population structure in MvSd, across an east-west gradient in Europe. Although this population structure had not been detected before, our analysis of genome-wide SNP data called against an MvSd reference genome partly recovered a similar population structure.

However, PCAs and F ST values indicated a stronger structure in gene presence-absence polymorphism than in SNPs in MvSd, which may reflect distinct population genetic properties of genetic markers in populations, such as mutation rate, or may result from an effect of selection on some of the gene gains/losses. The population structure based on gene presence-absence polymorphism in MvSd may thus reflect, at least partly, adaptive events, possibly in response to a population subdivision of the host in terms of resistance genes. Indeed, a population subdivision of the host plant species S. dioica separates the Western and Eastern populations [START_REF] Hathaway | Geographically congruent large-scale patterns of plastid haplotype variation in the European herbs Silene dioica and S. latifolia (Caryophyllaceae)[END_REF][START_REF] Rautenberg | Geographic and phylogenetic patterns in Silene section Melandrium (Caryophyllaceae) as inferred from chloroplast and nuclear DNA sequences[END_REF]). However, here again, functional validation studies are needed to test this hypothesis.

In conclusion, our findings show that gene presence-absence polymorphism contributes to intraspecific genetic variation in anther-smut fungi, although to a lower extent in terms of total gene content than previously reported for fungal crop pathogens. Few genes encoding secreted proteins showed presence-absence polymorphisms, as expected given the lack of gene-for-gene relationships in anther-smut fungi. Gene presence-absence polymorphism mostly affected recently gained genes, found in a single species, which were likely mostly neutral, with a few cases perhaps corresponding to innovations allowing adaptation to the host or interaction with the

Fig. 1 :

 1 Fig. 1: Genome-wide location of gene presence-absence polymorphism events and their

Fig. 2 :

 2 Fig. 2: Evolutionary history of genes displaying presence-absence polymorphism. A. Venn

Fig. 3 :

 3 Fig. 3: Genome features associated with gene presence-absence polymorphism in

Fig. 4 :

 4 Fig. 4: Population structure of gene presence-absence polymorphism in Microbotryum

Fig. 5 :

 5 Fig. 5: Gene presence-absence polymorphism within individual genetic clusters in

Fig. 6 :

 6 Fig. 6: Identification of strain-specific genes in three Microbotryum lychnidis-dioicae (MvSl)

  [START_REF] Mckenna | The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data[END_REF]. Mean genome-wide mapping coverage of MvSl strains ranged from 54X to 350X and 72X to 156X for the MvSl-1064-A1-R4 and MvSd-1303-D

very-sensitive-local -phred33 -X 1000. We removed PCR duplicates with the MarkDuplicates tool of Picard tools version 1.88 (http://broadinstitute.github.io/picard). We improved alignment accuracy in indel regions, by locally realigning the mapped reads with the RealignerTargetCreator and IndelRealigner tools of the Genome Analysis Toolkit (GATK) Downloaded from https://academic.oup.com/gbe/advance-article-abstract/doi/10.1093/gbe/evy089/4990910 by guest on 21 May 2018 version 3.7 (

  to belong to the MvSl Eastern genetic cluster), with P6/C4 Pacific Biosciences SMRT technology. The sequencing and assembly of the MvSl-1318 genomes were generated as described in[START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]. Briefly, we assembled the genome using the wgs-8.2 version of the PBcR assembler and polished the assembly with the quiver software (https://github.com/PacificBiosciences/GenomicConsensus). The gene models were predicted with EuGene using a training dataset as described in[START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]. The assembly of the haploid genome of a 1 mating type (MvSl-1318-D) and the assembly of the haploid genome of a 2

assembled the Illumina genomes of the MvSl-1005 and the MvSl-140-01 strains that belonged to the MvSl Northwestern and MvSl Southwester genetic clusters, respectively

[START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]

, and that showed high read depth (233X and 89X, respectively). The assemblies of the MvSl-1005 and the MvSl-140-01 Illumina genomes were performed in SpAdes v3.11.0

[START_REF] Bankevich | SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing[END_REF]

with the following options: --only-assembler --careful -covcutoff auto using the whole set of trimmed Illumina reads of each genome. KmerGenie

[START_REF] Chikhi | Informed and automated k-mer size selection for genome assembly[END_REF]

was used to estimate the best k-mer length for genome de novo assembly.

The best predicted k values were 61 and 87 for the MvSl-1005 genome and the MvSl-140-01 genome, respectively. These assemblies had lower quality than the MvSl-1064 genome assembly (see summary in Supplementary Table

S2

) and they are available upon request. In order to obtain a high-quality genome assembly of a second MvSl strain, we sequenced separately the two haploid genomes (corresponding to the a 1 and a 2 mating types) of the diploid MvSl-1318 strain (collected on S. latifolia, Olomouc, Czech Republic, Coord. GPS : 49.569172 and 17.287057, previously shown mating type (MvSl-1318-T) are available at GenBank under the BioProject accession number PRJNA437556 (BioSample IDs SAMN08667584 and SAMN08667587, respectively). The raw reads are available at NCBI SRA under the BioProject accession number PRJNA437556 (BioSample IDs SRR6825078 and SRR6825079 for MvSl-1318-D and MvSl-1318-T, respectively

  [START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]. Read mapping identified genes present in the reference genomes that were absent from the Illumina genomes. We called missing fragments of more than 500 bp by two different methods, to increase the sensitivity of copy number variation calls. The read depth-based method involved scanning the genome for local variation in read depth and calls missing fragments in genomic regions with low read-depth values[START_REF] Abyzov | CNVnator: an approach to discover, genotype, and characterize typical and atypical CNVs from family and population genome sequencing[END_REF]. The split read-based method is based on the identification of read pairs for which only one read is mapped onto the genome (Supplementary figureS1;[START_REF] Ye | Pindel: a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads[END_REF][START_REF] Lin | Making the difference: integrating structural variation detection tools[END_REF]. We focused on missing fragments affecting more than 90% of gene lengths on autosomes, after masking gene models for repeated elements. Merging the missing fragments detected by the two methods resulted in the calling of 810 missing fragments in MvSl and 99 in MvSd. Most missing fragments were small (median length of 2.0 kb in MvSl and 3.1 kb in MvSd; Supplementary figureS2A). In both MvSl and MvSd, there was little overlap between methods concerning the missing fragments detected (18.7% and 12.5%, respectively). This was expected, as the two methods were based on different read mapping approaches. The overlap between missing fragments detected by the two methods covered a mean of 68 kb per individual in MvSl, and 19 kb in MvSd (Supplementary figureS2B). Most missing fragments affected a single gene (Supplementary figure

coverage Illumina genomes of 38 MvSl strains and 19 MvSd strains sampled across Europe (Supplementary Table

S1

;

[START_REF] Whittle | Recent and massive expansion of the mating-type-specific region in the smut fungus Microbotryum[END_REF][START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]

with high-quality reference genomes for these two species: MvSl-1064 and MvSd-1303 (Supplementary Table

S2

;

[START_REF] Badouin | Chaos of rearrangements in the mating-type chromosomes of the anthersmut fungus Microbotryum lychnidis-dioicae[END_REF]

  and 88% of the genes displaying presence-absence polymorphism encoded proteins of unknown function compared to 44% and 45% for all autosomal genes in MvSl and MvSd, respectively. Based on the available expression data for MvSl-1064 in three types of conditions, poor and rich in vitro medium and in planta[START_REF] Perlin | Sex and parasites: genomic and transcriptomic analysis of Microbotryum lychnidis-dioicae, the biotrophic and plant-castrating anther smut fungus[END_REF], we showed that the genes displaying presence-absence polymorphism in MvSl had lower levels of expression than other genes (comparisons of reads per kilobase of transcript per million mapped reads (RPKM) values; Kruskal-Wallis rank sum tests p-value < -tailed Fisher's exact test p-value < 2.2e-16). In both species, the genes displaying presence-absence polymorphism were also significantly shorter than other genes (difference in median gene length in MvSl strains= 928 bp; Kruskal-Wallis rank sum test p-value < 2.2e-16; difference in median gene length in MvSd strains = 852 bp; Kruskal-Wallis rank sum test p-value = 1.5e-07). Finally, the genes displaying presence-absence polymorphism included a higher proportion of genes from multiple-copy gene families (24% and 25% of all genes displaying presence-absence polymorphism in MvSl and MvSd, respectively) than other genes (7% and 9% in MvSl and MvSd, respectively).

2.2e-16 in all conditions tested). The genes displaying presence-absence polymorphism in MvSl included a high proportion of genes not expressed in planta (RPKM<1; 37%, versus 10% for other genes; one

The genes displaying presence-absence polymorphism and encoding proteins with predicted functions included a large proportion of genes associated with biological processes relating to DNA repair mechanisms and the cellular response to stress in MvSl and to glycoprotein metabolism and glycosylation in MvSd. Corresponding gene ontology terms were not significantly enriched after Bonferroni multiple testing correction, but were significantly enriched after Benjamini & Hochberg multiple testing correction (Supplementary Table

Downloaded from https://academic.oup.com/gbe/advance-article-abstract/doi/10.1093/gbe/evy089/4990910 by guest on 21 May 2018

Acknowledgments

PacBio sequencing was conducted at the IGM Genomics Center, University of California, San Diego, La Jolla, CA. This work was supported by the European Research Council (ERC) [starting grant GenomeFun 309403 to TG]; the French National Research Agency (ANR) [Gandalf ANR Grant ANR-12-ADAP-0009 to TG]; and a Marie Curie European grant [PRESTIGE-2016-4-0013 to FEH]. FEH received the Young Biological Researcher Prize from the Fondation des Treilles, created by Anne Gruner Schlumberger, which supports research in Science and Art

(http://www.les-treilles.com). JTB was supported by an IDEEV funding. (Branco et al. 2017). The strains MvSl-1318, MvSl-1005, and MvSl-140-01 were de novo assembled in this study (see Supplementary Material and Methods for details).

environment, which should be validated using functional studies in the near future.

Supplementary Material legends

Supplementary Text legends: Supplementary Text S1: Quality control procedure of detected gene presence-absence polymorphism using a comparison of high-quality genome assemblies of two Microbotryum lychnidis-dioicae (MvSl) strains. Supplementary Text S2: Detection of strain-specific genes and study of their population structure.

Supplementary Table legends:

Supplementary Table S1: Whole-genome sequencing data of the Microbotryum lychnidisdioicae (MvSl) populations and Microbotryum silenes-dioicae (MvSd) populations used in this study. Information were retrieved from [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF]) and [START_REF] Whittle | Recent and massive expansion of the mating-type-specific region in the smut fungus Microbotryum[END_REF]. [START_REF] Badouin | Chaos of rearrangements in the mating-type chromosomes of the anthersmut fungus Microbotryum lychnidis-dioicae[END_REF] and [START_REF] Branco | Evolutionary strata on young mating-type chromosomes despite the lack of sexual antagonism[END_REF]. Information on the genome assembly of the MvSd-1303 strain were retrieved from Supplementary Table S3: List of the genes affected by presence-absence polymorphism in Microbotryum lychnidis-dioicae (MvSl) and Microbotryum silenes-dioicae (MvSd).

Supplementary Table S2: Characteristics of the genome assemblies of the

Frequencies of the gene absence allele among all strains and functional annotation are shown. For the genes affected by presence-absence polymorphism in MvSl, the overlaps of genes with selective sweep regions identified in the Northwestern MvSl genetic cluster in [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF] are shown. S4: Significantly over-represented biological process gene ontology (GO) terms for genes affected by presence-absence polymorphism in Microbotryum lychnidis-dioicae (MvSl) and Microbotryum silenes-dioicae (MvSd). The "Effective GO term count" column is indicating the number of genes with the given GO term. GO term enrichment pvalues corrected for multiple testing with Bonferroni correction and with Benjamini & Hochberg correction [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] is presented.

Supplementary Table

Supplementary Figure legends

Supplementary figure S4: Distribution of allele counts for gene absence in MvSl isolates

belonging to the Northwestern genetic cluster according to their genomic location, within a selective sweep region or not. We used selective sweep regions identified in the Northwestern MvSl genetic cluster in [START_REF] Badouin | Widespread selective sweeps throughout the genome of model plant pathogenic fungi and identification of effector candidates[END_REF].