# Genome compartmentalization predates species divergence in the plant pathogen genus Zymoseptoria

Alice Feurtey, Cécile Lorrain, Daniel Croll, Christoph Eschenbrenner, Michael Freitag, Michael Habig, Janine Haueisen, Mareike Möller, Klaas Schotanus, Eva H Stukenbrock

► To cite this version:

Alice Feurtey, Cécile Lorrain, Daniel Croll, Christoph Eschenbrenner, Michael Freitag, et al.. Genome compartmentalization predates species divergence in the plant pathogen genus Zymoseptoria. BMC Genomics, 2020, 22 (1), 15 p. 10.1186/s12864-020-06871-w. hal-03361059

HAL Id: hal-03361059

https://hal.inrae.fr/hal-03361059

Submitted on 1 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Genome compartmentalization predates species divergence in the plant pathogen genus Zymoseptoria

Alice Feurtey<fnref n="1" />,<fnref n="2" />, Cécile Lorrain<fnref n="1" />,<fnref n="2" />,³* , Daniel Croll⁴, Christoph Eschenbrenner<fnref n="1" />,<fnref n="2" />, Michael Freitag⁵, Michael Habig<fnref n="1" />,<fnref n="2" />, Janine Haueisen<fnref n="1" />,<fnref n="2" />, Mareike Möller<fnref n="1" />,<fnref n="2" />,⁵, Klaas Schotanus<fnref n="1" />,<fnref n="2" />,⁶ and Eva H. Stukenbrock<fnref n="1" />,<fnref n="2" />

## Abstract

**Background:** Antagonistic co-evolution can drive rapid adaptation in pathogens and shape genome architecture. Comparative genome analyses of several fungal pathogens revealed highly variable genomes, for many species characterized by specific repeat-rich genome compartments with exceptionally high sequence variability. Dynamic genome structure may enable fast adaptation to host genetics. The wheat pathogen *Zymoseptoria tritici* with its highly variable genome, has emerged as a model organism to study genome evolution of plant pathogens. Here, we compared genomes of *Z. tritici* isolates and of sister species infecting wild grasses to address the evolution of genome composition and structure.

**Results:** Using long-read technology, we sequenced and assembled genomes of *Z. ardabiliae*, *Z. brevis*, *Z. pseudotritici* and *Z. passerinii*, together with two isolates of *Z. tritici*. We report a high extent of genome collinearity among *Zymoseptoria* species and high conservation of genomic, transcriptomic and epigenomic signatures of compartmentalization. We identify high gene content variability both within and between species. In addition, such variability is mainly limited to the accessory chromosomes and accessory compartments. Despite strong host specificity and non-overlapping host-range between species, predicted effectors are mainly shared among *Zymoseptoria* species, yet exhibiting a high level of presence-absence polymorphism within *Z. tritici*. Using *in planta* transcriptomic data from *Z. tritici*, we suggest different roles for the shared orthologs and for the accessory genes during infection of their hosts.

**Conclusion:** Despite previous reports of high genomic plasticity in *Z. tritici*, we describe here a high level of conservation in genomic, epigenomic and transcriptomic composition and structure across the genus *Zymoseptoria*. The compartmentalized genome allows the maintenance of a functional core genome co-occurring with a highly variable accessory genome.

**Keywords:** Genome evolution, Orphan genes, Effectors, Genome architecture, Accessory genes

<footnote>* Correspondence: lorrain@evolbio.mpg.de</footnote>

<footnote n="1">¹Environmental Genomics, Max Planck Institute for Evolutionary Biology, 24306 Plön, Germany</footnote>

<footnote n="2">²Environmental Genomics, Christian-Albrechts University of Kiel, 24118 Kiel, Germany</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
## Background

Co-evolution between plants and pathogens can drive rapid evolution of genes involved in antagonistic interactions [1]. In filamentous plant pathogens, rapid evolution may be fueled by highly dynamic genome architecture involving repeat-rich compartments such as gene-sparse islands of repetitive DNA and accessory chromosomes [2, 3]. These compartments can show a high plasticity revealed by a high extent of gene and/or chromosome presence-absence variation and structural variants, such as inversions, insertions and deletions [4, 5]. Several plant pathogenic fungi have isolate-specific chromosomes, so-called accessory chromosomes.

Accessory chromosomes are characterized by intra-species presence-absence polymorphism, low gene density, an enrichment of repetitive sequences and, in some species, a different histone methylation pattern [6, 7]. It has been shown that accessory chromosomes encode genes involved in virulence such as in the species *Fusarium solani*, *Fusarium oxysporum* and *Leptosphaeria maculans* [8–11]. Little is known about the evolutionary origin of accessory chromosomes although experimental evidence from the asexual species *F. oxysporum* shows that accessory chromosomes may be acquired horizontally as chromosomes can be transferred between distinct isolates by hyphal fusion [10]. Through such transfers, virulence determinants may be exchanged between clonal lineages as accessory chromosomes in this species were shown to encode host specific virulence determinants and transcription factors regulating their expression [12].

Genes involved in plant-pathogen interactions may diversify at a higher rate in repeat-rich genome compartments and thereby evolve new virulence specificity faster [3]. These genes encode secreted proteins, so-called effectors [1]. Most known effectors target diverse cellular compartments and molecular pathways, including immune response-related pathways [13, 14]. Genes encoding Carbohydrate-active enzymes (CAZymes) have also been associated to the pathogenic lifestyle of fungal plant pathogens, particularly through their role in plant-cell wall degradation [15]. Thus, some secreted CAZymes may be essential from the early infection stage, like penetration of plant tissue, to later stages such as the necrotrophic phase where the pathogen feeds from dead plant tissue [16]. Likewise, secondary metabolites are known to be involved in plant infection and contribute to virulence and the interaction with other plant-associated microorganisms [17, 18]. Many of these genes can be predicted either according to their composition and known protein domains or through machine learning methods [19]. Thereby, in-depth genome annotations have proven important to predict and compare the content of pathogenicity-related genes in plant pathogens, as well as their genomic localization for example in rapidly evolving genome compartments.

The ascomycete pathogen *Zymoseptoria tritici* has emerged as a model organism in evolutionary genomics of pathogens. This species originated in the Fertile Crescent during the domestication of its host, wheat [20]. Closely related species of *Z. tritici* have been collected from wild grasses in the Middle East providing an excellent resource for comparative genome analyses of closely related and recently diverged pathogen species. Comparative analyses of genome organization and gene content within and among *Zymoseptoria* species have previously revealed a wide distribution of accessory chromosomes and dynamic gene content [21, 22]. The haploid genome of the reference isolate IPO323 comprises thirteen core and eight accessory chromosomes [23]. Some of these accessory chromosomes may encode traits that impact virulence of the fungus, however no gene encoded on an accessory chromosome has so far been described as a virulence or avirulence determinant [24–29]. Interestingly, the accessory chromosomes in *Z. tritici* show a low transcriptional activity in vitro as well as *in planta* [30, 31]. This suppression of gene expression correlates with an enrichment of heterochromatin associated with the histone modification H3K27me3 on the accessory chromosomes [6, 32].

In the reference isolate IPO323, the accessory chromosomes comprise more than 11% of the entire genome assembly. To which extent such a high amount of accessory DNA is also found in genomes of other members of the *Zymoseptoria* genus has so far been unknown due to the lack of high-quality genome assemblies and large-scale population sequencing. Assemblies based on short-read data failed to recover complete sequence of accessory chromosomes and “orphan regions” due to their high repeat content [22]. The asset of genome assemblies based on long-read sequencing was demonstrated in detailed genome comparisons of *Z. tritici* isolates sequenced with PacBio long-read sequencing [28, 33]. Comparison these *Z. tritici* high-quality chromosome assemblies revealed the occurrence of “orphan regions” enriched with transposable elements and encoding putative virulence-related genes [28, 34].

In this study, we investigate the genomic architecture and variability among five *Zymoseptoria* species. Beside presenting a new and significantly improved resource for future genomic studies of these fungal pathogens, we specifically ask: 1) how conserved is the genome architecture among *Zymoseptoria* species? 2) can we identify accessory compartments in other *Zymoseptoria* isolates? 3) to which extent does variation in genome architecture reflect variation in gene content?

To answer these questions, we used high-quality assemblies based on long-read sequence data and new gene predictions in two isolates of *Z. tritici* (Zt05 and Zt10) and one isolate of each of the sister species, *Z. ardabiliae*, *Z. brevis*, *Z. passerinii*, and *Z. pseudotritici*.
We explore the core and non-core genome architecture of *Zymoseptoria* spp. combining genomic data with transcriptome and histone methylation data and relate this to core and accessory genome compartments. Furthermore, we compare the distribution of orthologous and non-orthologous genes in the *Zymoseptoria* genomes and one additional Dothideomycete species. Our analyses reveal an overall conserved genome architecture characterized by gene-rich core compartments and accessory compartments enriched in species-specific genes. Finally, we report a remarkably high extent of variation in presence-absence of protein coding genes in a eukaryote genome.

## Results

### De novo assemblies using long-read sequencing for six Zymoseptoria spp.

We sequenced and assembled the genome of the reference isolates of *Z. ardabiliae*, *Z. brevis*, *Z. pseudotritici* and *Z. passerinii* and the genomes of two *Z. tritici* isolates sampled in Denmark and Iran [30]. The obtained contigs were filtered based on base-quality confidence and read depth to ensure high quality of the final assemblies (see Methods). This filter removed a high number of contigs (between 17 and 58% of the total), but little overall length (between 0.4 and 2.6% of the total assemblies), indicating that most of the excluded contigs were of small size (Table S1 and S5). The best assemblies were of the two *Z. tritici* isolates comprising 19 and 30 contigs and the most fragmented was of *Z. passerinii* comprising 103 contigs (Fig. 1). The resulting

assembly lengths ranged from 38.1 Mb for *Z. ardabiliae* to 41.6 Mb for *Z. brevis*, which is comparable to the reference assembly length of *Z. tritici* (39.7 Mb) but larger than previous short-read based assemblies (Table 1; previous assemblies ranged from 31.5 Mb for *Z. ardabiliae* to 32.7 Mb for *Z. pseudotritici* [22, 23]). The assembly of the Iranian *Z. tritici* isolate Zt10 has telomeric repeats at the end of all contigs, indicating that each chromosome is completely assembled, comprising six accessory and thirteen core chromosomes. The assemblies for the Danish *Z. tritici* isolate (Zt05), *Z. brevis* (Zb87) and *Z. pseudotritici* (Zp13) contained, respectively, twelve, nine and five fully assembled chromosomes including both core and accessory chromosomes (Figure S1). The assemblies of the *Z. ardabiliae* (Za17) and *Z. passerinii* (Zpa63) genomes included no fully assembled chromosomes, but twelve and ten contigs respectively with telomeres at one of the ends (Table 1; Figure S1).

The transcriptome-based gene predictions for these new assemblies include between 10,528 and 12,386 protein-coding genes (Table S2). This range is consistent with the annotation of the reference genome IPO323 reporting 11,839 protein-coding genes [22]. We used Benchmarking Universal Single-Copy Orthologs (BUSCO) from the lineage dataset *Pezizomycotina* to evaluate the completeness of the assemblies and gene predictions [37]. The proportion of complete BUSCO genes identified in our assemblies were comparable to the one obtained for the reference genome of *Z. tritici*

<figure><img src="image_3.png" /></figure>
<table><caption>**Table 1** Metrics of genome assemblies and annotation</caption><thead><tr><th>Species</th><th>Zymoseptoria tritici</th><th>Zymoseptoria pseudotritici</th><th>Zymoseptoria brevis</th><th>Zymoseptoria ardabiliae</th><th>Zymoseptoria passerinii</th></tr></thead><tbody><tr><td>Isolate</td><td>Zt05</td><td>Zt10</td><td>Zp13</td><td>Zb87</td><td>Za17</td></tr><tr><td>Origin</td><td>Denmark</td><td>Iran, llam province</td><td>Iran, Ardabil province</td><td>Iran</td><td>Iran, Ardabil province</td></tr><tr><td>Host</td><td>Triticum aestivum</td><td>Triticum aestivum</td><td>Dactylis glomerata</td><td>Phalaris paradoxa</td><td>Lolium perenne</td></tr><tr><td>Year (of isolation)</td><td>2004</td><td>2001</td><td>2004</td><td>2004</td><td>2004</td></tr><tr><td>Contig number</td><td>30</td><td>19</td><td>42</td><td>29</td><td>50</td></tr><tr><td>Total length (bp)</td><td>41,240,984</td><td>39,248,105</td><td>40,312,446</td><td>41,586,671</td><td>38,100,668</td></tr><tr><td>Mean contig size (bp)</td><td>1,374,699</td><td>2,065,690</td><td>59,820</td><td>1,434,023</td><td>762,013</td></tr><tr><td>N50</td><td>2,454,671</td><td>2,925,395</td><td>2,115,121</td><td>2,744,794</td><td>1,156,695</td></tr><tr><td>L50</td><td>6</td><td>5</td><td>7</td><td>7</td><td>11</td></tr><tr><td>Contigs with telomeric repeats on both ends</td><td>12</td><td>19</td><td>5</td><td>9</td><td>0</td></tr><tr><td>Complete BUSCO genes (%)</td><td>98.4</td><td>98.7</td><td>98.5</td><td>97.0</td><td>98.2</td></tr><tr><td>Number of genes</td><td>12,386</td><td>11,991</td><td>11,661</td><td>11,480</td><td>11,463</td></tr><tr><td>Repeat content (%)</td><td>19.9</td><td>16.5</td><td>20.8</td><td>29.2</td><td>18.2</td></tr><tr><td>Percentage of non-core sequences compared to the IPO323 reference</td><td>12.49</td><td>8.69</td><td>4.27</td><td>4.6</td><td>3.45</td></tr><tr><td>NCBI Biosample</td><td>SAMN04494882</td><td>SAMN02981321</td><td>SAMN02981322</td><td>SAMN03294124</td><td>SAMN02981326</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>SAMN02981330</td></tr></tbody></table>

(97.8%, Table 1). The assessment of gene content completeness (see Methods) indicates that, despite more fragmented assemblies of *Z. ardabiliae* and *Z. passerinii*, the genomes are complete in terms of gene content and that the unassembled fragments are more likely to comprise repeats and not protein-coding genes.

Based on the whole-genome sequences and the predicted genes we reconstructed the phylogeny of the *Zymoseptoria* genus using the publicly available genome of *Cercospora beticola* as an outgroup [38, 39]. For both trees, the phylogenetic relationship of the *Zymoseptoria* species is in accordance with previously published phylogeny based on seven loci sequenced in multiples isolates (Fig. 1) [21].

### Genomes of Zymoseptoria spp. comprise accessory chromosomes and compartments but show overall high synteny

Next we addressed the extent of co-linearity of the *Zymoseptoria* genomes. Using coordinates of orthologous genes, we were able to reveal a high extent of synteny conservation among the five *Zymoseptoria* species and between the three isolates of *Z. tritici*, as depicted in Fig. 2 and S2. Based on this high extent of synteny and the prediction of telomeric repeats, we identified the correspondence of chromosomes between the reference genome of *Z. tritici* IPO323 and the other *Zymoseptoria* genomes (Fig. 2 and S2). *Z. brevis* and *Z. pseudotritici* share a near perfect synteny in their core

chromosomes, however, when compared to *Z. tritici*, *Z. brevis* and *Z. pseudotritici* have two large-scale inversions comprising roughly ~ 900 kb and ~ 1.2 Mb of chromosomes 2 and 6, respectively (Fig. 2, S2 and S3). Based on the phylogeny in Fig. 1, it is likely that these two events occurred after the divergence of *Z. tritici* from *Z. brevis* and *Z. pseudotritici*. Overall, we observe a higher extent of synteny conservation between *Z. brevis* and *Z. pseudotritici* compared to *Z. tritici* IPO323 (Fig. 2; S2 and S3).

In *Z. tritici*, core and accessory compartments have very distinct genomic features. It was previously shown that hallmarks of accessory regions in the reference isolate IPO323 include lower gene density, lower levels of H3K4me2 methylation and reduced gene expression [6, 31]. In the reference genome of IPO323, compartments with these genomic and epigenomic hallmarks represent either accessory chromosomes or specific regions of the core chromosomes. Here we find that the specific accessory hallmarks including low gene density, low expression, low H3K4me2 methylation and significant enrichment of species-specific genes (see description below) on the non-core contigs are found in genomic compartments throughout the genus (Table S3, Fig. 3 and S4).

In the genome of the reference *Z. tritici* strain, the compartments that exhibit the hallmark of accessory chromosomes includes a particular region of the core chromosome 7 of ~ 0.6 Mb (Fig. 3a) [6]. As previously suggested, we also here define this region as “accessory-
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Intra- and inter-species synteny conservation in *Zymoseptoria* genus. a) Intra-species synteny between the reference genome of *Z. tritici* IPO323 and the genome of the Iranian *Z. tritici* isolate Zt10. Each color represents a different chromosome as defined in the reference *Z. tritici* IPO323 genome, except for accessory chromosomes, which are in grey. The links represent a subsample of orthologous genes (subsampled 1:2 for the accessory chromosomes and 1:10 for the core chromosomes for clarity of the visual representation). Contigs are ordered according to their synteny to the reference genome IPO323. Telomeric repeats are indicated in orange b) Inter-species synteny between the reference genome of *Z. tritici* and the genome of *Z. brevis* Zb87. The arrows represent the large-scale inversions identified between the genomes of these two species</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Genome architecture of the reference genome *Z. tritici* IPO323 (a) and *Z. pseudotritici* Zp13 (b). The segments constituting the first circle represents the chromosomes of IPO323 (a) and contigs of Zp13 (b) ordered according to the synteny with the chromosomes of the reference genome. Tracks from the outside to the inside are heatmaps representing respectively: gene density along chromosomes/contigs; gene expression in vitro (TPM); H3K4me2 levels in vitro and species-specific gene density per 100 kb windows. The arrows indicate the location of the region on chromosome 7 (and the corresponding syntenic region in *Z. pseudotritici*) displaying accessory-like genomic and regulatory hallmarks. Telomeric repeats are indicated in orange</figcaption></figure>
like” because the region exhibits clearly two of the three above-mentioned criteria [6, 31]. The region has low gene expression, low H3K4me2 levels but unlike accessory chromosomes exhibits high gene density. The encoded genes are mostly species-specific. Interestingly, this particularly large accessory region is observed in several *Zymoseptoria* spp. (Fig. 3b; S4). Based on synteny plots, we recognize ~ 0.7 Mb of the contig 28 in *Z. pseudotritici* and ~ 0.6 Mb of the contig 17 for *Z. brevis*, corresponding to chromosome 7 of *Z. tritici* (Fig. 2 and S2) and sharing the same hallmarks of accessory chromosomes (Fig. 3 and S4, Table S3). For the two remaining sister species, the fragmentation of the assembly does not allow the identification of such pattern although we observe a similar tendency with respect to transcription and species-specific gene enrichment on contig 19 (~ 0.6 Mb) of *Z. ardabiliae* corresponding to a fragment of chromosome 7 in IPO323 (Figure S5).

We also identified other regions enriched in isolate-specific genes, thus defining orphan loci in the core chromosomes of both Zt10 and Zt05 (Fig. 3 and S4). We observed a region of ~ 0.2 Mb of contig 1 in the Iranian isolate Zt10 corresponding to the core chromosome 3 in the IPO323 genome with high content of isolate-specific genes (Figure S4; Table S3). We furthermore identified small segments with species-specific genes on the core chromosomes of the wild-grass infecting sister species including a ~ 0.3 Mb region of the contig 26 in *Z. brevis* and ~ 0.1 Mb of contig 30. Overall, we show that genome compartmentalization in core and accessory regions is an ancestral and shared trait among the *Zymoseptoria* species. This phenomenon generates highly variable compartments and defines loci that deviate from genome averages in terms of gene content, sequence composition and synteny conservation.

## Variable repertoires of effector candidate genes

To obtain gene annotations for the *Zymoseptoria* genome assemblies, we established a custom pipeline adapted from Lorrain and co-workers (Figure S5) [40]. Briefly, we use the consensus of three methods to predict gene product localization, then extract secreted proteins to further identify predicted effectors. This detailed functional annotation provided a catalog of predicted gene functions and cellular localizations (Figure S6). For each genome, a large proportion of genes could not be assigned to a protein function. 49.6% of genes in *Z. tritici* (*N* = 5953) and up to 71.8% of genes in *Z. pseudotritici* (*N* = 8373) lack a predicted function (i.e. proteins of unknown function; Figure S6A). A relatively consistent number of genes are predicted for each functional category among *Zymoseptoria* spp. (Figure S6A and B). Likewise, the numbers of gene products predicted to belong to the different subcellular

localizations are very similar (Figure S6C) across the whole genus, including secreted proteins. The difference between the minimal and maximal gene number for the different categories of subcellular localizations does not exceed 1.6X between species (Figure S6C). Overall, secretomes range from 7% of the genes predicted in *Z. passerinii* (*N* = 828) to 11% of genes in *Z. ardabiliae* (*N* = 1328 genes, Figure S6B).

We further investigated the number and distribution of genes predicted to encode proteins with a pathogenicity-related function, such as secondary metabolites, CAZymes and predicted effectors (Figure S1 and S6B). Genes involved in the synthesis of secondary metabolites are typically organized in clusters, with genes participating in the same biosynthetic pathway grouping together at a genomic locus (Shi-Kunne et al. 2019). The number of biosynthetic gene clusters (BGC) ranges from 25 in *Z. ardabiliae* and *Z. passerinii* to 33 in the IPO323 reference genome and includes from 305 to 471 predicted genes (Figure S1). The only BGC identified in a non-core contig is a non-ribosomal peptide synthetase BGC found on the contig 38 of *Z. brevis* which has no orthologous cluster detected in any of the other *Zymoseptoria* genomes (Figure S1). We identified between 454 and 515 CAZyme genes in the *Zymoseptoria* species. Both BGCs and CAZymes are almost exclusively found on the core chromosomes (Figure S1). The only exceptions are a CAZyme encoding gene found on chromosome 14 in *Z. tritici* IPO323 and Zt05, and a CAZyme encoding gene on the putative accessory contig 38 of *Z. brevis* (Figure S1). These two genes encode for a beta-glucosidase and a carboxylic-ester hydrolase, respectively.

In contrast to the high conservation of CAZyme and BGC gene content among the *Zymoseptoria* genomes, we find that predicted effector genes exhibit a large variation in gene numbers between genomes (Figure S6B). In fact, the predicted effector gene repertoire in *Z. ardabiliae* (*N* = 637) is three times higher compared to *Z. brevis* (*N* = 206). Interestingly, the three *Z. tritici* isolates also vary considerably in their predicted effector repertoires. The reference isolate IPO323 has a reduced set of predicted effector genes (*N* = 274) compared to Zt05 and Zt10 that encode approximately 30% more predicted effector genes (*N* = 417 and *N* = 403, respectively, Figure S6B). Despite the high variability, the predicted effector genes are mostly located on core chromosomes and none of the five *Zymoseptoria* species have more than ten predicted effector genes located on accessory chromosomes (Figure S1).

## The accessory genes of *Z. tritici* are shared with the closely related wild-grass infecting species

To further characterize variation in gene content among the five *Zymoseptoria* species, we identified orthologous
genes (i.e. orthogroups) from the gene predictions. We categorized 22,341 gene orthogroups identified in the seven *Zymoseptoria* genomes and in *C. beticola* according to their distribution among fungal genomes (Fig. 4a). The core orthogroups, which are genes present in all eight genomes, represent around 30% of all orthogroups (N = 6698). The *genus-specific* orthogroups, shared between several *Zymoseptoria* spp. but not found in the *C. beticola* genome, represent 45% of the orthogroups (N = 9955; ranging from 2066 to 3212 per species). Among the *genus-specific* orthogroups, 1100 are found in all *Zymoseptoria* genomes (Fig. 4a), whereas all others show presence-absence polymorphisms within the genus. A total of 2476 *species-specific* orthogroups (ranging from 552 to 1191 per species) are found only in individual species. Among the *species-specific* genes, 205 orthogroups (N<sub>genes</sub> = 414 to 562) are found in all three *Z. tritici* genomes while the *isolate-specific* genes in *Z. tritici* represent 391 (Zt10) to 792 (IPO323) genes.

Comparing the three *Z. tritici* isolates independently from the other species, we observe extensive gene presence-absence polymorphisms between the three isolates: 1540 orthogroups are identified in only two

strains and 2522 are found in only one (Fig. 4b). The number of genes showing presence-absence variation is striking compared to the 10,098 core genes in *Z. tritici* as these genes comprise almost 30% of all predicted genes. Interestingly, we show that the number of orthogroups detected as *isolate-specific* is much larger when the comparison includes only members of the same species than when the other species are included (1035, 849 and 638 vs 792, 659 and 391 genes for IPO323, Zt05 and Zt10 respectively; Fig. 4a and b). This indicates that a large part of the accessory gene content in *Z. tritici* is shared among the sister species, and highlights the importance of including sister species when establishing core and accessory gene content.

Interestingly, we show that predicted effectors are enriched among the *genus-specific* genes but not among the *species-specific* or *isolate-specific* gene categories (with the exception of *Z. ardabiliae*). Fifty-six percent of predicted effectors in *Z. ardabiliae* and up to 78% of predicted effectors in *Z. pseudotritici* are shared with at least one of the other five *Zymoseptoria* species (Fig. 4c). Indeed, 427 predicted effector orthogroups are found in at least two genomes. However, only 47 (10% of the total

<figure><img src="image_6.png" /></figure>
predicted effector orthogroups N = 474) are found in all seven *Zymoseptoria* genomes. Among the predicted effectors shared by *Z. tritici*, and at least one other *Zymoseptoria* species, 32% (N = 112 of 352) are present in all three *Z. tritici* isolates while 68% (N = 240 of 352) show presence-absence polymorphisms in at least one of the three isolates. These results indicate that the majority of these shared predicted effectors are actually accessory (i.e. presence-absence polymorphism) in *Z. tritici*.

**Among in planta differentially expressed genes, species-specific are more expressed than core genes**

Finally, we addressed the functional relevance of accessory and orphan genes in *Z. tritici* by analyzing gene expression patterns. We used previously published *in planta* expression data of three *Z. tritici* isolates [30]. The expression profiling was obtained from four subsequent infection stages including infection establishment (stage A), biotrophic colonization (stage B), the transition from biotrophic to necrotrophic phase (stage C) and necrotrophic colonization (stage D) [30]. We sorted *in planta* expression data into two different

infection phases: the biotrophic phase and the necrotrophic phase, a separation supported by principal component analysis of normalized DESeq2 counts (Figure S7). Furthermore, we distinguished gene expression of the above-defined categories (*core genes*, *genus-specific*, *species-specific*, and *isolate-specific*). We compared expression levels by mapping RNA-seq reads to the genomes of IPO323, Zt05 and Zt10, using normalized read mappings to transcript per million. We tested differences among gene categories using pairwise comparisons with a Kruskal-Wallis test (Fig. 5). Overall, we find that gene expression of the *species-specific* and *isolate-specific* genes is significantly lower in IPO323 and Zt10, but not in Zt05 (Kruskal-Wallis *p*-value < 0.05). *Species-specific* and *isolate-specific* gene median expression ranges from 3.2 to 5.6 TPM in IPO323 and Zt10 while median expression of core genes is 12.1 and 10.9, respectively. The Zt05 expression profile does not follow the same trend: the *core genes* are the lowest expressed gene category (8.9 median TPM), while *genus-; species-* and *isolate-specific* genes showed higher transcription levels (12.0; 14.4 and 13.5 median TPM respectively, Kruskal-Wallis *p*-value < 0.05).

<figure><img src="image_7.png" /></figure>
In contrast, we observe a significantly higher expression of the *species-specific* and *isolate-specific* genes for all three isolates (Table S4; Kruskal-Wallis *p*-value < 0.05) when comparing the expression of genes that are differentially expressed (DEGs; DESeq2 *p*-adjusted < 0.05) between the biotrophic and necrotrophic phases. *Species-specific* and *isolate-specific* DEGs are higher expressed *in planta* than the *core* and *genus-specific* genes (Fig. 5b; Kruskal-Wallis *p*-value < 0.05). The expression patterns of DEGs with different levels of specificity present a consistent pattern in all three isolates (Table S4). Overall, this comparison reveals a potential functional relevance of accessory genes, which are up-regulated during infection of *Z. tritici*.

## Discussion

In this study we present a new resource of high-quality whole genome assemblies and gene annotations for the fungal grass pathogens *Z. ardabiliae*, *Z. brevis*, *Z. passerinii*, *Z. pseudotritici*, and two isolates of the wheat pathgoen *Z. tritici*. This new dataset provides a valuable resource for detailed analyses of genome architecture and evolutionary trajectories in this group of plant pathogens. Here, we conduct some detailed comparative analyses of genome architecture and show a considerable extent of variation in sequence composition during the recent evolution of the *Zymoseptoria* lineages. We show that genome compartmentalization and accessory chromosomes represent shared ancestral traits among these pathogen species.

We identify extensive presence-absence variation of protein coding genes in genomes of the five *Zymoseptoria* species consistent with the variable gene repertoire already reported for one of the species, *Z. tritici* [28]. Furthermore, the different species, share a particular genomic architecture that comprises specific accessory genome compartments. In spite of this variation, we observe an overall conserved synteny of the core chromosomes. In the *Zymoseptoria* genomes, we observe gene-dense, actively transcribed and H3K4me2-enriched compartments associated with most of the core chromosomes. These compartments are clearly distinguishable from gene-sparse, non-transcribed and H3K4me2-deprived compartments. Based on previous analyses of accessory chromosomes in *Z. tritici*, we here consider this pattern as a specific hallmark of accessory genome compartments in the genus *Zymoseptoria* beyond only in *Z. tritici* [6]. We hypothesize that these compartments likely represent accessory chromosomes in the different *Zymoseptoria* species.

We also identify accessory signatures in core chromosomes, including the previously described right arm of chromosome 7 [6]. Although this region has not been reported to share the same extent of presence-absence polymorphism as the accessory chromosomes, a considerably smaller chromosome 7 was reported in a single *Z. tritici* isolate originating from Yemen [33]. Here we

show that the region homologous to chromosome 7 in the other *Zymoseptoria* species also exhibits accessory compartment hallmarks. Our results support the occurrence of a past chromosome fusion, but hereby show that it very likely occurred prior to the divergence of the species (estimated to date tens of thousands of years [21]). The specific genomic and epigenetic features have remained stable through speciation and evolutionary time.

In this study, we confirm previously reported genome comparisons showing that gene content in *Z. tritici* is highly variable [28]. We further extended the identification of orthologs throughout the whole *Zymoseptoria* genus. Thereby, we show that more than 25% of the genes identified as *isolate-specific* in a comparison including only *Z. tritici* isolates are actually present in the wild-grass infecting sister species. This observation suggests that a large proportion of the accessory genome of *Z. tritici* is not specific to this species. Instead, the accessory genome content is shared among *Zymoseptoria* species. The proportion of accessory *Z. tritici* genes shared with other *Zymoseptoria* species was found to be the highest in the Iranian isolate, which is the only isolate sympatric with the four sister-species. A likely explanation for this observation would be inter-specific gene flow, which would allow the different wild species to exchange genes with sympatric *Z. tritici* isolates. This new finding is consistent with recent findings from population genomic data studies revealing extensive introgression between *Zymoseptoria* species [41, 42]. Our observation opens new perspectives for further analysis to understand how inter-specific gene flow has affected the evolution of the accessory genome of *Z. tritici*.

The genes with predicted functions and, in particular, functions related to pathogenicity are largely shared among species in the *Zymoseptoria* genus. Although the lifestyles of the wild-grass infecting *Zymoseptoria* are poorly understood, the species share major features of their lifestyles. Thus, as expected, we find similar CAZymes and BGC contents across the genomes studied here. In *Zymoseptoria*, most of the predicted effectors are shared among all species, although they show presence-absence variation. In the *Botrytis* genus (Dothideomycetes), sister species infecting different hosts share effectors with confirmed functions [43]. Likewise, in *Microbotrym*, a fungal plant pathogen including several species specialized on different hosts a large set of genes encoding conserved and shared secreted proteins was identified. This repertoire of effectors is hypothesized to include traits relevant for pathogenicity among *Microbotrym* lineages [44]. We hypothesize that the different specificity levels reflect functional differences in the effector repertoire of *Zymoseptoria*. Predicted effector
genes conserved across the *Zymoseptoria* genus are likely core pathogenicity factors potentially targeting key plant defense mechanisms common to all of the grass hosts [45]. Variation in the composition of effector genes in plant pathogen genomes, including the presence of *species-specific* and *isolate-specific* may reflect different host specificities and rapid evolution of these genes [45]. Here we find that only a fraction of the *genus-specific* predicted effector genes is shared among all *Zymoseptoria* species; the majority shows presence-absence polymorphisms suggesting that a variable effector repertoire is an ancestral trait in these plant pathogens. *Z. ardabiliae* has been isolated from leaves of distantly related grass species in Iran, including *Lolium* spp., *Elymus repens*, and *Dactylis glomerata* and potentially resulting in a broader *species-specific* effector repertoire of *Z. ardabiliae* compared to the other species [21].

Consistent with a previous study [28], we found that *Z. tritici* core genes are generally more expressed compared to accessory genes *in planta*. Core genes are more likely to encode essential functions, which could explain higher expression pattern during infection. Differentially expressed genes that are specifically induced during the course of the infection are very likely to have functions essential to pathogenicity of the fungus. Interestingly, here we show that within the genes differentially expressed between the biotrophic and the necrotrophic phases of infection, *isolate- and species-specific* genes have higher expression levels than core genes. These *isolate- and species-specific* genes could be functionally important and regulate functions linked to infection success in the biotrophic phase or to leaf colonization in the necrotrophic phase. Since these genes show presence-absence polymorphisms in the genus and in the *Z. tritici* species, they could represent a reservoir for possible adaptations to host species, host cultivars or local environments.

## Conclusion

We investigated the genomic architecture in a genus of plant pathogens, including the economically relevant wheat pathogen *Z. tritici*. Comparing genome content and genome structure, we identified a large shared predicted effector repertoire characterized by inter- and intraspecies presence-absence polymorphisms. Major features of genomic, transcriptomic and epigenetic compartmentalization, distinguishing accessory and core compartments, were shared among wheat and wild-grass infecting *Zymoseptoria* species. We conclude that compartmentalization of genomes is an ancestral trait in the *Zymoseptoria* genus.

## Methods

### Fungal material, DNA extraction and sequencing

Details regarding the individual *Zymoseptoria* isolates can be found in Table 1. For genomic data we used the

three *Z. tritici* isolates IPO323 (reference), Zt05 and Zt10, one *Z. ardabiliae* isolate (Za17), one *Z. brevis* isolate (Zb87), one *Z. passerinii* isolate (Zpa63), and the *Z. pseudotritici* isolate (Zp13). For transcriptomic and epigenomic data we used *Z. tritici* Zt09 (IPO323 ΔChr18) a derivate of the reference isolate IPO323 deleted with the chromosome 18 [31].

Long read assemblies of the *Z. tritici* isolates Zt05 and Zt10 were described and published previously [30]. For DNA extraction and long read sequencing cultures of *Z. pseudotritici*, *Z. ardabiliae*, *Z. brevis* and *Z. passerinii* were maintained in liquid YMS medium (4 g/L yeast extract, 4 g/L malt extract, 4 g/L sucrose) at 200 rpm and 18 °C. DNA extraction was conducted as previously described [26]. PacBio SMRTbell libraries were prepared using DNA extracted from single cells based on a CTAB extraction protocol [30, 46]. The libraries were size selected with an 8-kb cutoff on a BluePippin system (Sage Science).

After selection, the average fragment length was 15 kb. Sequencing of the isolates Za17, Zb87, and Zp13 was run on a PacBio RS II instrument at the Functional Genomics Center, Zurich, Switzerland. Sequencing of the Zpa63 isolate was performed at the Max Planck-Genome-Centre, Cologne, Germany.

### Genome assembly, and repeat and gene predictions

For each isolate, we assembled the genome de novo using SMRT Analysis software v.5 (Pacific Bioscience) with two sets of parameters: default parameters and “fungal” parameters. We chose the best assemblies generated by comparison of all assembly statistics produced by the software Quast such as the number of finished contigs, the size of the assembly and the N50 [47]. Summary statistics for each assembly can be found in Table 1. In order to exclude poor quality contigs from the raw assemblies, we filtered out the contigs with less than 1.5X and more than 2X median read coverage as these might be unreliable from lack of data or because they contain only repeated DNA [34]. This filter removed a high number of contigs, i.e., between 58 and 17% of contigs, but a only a small quantity of base pairs as compared to the genome size (Table S1; S5). In order to identify the number of fully assembled chromosomes or chromosome arms, we investigated the presence of telomeric repeats (“CCCTAA”) in the assembled contigs using bowtie2 and recorded the presence of more than six repeats at the contig extremities [48–50]. We reported for each contig the number of such blocks of telomeric repeats and considered a contig flanked by these repeats on both sides to be a fully assembled chromosome.

We next used the REPET package to annotate the repeat regions of *Z. ardabiliae* Za17, *Z. brevis* Zb87, *Z.*
pseudotritici Zp13, Z. passerinii Zpa63, and the three Z. tritici isolates IPO323, Zt05 and Zt10 (https://urgi.ver sailles.inra.fr/Tools/REPET [51, 52]). For each genome, we annotated the repetitive regions as follows: we first identified repetitive elements in each genome using TEdenovo following the developer's recommendations and default parameters. The library of identified consensus repeats was then used to annotate the respective genomes using TEannot with default parameters.

We used previously published RNA sequencing data to increase the quality of the gene prediction and combined three distinct methodologies [22, 30, 32]. As a first approach, we used GeneMark-ES for an ab initio prediction using the option “--fungus” [53]. Our second and third approaches both used RNA-seq data. For this, we first trimmed the reads using Trimmomatic [54]. We mapped the filtered and trimmed reads to the newly assembled genomes using hisat2 [50] and used the BRAKER1 pipeline to predict genes for each genome using the fungus flag [55]. BRAKER applies GeneMark-ET and Augustus to create the first step of gene predictions based on spliced alignments and to produce a final gene prediction based on the best prediction of the first set [56, 57]. For our third approach, the RNA-seq reads were separately assembled into gene transcripts using Trinity [58]. These were aligned using PASA and EEvidence Modeler to produce consensus gene models from the two independent predictions and the de novo assembled transcripts [59]. Gene counts, length and other summary statistics presented in Table S1 and S3 were obtained using GenomeTools [60] and customs scripts (https://gitlab.gwdg.de/alice.feurtey/genome_archi tecture_zymoseptoria).

The predicted gene sequences were the basis for an evaluation of the completeness of the assembly and gene prediction by the program BUSCO v.3 [37]. We used this method with the lineage dataset *Pezizomycotina*. The predicted genes were also used to create a phyl- ogeny with the online implementation of CVtree3, using kmer sizes of 6 and 7 as recommend for fungi [36]. We generated a second tree with the whole assemblies, esti- mating a distance matrix using the andi software [35].

We predicted orthologs between the newly assembled genomes, the reference *Z. tritici* genome and the refer- ence *Cercospora beticola* genome, a related Dothideomy- cete, which we used as outgroup to identify genes with orthologs restricted to the *Zymoseptoria* genus [21]. For this, we used the software PoFF [39, 61] which takes into account synteny information in the analyses of similarity inferred by the program Proteinortho [39]. These orthogroups were used to visualize synteny between genomes using Circos [62].

The whole-genome assemblies were used to create a matrix distance with the software andi, from which we generated a tree [35]. A second tree was generated from

the gene prediction with the online implementation of CVtree3 [36].

## Functional annotations

We used several tools to predict the putative functions for the gene models. First, we used the eggnog-mapper which provide COG, GO and KEGG annotations [63]. The online resource dbCAN2 was run to identify carbohydrate-active enzymes (CAZymes) [64]. Finally, for each genome, we used Antismash v3 (fungal version) to detect biosynthetic gene clusters (Figure S1 [65]).

Additionally, we designed a pipeline to predict protein cell localization and to identify effector candidates. The pipeline for effector prediction is outlined in Figure S1 and includes the software DeepLoc [66], SignalP [67], TargetP [68], phobius [69] and TMHMM [70, 71], which predict the cellular location, the peptide signals and whether proteins are transmembrane. Effector candi- dates were identified with EffectorP v2 which uses both a new machine learning approach and more complete databases to improve effector prediction compared to the previous version [19]. The pipeline also includes software which are specifically targeted to annotate plant pathogenic functions, namely the program ApoplastP [72] and LOCALIZER [73]. We wrote wrappers scripts, which run the software and create consensus between the different prediction tools providing one command line from the user. These scripts are available at https:// gitlab.gwdg.de/alice.feurtey/genome_architecture_zymo septoria. Briefly, we gathered outputs of several software to predict the cellular location, transmembrane domain and secretion and created a consensus based on the different output to prevent the pitfalls of any one of these methods. From this consensus, we extracted the gene products predicted to be secreted and without a transmembrane domain. The comparisons of genes functions repartition were done by combining predic- tions of COG categories, secondary metabolite genes with pathogenicity-related gene functional categories such as CAZymes and effector predictions.

## Gene expression analyses

To update expression profiles on the new genome assemblies and new gene predictions of the three *Z. tritici* isolates IPO323, Zt05 and Zt10, we used previ- ously generated RNA-seq data from *in planta* and *in vitro* growth [30, 32]. The *in planta* RNA-seq data was obtained from infected leaves at four different stages corresponding to early and late biotrophic and necro- trophic stages of the three *Z. tritici* isolates [30]. Strand- specific RNA-libraries were sequenced using Illumina HISeq2500, with 100pb single-end reads for a total read number ranging from 89.5 to 147.5 million reads per sample. This data was previously analyzed [30], using
gene predictions generated from an Illumina-based assembly [22]. The reads were here mapped on the new assemblies of Zt05 and Zt10 and the reference genome of IPO323 after trimming. We used the DESeq2 R package to determine differential gene expression during *in planta* infection, considering only two infection stages; biotrophic and necrotrophic [74]. Gene expression was assessed as Transcript per Million (TPM). Briefly, TPM is calculated by normalizing read counts with coding region length resulting in the number of reads per kilobase (RPK). RPK total counts per sample are then divided by 1 million to generate a “per million” scaling factor. We calculated the coding region length of each gene with GenomicFeatures R package using the function called “exonsBy” [75]. For gene expression analyses, we further filtered our gene predictions to remove any predicted transposases and other TE-related annotations based on the Eggnog mapper annotations.

## ChIP-sequencing and data analysis

Z. ardabiliae (Za17) and Z. pseudotritici (Zp13) cells were grown in liquid YMS medium for 2 days at 18 °C until an OD$_{600}$ of ~1 was reached. Chromatin immunoprecipitation and library preparation were performed as previously described [76]. We sequenced two biological and two technical replicates per isolate and used antibodies against the euchromatin histone mark H3K4me2 (#07–030, Merck Millipore). Sequencing was performed at the OSU Center for Genome Research and Biocomputing (Oregon State University, Corvallis, USA) on an Illumina HiSeq2000 to obtain 50-nt reads. The data was quality-filtered using the FastX toolkit (http://hannonlab.cshl.edu/fastx_toolkit/), mapping was performed using bowtie2 [77] and peaks were called using HOMER [78]. Peaks were called individually for each replicate, but only peaks that were detected in all replicates were considered and merged for further analysis. Merging of peaks and genome wide sequence coverage with enriched regions was assessed using bedtools [49].

**Additional file 5: Table S5.** Summary metrics of *Zymoseptoria* genome sequencing coverage, assembly size and median coverage per contig.

**Additional file 6: Figure S1.** Plant-associated genes compartmentalization along the chromosomes. The first track represents core (dark grey) and accessory (light grey) chromosomes/contigs. Telomeric repeats are indicated in orange. Circles from outside to inside represent the position of: predicted effector genes (blue), biosynthetic gene clusters (BGC, green) and CAZymes (yellow) per 100 kb windows, respectively.

**Additional file 7: Figure S2.** A) Intra-species synteny between the reference genome of Z. tritici and the genome of the Z. tritici isolate Zt05. Each color represents a different chromosome as based on the reference Z. tritici genome, except for accessory chromosomes, which are in grey. The connecting lines represent orthologs between each genome. The track between the chromosomes and connecting lines are predicted effector genes. Telomeric repeats are indicated in orange. B) Inter-species synteny between reference genome of Z. tritici and Z. pseudotritici. The arrows represent the large-scale inversions identified between the genomes of these two species.

**Additional file 8: Figure S3.** Inter-species synteny between genomes of Z. pseudotritici and Z. brevis. A) High level of inter-species synteny at the whole genome level between Z. pseudotritici (dark grey) and Z. brevis (blue). The arrows indicate the contigs identified as carrying rearrangement with Z. tritici but not between Z. pseudotritici and Z. tritici. B – D) Zoom-in on the contigs highlighted with the arrow in part A, showing the details of the pair-wise synteny of both Z. pseudotritici (Zp13) and Z. brevis (Zb87) with the Z. tritici reference (IPO323) chromosomes.

**Additional file 9: Figure S4.** Genome architecture in A) Z. tritici Zt05, B) Z. tritici Zt10. C) Z. brevis Zb87 and D) Z. ardabiliae Za17. Circles from the outside to the inside represent respectively: gene density along chromosomes/contigs; gene expression in vitro (TPM); H3K4me2 distribution in vitro (only for Za17) and species-specific gene density per 100 kb windows. Contigs are ordered based on synteny with chromosomes of the reference strain IPO323. Telomeric repeats are indicated in orange.

**Additional file 10: Figure S5.** Simplified diagram of the pipeline used to predict the functions and subcellular localization of gene model products.

**Additional file 11: Figure S6.** Functional gene categories in *Zymoseptoria spp.* genomes. A) The number of genes in Egg mapping COG categories in addition to Effectors and CAZymes. B) Pathogenicity-related genes of interest: secreted proteins, predicted effectors, secreted CAZymes and non-secreted CAZymes. C) Subcellular localization of predicted gene products.

**Additional file 12: Figure S7.** Principal component analysis of DESeq2 rlog transformed expression data. Principal component analysis (PCA) of *Zymoseptoria tritici* IPO323, Zt05 and Zt10 transcripts levels measured in biotrophic (red) and necrotrophic (blue) stages using RNA-sequencing. Reads detected per transcript (counts) were normalized using the size factor method used by DESeq2 package [73]. The PCA plot places biological replicates (four replicates per stage) along the two first PC axes explaining 58 to 65% of the variance (x-axis) and 15 to 25% (y-axis) of the variance within samples, respectively.

## Supplementary information

Supplementary information accompanies this paper at https://doi.org/10.1186/s12864-020-06871-w.

**Additional file 1: Table S1.** Summary metrics of assemblies filtering process in function of contig coverage quality, presence of telomeric repeats and mitochondrial DNA.

**Additional file 2: Table S2.** Summary of gene expression in Z. tritici isolates during infection stages.

**Additional file 3: Table S3.** Enrichment of species-specific and isolate-specific genes per chromosome/contig.

**Additional file 4: Table S4.** Summary metrics of Z. tritici IPO323, Zt10 and Zt05 expressed genes during *in planta* infection stages: biotrophic versus necrotrophic.

## Abbreviations

BGC: Biosynthetic gene cluster; BUSCO: Benchmarking universal single-copy orthologs; CAZymes: Genes encoding carbohydrate-active enzymes; COG: Clusters of orthologous groups; DEG: Differentially expressed gene; GO: Gene ontology; H3K27me3: Histone H3, 27th lysine residue trimethylation; H3K4me2: Histone H3, 4th lysine residue di-methylation; KEGG: Kyoto encyclopedia of genes and genomes; PacBio: Pacific biosciences; RPK: Reads per kilobase; TPM: Transcripts per million

## Acknowledgements

Not applicable.

## Authors' contributions

AF, CL and MM performed data and results analyses. EHS, AF and CL contributed to the design and implementation of the research. AF, CL, CE and DC performed genome assemblies. MH, JH, MM, MF and KS performed
the experimental procedures. All authors contributed to the writing of the manuscript. All authors read and approved the manuscript.

## Authors' information

Present Address: Pathogen Evolutionary Ecology, ETH Zurich, 8092 Zürich, Switzerland.

Present Address: Laboratory of Evolutionary Genetics, Institute of Biology, University of Neuchâtel, 2000 Neuchâtel, Switzerland.

## Funding

CL is funded by the Institut national de la recherche agronomique (INRA) in the framework of a "Contrat Jeune Scientifique" and by the Labex ARBRE (Lab of Excellence ARBRE). Genome research in the group of EHS is funded by the DFG Priority Program SPP1819 and the Canadian Institute for Advanced Research (CIFAR). AF is funded by the DFG Priority Program SPP1819.

## Availability of data and materials

The datasets generated and analyzed during the current study are available in the NCBI Short Read Archive (https://www.ncbi.nlm.nih.gov/bioproject/?term=) under the BioProject accession numbers: PRJNA638605, PRJNA639021 (Zpa63), PRJNA638553 (Zb87), PRJNA638515 (Zp13), PRJNA638382 (Za17), PRJNA414407 (Zt05 and Zt10). All the data supporting the findings of this study are openly available at https://doi.org/10.5281/zenodo.3820378. The assembled genomes can be found at https://doi.org/10.5281/zenodo.3820378. The gene annotations are deposited at https://doi.org/10.5281/zenodo.3820378. The functional annotation pipeline, additional scripts and command lines used to create the results presented in this manuscript can be found at https://gitlab.gwdg.de/alice.feurtey/genome_architecture_zymoseptoria. In planta RNA-seq datasets for Z. tritici IPO323, Zt05, Zt10 from [30] are available at the accession number NCBI Gene Expression Omnibus GSE106136 and in vitro RNA-seq data for Za17, Zp13 and Zb87 from [22] are available at the accession number NCBI BioProjects PRJNA277173, PRJNA277174, and PRJNA277175. The genome sequence of the reference isolate IPO323 available at: http://genome.jgi.doe.gov/Mycgr3/Mycgr3.home.html. The ChIP-seq dataset are available under the BioProject accession PRJNA640982 for Za17 and Zp13. The ChIP-seq datasets for the reference genome IPO323 from [6] are available under the accession number PRJNA286790.

## Ethics approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Competing interests

Pr. Eva Stukenbrock is a member of the editorial board for the BMC Genomics journal.

## Author details

¹Environmental Genomics, Max Planck Institute for Evolutionary Biology, 24306 Plön, Germany. ²Environmental Genomics, Christian-Albrechts University of Kiel, 24118 Kiel, Germany. ³INRA Centre Grand Est - Nancy, UMR 1136 INRA/Université de Lorraine Interactions Arbres/Microorganismes, 54280 Champenoux, France. ⁴Laboratory of Evolutionary Genetics, Institute of Biology, University of Neuchâtel, 2000 Neuchâtel, Switzerland. ⁵Department of Biochemistry and Biophysics, Oregon State University, Corvallis, OR, USA.

⁶Department of Molecular Genetics and Microbiology, Duke University, Duke University Medical Center, Durham, NC 27710, USA.

Received: 13 December 2019 Accepted: 26 June 2020

Published online: 26 August 2020

## References

1. Cook DE, Mesarich CH, Thomma BPHJ. Understanding plant immunity as a surveillance system to detect invasion. Annu Rev Phytopathol. 2015;53:541–63.

2. Rouxel T, Grandaubert J, Hane JK, Hoede C, van de Wouw AP, Couloux A, et al. Effector diversification within compartments of the Leptosphaeria maculans genome affected by repeat-induced point mutations. Nat Commun. 2011;2:202.

3. Dong S, Raffaele S, Kamoun S. The two-speed genomes of filamentous pathogens: waltz with plants. Curr Opin Genet Dev. 2015;35:57–65.

4. Todd RT, Wikoff TD, Forche A, Selmecki A. Genome plasticity in *Candida albicans* is driven by long repeat sequences. Elife. 2019;8.

5. Mehrabi R, Mirzadi Gohari A, Kema GHJ. Karyotype variability in plant-pathogenic Fungi. Annu Rev Phytopathol. 2017;55:483–503.

6. Schotanus K, Soyer JL, Connolly LR, Grandaubert J, Happel P, Smith KM, et al. Histone modifications rather than the novel regional centromeres of *Zymoseptoria tritici* distinguish core and accessory chromosomes. Epigenetics Chromatin. 2015;8:41.

7. Fokkens L, Shahi S, Connolly LR, Stam R, Schmidt SM, Smith KM, et al. The multi-speed genome of Fusarium oxysporum reveals association of histone modifications with sequence divergence and footprints of past horizontal chromosome transfer events. bioRxiv. 2018:465070.

8. Miao VP, Covert SF, Vanetrent HD. A Fungal Gene for Antibiotic Resistance on a Dispensable (" B ") Chromosome. Science 1991;254.

9. Temporini E, VanEtten H. An analysis of the phylogenetic distribution of the pea pathogenicity genes of Nectria haematococca MPV supports the hypothesis of their origin by horizontal transfer and uncovers a potentially new pathogen of garden pea: Neocosmospora boniensis. Curr Genet. 2004;46:29–36.

10. Ma L-J, van der Does HC, Borkovich KA, Coleman JJ, Daboussi M-JM-J, Di Pietro A, et al. Comparative genomics reveals mobile pathogenicity chromosomes in Fusarium. Nature. 2010;464:367–73.

11. Balesdent M-H, Fudal I, Ollivier B, Bally P, Grandaubert J, Eber F, et al. The dispensable chromosome of *Leptosphaeria maculans* shelters an effector gene conferring avirulence towards *Brassica rapa*. New Phytol. 2013;198:887–98.

12. van der Does HC, Rep M. Adaptation to the host environment by plant-pathogenic Fungi. Annu Rev Phytopathol. 2017;55:427–50.

13. Lorrain C, Petre B, Duplessis S. Show me the way: rust effector targets in heterologous plant systems. Curr Opin Microbiol. 2018;46:19–25. https://doi.org/10.1016/j.mib.2018.01.016.

14. Toruño TY, Stergiopoulos I, Coaker G. Plant-pathogen effectors: cellular probes interfering with plant defenses in spatial and temporal manners. Annu Rev Phytopathol. 2016;54:419–41.

15. Zhao Z, Liu H, Wang C, Xu J-R. Comparative analysis of fungal genomes reveals different plant cell wall degrading capacity in fungi. BMC Genomics. 2013;14:274.

16. Lo Presti L, Lanver D, Schweizer G, Tanaka S, Liang L, Tollot M, et al. Fungal effectors and plant susceptibility. Annu Rev Plant Biol. 2015;66:513–45. https://doi.org/10.1146/annurev-arplant-043014-114623.

17. Snelders NC, Kettles GJ, Rudd JJ, Thomma BPHJ. Plant pathogen effector proteins as manipulators of host microbiomes? Mol Plant Pathol. 2018;19:257–9. https://doi.org/10.1111/mpp.12628.

18. Shi-Kunne X, Jové R de P, Depotter JRL, Ebert MK, Seidl MF, Thomma BPHJ. In silico prediction and characterisation of secondary metabolite clusters in the plant pathogenic fungus Verticillium dahliae. FEMS Microbiol Lett 2019; 366. doi:https://doi.org/10.1093/femsle/fnz081.

19. Sperschneider J, Dodds PN, Gardiner DM, Singh KB, Taylor JM. Improved prediction of fungal effector proteins from secretomes with EffectorP 2.0. Mol Plant Pathol. 2018;19:2094–110.

20. Stukenbrock EH, Banke S, Javan-Nikkhah M, McDonald BA. Origin and domestication of the fungal wheat pathogen Mycosphaerella graminicola via sympatric speciation. Mol Biol Evol. 2007;24:398–411.

21. Stukenbrock EH, Quaedvlieg W, Javan-Nikhah M, Zala M, Crous PW, McDonald BA. Zymoseptoria ardabiliae and Z. pseudotritici, two progenitor species of the septoria tritici leaf blotch fungus Z. tritici (synonym: *Mycosphaerella graminicola*). Mycologia. 2012;104:1397–407.

22. Grandaubert J, Bhattacharyya A, Stukenbrock EH. RNA-seq-based gene annotation and comparative genomics of four fungal grass pathogens in the genus Zymoseptoria identify novel orphan genes and species-specific invasions of transposable elements. G3 Genes, Genomes Genet. 2015;5:1323–33.

23. Goodwin SB, Ben MBarek S, Dhillon B, Wittenberg AHJ, Crane CF, Hane JK, et al. Finished genome of the fungal wheat pathogen Mycosphaerella graminicola reveals Dispensable structure, chromosome plasticity, and stealth pathogenesis. PLoS Genet. 2017;17:e100270.

24. Habig M, Quade J, Stukenbrock EH. Forward genetics approach reveals host genotype-dependent importance of accessory chromosomes in the fungal wheat pathogen Zymoseptoria tritici. MBio. 2017;8:e01919–7.
25. Habig M, Kema G, Holtgrewe Stukenbrock E. Meiotic drive of female-inherited supernumerary chromosomes in a pathogenic fungus. *Elife*. 2018;7:e40251.

26. Möller M, Habig M, Freitag M, Stukenbrock EH. Extraordinary genome instability and widespread chromosome rearrangements during vegetative growth. *Genetics*. 2018;210:517–29.

27. Fouché S, Plissonneau C, McDonald BA, Croll D. Meiosis leads to pervasive copy-number variation and distorted inheritance of accessory chromosomes of the wheat pathogen Zymoseptoria tritici. *Genome Biol* Evol. 2018;10:1416–29.

28. Plissonneau C, Hartmann FE, Croll D. Pangenome analyses of the wheat pathogen Zymoseptoria tritici reveal the structural basis of a highly plastic eukaryotic genome. *BMC Biol*. 2018;16:5.

29. Hartmann FE, McDonald BA, Croll D. Genome-wide evidence for divergent selection between populations of a major agricultural pathogen. *Mol Ecol*. 2018;27:2725–41. https://doi.org/10.1111/mec.14711.

30. Haueisen J, Möller M, Eschenbrenner CJ, Grandaubert J, Seybold H, Adamiak H, et al. Highly flexible infection programs in a specialized wheat pathogen. *Ecol Evol*. 2019;9:275–94.

31. Kellner R, Bhattacharyya A, Poppe S, Hsu TY, Brem RB, Stukenbrock EH. Expression profiling of the wheat pathogen Zymoseptoria tritici reveals genomic patterns of transcription and host-specific regulatory programs. *Genome Biol Evol*. 2014;6:1353–65. https://doi.org/10.1093/gbe/evu101.

32. Möller M, Schotanus K, Soyer JL, Haueisen J, Happ K, Stralucke M, et al. Destabilization of chromosome structure by histone H3 lysine 27 methylation. *PLoS Genet*. 2019;15:e1008093.

33. Badet T, Oggenfuss U, Abraham L, McDonald BA, Croll D. A 19-isolate reference-quality global pangenome for the fungal wheat pathogen Zymoseptoria tritici. *BMC Biol*. 2020;18:12. https://doi.org/10.1186/s12915-020-0744-3.

34. Plissonneau C, Stürchler A, Croll D. The evolution of orphan regions in genomes of a fungal pathogen of wheat. *MBio*. 2016;7:e01231–16. https://doi.org/10.1128/mBio.01231-16.

35. Haubold B, Klötzl F, Pfaffelhuber P and: Fast and accurate estimation of evolutionary distances between closely related genomes Bioinformatics 2015;31:1169–1175.

36. Zuo G, Hao B. CVTree3 web server for whole-genome-based and alignment-free prokaryotic phylogeny and taxonomy. *Genomics Proteomics Bioinformatics*. 2015;13:321–31.

37. Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, Kriventseva EV, Zdobnov EM. BUSCO Applications from Quality Assessments to Gene Prediction and Phylogenomics. *Molecular Biology and Evolution*. 2018;35(3):543–48.

38. Vaghefi N, Kikkert JR, Bolton MD, Hanson LE, Secor GA, Pethybridge SJ. De novo genome assembly of Circospora beticola for microsatellite marker development and validation. *Fungal Ecol*. 2017;26:125–34. https://doi.org/10.1016/j.funeco.2017.01.006.

39. Lechner M, Findeiß S, Steiner L, Marz M, Stadler PF, Prohaska SJ. Proteinortho: detection of (co-)orthologs in large-scale analysis. *BMC Bioinformatics*. 2011;12:124.

40. Lorrain C, Hecker A, Duplessis S. Effector-Mining in the Poplar Rust Fungus Melampsora larici-Populina Secretome. *Front Plant Sci*. 2015;6:1051.

41. Feurtey A, Stevens DM, Stephan W, Stukenbrock EH. Interspecific gene exchange introduces high genetic variability in crop pathogen. *Genome Biol Evol*. 2019;11:3095–105.

42. Wu B, Macielog AI, Hao W. Origin and spread of Spliceosomal introns: insights from the fungal clade Zymoseptoria. *Genome Biol Evol*. 2017;9:2658–67. https://doi.org/10.1093/gbe/evx211.

43. Valero-Jiménez CA, Veloso J, Staats M, van Kan JAL. Comparative genomics of plant pathogenic Botrytis species with distinct host specificity. *BMC Genomics*. 2019;20:203. https://doi.org/10.1186/s12864-019-5580-x.

44. Beckerson WC, Rodríguez De La Vega RC, Hartmann FE, Duhamel M, Giraud T, Perlin MH. Cause and effectors: Whole-genome comparisons reveal shared but rapidly evolving effector sets among host-specific plant-castrating fungi. *MBio*. 2019;10(6):e02391–19.

45. Thines M. An evolutionary framework for host shifts – jumping ships for survival. *New Phytol*. 2019;224:605–17. https://doi.org/10.1111/nph.16092.

46. Allen GC, Flores-Vergara MA, Krasynanski S, Kumar S, Thompson WF. A modified protocol for rapid DNA isolation from plant tissues using cetyltrimethylammonium bromide. *Nat Protoc*. 2006;1:2320–5.

47. Gurevich A, Saveliev V, Vyahii N, Tesler G. QUAST: quality assessment tool for genome assemblies. *Bioinformatics*. 2013;29:1072–5.

48. Fulnecková J, Sevcíková T, Fajkus J, Lukesová A, Lukes M, Vlcek C, et al. A broad phylogenetic survey unveils the diversity and evolution of telomeres in eukaryotes. *Genome Biol Evol*. 2013;5:468–83.

49. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. *Bioinformatics*. 2010;26:841–2.

50. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. *Nat Methods*. 2015;12:357–60.

51. Flutte T, Duprat E, Feuillet C, Quesneville H. Considering transposable element diversification in De novo annotation approaches. *PLoS One*. 2011;6:e16526.

52. Quesneville H, Bergman CM, Andrieu O, Autard D, Nouaud D, Ashburner M, et al. Combined evidence annotation of transposable elements in genome sequences. *PLoS Comput Biol*. 2005;1:e22.

53. Ter-Hovhannisyan V, Lomsadze A, Chernoff YQ, Borodovsky M. Gene prediction in novel fungal genomes using an ab initio algorithm with unsupervised training. *Genome Res*. 2008;18:1979–90.

54. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics*. 2014;30:2114–20.

55. Hoff KJ, Lange S, Lomsadze A, Borodovsky M, Stanke M. BRAKER1: unsupervised RNA-Seq-based genome annotation with GeneMark-ET and AUGUSTUS: table 1. *Bioinformatics*. 2016;32:767–9.

56. Stanke M, Diekhans M, Baertsch R, Haussler D. Using native and syntenically mapped cDNA alignments to improve de novo gene finding. *Bioinformatics*. 2008;24:637–44.

57. Lomsadze A, Burns PD, Borodovsky M. Integration of mapped RNA-Seq reads into automatic training of eukaryotic gene finding algorithm. *Nucleic Acids Res*. 2014;42:e119.

58. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. *Nat Biotechnol*. 2011;29:644–52.

59. Haas BJ, Salzberg SL, Zhu W, Pertea M, Allen JE, Orvis J, et al. Automated eukaryotic gene structure annotation using EvidenceModeler and the program to assemble spliced alignments. *Genome Biol*. 2008;9:R7.

60. Gremme G, Steinbiss S, Kurtz S. GenomeTools: a comprehensive software library for efficient processing of structured genome annotations. *IEEE/ACM Trans Comput Biol Bioinforma*. 2013;10:645–56.

61. Lechner M, Hernandez-Rosales M, Doerr D, Wieseke N, Thévenin A, Stoye J, et al. Orthology detection combining clustering and synteny for very large datasets. *PLoS One*. 2014;9:e105015.

62. Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, et al. Circos: an information aesthetic for comparative genomics. *Genome Res*. 2009;19:1639–45. https://doi.org/10.1101/gr.092759.109.

63. Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, von Mering C, et al. Fast genome-wide functional annotation through Orthology assignment by eggNOG-mapper. *Mol Biol Evol*. 2017;34:2115–22.

64. Zhang H, Yohe T, Huang L, Entwistle S, Wu P, Yang Z, et al. dbCAN2: a meta server for automated carbohydrate-active enzyme annotation. *Nucleic Acids Res*. 2018;46:W95–101.

65. Weber T, Blin K, Duddela S, Krug D, Kim HU, Bruccoleri R, et al. antiSMASH 3.0—a comprehensive resource for the genome mining of biosynthetic gene clusters. *Nucleic Acids Res*. 2015;43:W237–43.

66. Almagro Armenteros JJ, Sønderby CK, Sønderby SK, Nielsen H, Winther O. DeepLoc: prediction of protein subcellular localization using deep learning. *Bioinformatics*. 2017;33:3387–95.

67. Petersen TN, Brunak S, von Heijne G, Nielsen H. SignalP 4.0: discriminating signal peptides from transmembrane regions. *Nat Methods*. 2011;8:785–6.

68. Emanuelsson O, Nielsen H, Brunak S, von Heijne G. Predicting subcellular localization of proteins based on their N-terminal amino acid sequence. *J Mol Biol*. 2000;300:1005–16.

69. Käll L, Krogh A, Sonnhammer EL. A combined Transmembrane topology and signal peptide prediction method. *J Mol Biol*. 2004;338:1027–36.

70. Krogh A, Larsson B, von Heijne G, Sonnhammer EL. Predicting transmembrane protein topology with a hidden markov model: application to complete genomes. *J Mol Biol*. 2001;305:567–80.

71. Sonnhammer EL, von Heijne G, Krogh A. A hidden Markov model for predicting transmembrane helices in protein sequences. *Proc Int Conf Intell Syst Mol Biol*. 1998;6:175–82.

72. Sperschneider J, Dodds PN, Singh KB, Taylor JM. ApoplastP: prediction of effectors and plant proteins in the apoplast using machine learning. *New Phytol*. 2018;217:1764–78.
73. Sperschneider J, Catanzariti A-M, DeBoer K, Petre B, Gardiner DM, Singh KB, et al. LOCALIZER: subcellular localization prediction of both plant and effector proteins in the plant cell. Sci Rep. 2017;7:44598.

74. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

75. Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, et al. Software for computing and annotating genomic ranges. PLoS Comput Biol. 2013;9:e1003118.

76. Soyer JL, Möller M, Schotanus K, Connolly LR, Galazka JM, Freitag M, et al. Chromatin analyses of Zymoseptoria tritici: methods for chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq). Fungal Genet Biol. 2015;79:63–70.

77. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357–9.

78. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38:576–89.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions