Nucleosome patterns in four plant pathogenic fungi with contrasted genome structures

Fungal pathogens represent a serious threat towards agriculture, health, and environment. Control of fungal diseases on crops necessitates a global understanding of fungal pathogenicity determinants and their expression during infection. Genomes of phytopathogenic fungi are often compartmentalized: the core genome contains housekeeping genes whereas the fast-evolving genome mainly contains transposable elements and species-specific genes. In this study, we analysed nucleosome landscapes of four phytopathogenic fungi with contrasted genome organizations to describe and compare nucleosome repartition patterns in relation with genome structure and gene expression level. We combined MNase-seq and RNA-seq analyses to concomitantly map nucleosome-rich and transcriptionally active regions during fungal growth in axenic culture; we developed the MNase-seq Tool Suite (MSTS) to analyse and visualise data obtained from MNase-seq experiments in combination with other genomic data and notably RNA-seq expression data. We observed different characteristics of nucleosome profiles between species, as well as between genomic regions within the same species. We further linked nucleosome repartition and gene expression. Our findings support that nucleosome positioning and occupancies are subjected to evolution, in relation with underlying genome sequence modifications. Understanding genomic organization and its role in expression regulation is the next gear to understand complex cellular mechanisms and their evolution.


Introduction
Fungi account for a huge part of the Earth biodiversity with a current estimate of 2.2 to 3.8 million species (Hawksworth and Lücking, 2017). Fungi are organisms of major environmental importance as they develop beneficial symbiotic associations with plants and are able to decay dead organic matter (Taylor and Osborn, 1996). Unfortunately, fungi are also very efficient pathogens causing important damages in agriculture, human health, and the environment (Fisher et al., 2012). Control of fungal diseases on crops necessitates a global understanding of fungal pathogenicity determinants and the control of their expression during infection. Among these pathogenicity determinants, fungi secrete an arsenal of molecules known as effectors, key elements of pathogenesis which modulate innate immunity of the plant and facilitate infection. Effectors can be small proteins, secondary metabolites and small RNAs (Weiberg et al., 2013;Lo Presti et al., 2015;Collemare, O'Connell and Lebrun, 2019). Upon plant infection, fungi undergo a tightly controlled transcriptional reprogramming, and different sets of effectors are expressed at specific stages of pathogen development and host colonization (Toruño, Stergiopoulos and Coaker, 2016;Zhao et al., 2021;Ding, Gardiner and Kazan, 2022). Plant-associated fungi generally show contrasted genomic landscapes including 'plastic' loci with a high prevalence of transposable elements (TE). These genomes either show an overall large proportion of TE evenly distributed throughout the genome, or TE clustered in specific regions such as long TE-rich blocks, accessory chromosomes or subtelomeric areas (Sánchez-Vallet et al., 2018). Effector genes are over-represented in these TE-rich regions. TE-rich compartments have heterochromatin properties contrary to TE-poor regions which have euchromatin properties. The location of effector genes in regions enriched in TEs has been shown to provide a tight control of their expression through chromatin remodeling. Indeed, several recent studies pointed out the potential role of chromatin remodeling in the regulation of effector-encoding genes and the control of secondary metabolism (reviewed in Soyer, Rouxel and Fudal, 2015;Collemare and Seidl, 2019).
Eukaryotic chromatin is packaged into nucleosomes, each composed of DNA wrapped around a histone octamer associated with various other proteins, and separated by linker DNA (Richmond and Davey, 2003). These histone proteins are composed of histone core where the DNA is wrapped and histone tails which can be chemically modified by specific enzymes changing the chromatin 3D-structure and DNA accessibility to polymerases and transcription factors (TF). Nucleosome assembly is further stabilized by the binding of a linker histone H1. Positioning of nucleosomes throughout the genome and post-translational modifications of histones have a significant regulatory function by modifying availability of binding sites to TF and to polymerases, affecting DNA-dependent processes such as transcription, DNA repair, replication and recombination (Radman-Livaja and Rando, 2010;Struhl and Segal, 2013). Nucleosome positioning (i.e., the position of the nucleosome along the DNA sequence) and occupancy (i.e., a measure of the actual level of occupation of a given position by a nucleosome in a pool of cells) are determined by a combination of DNA sequence features, TF, chromatin remodelers and histones modifiers (see (Singh and Mueller-Planitz, 2021) for a review). Genome-wide maps of nucleosome occupancy and positioning are still sparse in fungi and have only been developed in a few Hemiascomycota yeast species, including Saccharomyces cerevisiae (Yuan et al., 2005;Tsankov et al., 2010), in the ascomycete Aspergillus fumigatus (Nishida et al., 2009) and the basidiomycete Mixia osmundae (Nishida et al., 2012). The studies revealed that promoter, enhancer and terminator regions were depleted in nucleosomes, allowing access to TF, and that the nucleosomal DNA length distribution was similar in M. osmundae and A. fumigatus but differed from that of hemiascomycetous yeasts. No comparative genome-wide analyses of nucleosome positioning have been performed in ascomycetes and notably not in plant pathogenic fungi.
In the present study, we investigated genome-wide nucleosome localization in four different plant pathogenic ascomycetes showing different genomic organizations: i) Leptosphaeria maculans 'brassicae' (Lmb), a hemibiotrophic pathogen of Brassica species, including oilseed rape; ii) the most closely related species of Lmb, Leptosphaeria maculans 'lepidii' (Lml), a pathogen of Lepidium spp.; iii) Fusarium graminearum, a hemibiotrophic pathogen of cereals and iv) Botrytis cinerea, a polyphagous necrotrophic pathogen causing grey mould on more than 1,400 plant species. The genome of Lmb has been invaded by TE (which represent more than 30 % of its genome) and is composed of alternating compartments: generich GC-equilibrated and TE-rich AT-rich genomic regions (Rouxel et al., 2011;Dutreux et al., 2018). In contrast, the Lml genome presents only 4 % of repeats which are evenly distributed throughout the genome (Grandaubert et al., 2014). Genomes of F. graminearum and B. cinerea have a very low to low TEcontent (Cuomo et al., 2007;Amselem et al., 2011;King et al., 2015;King, Urban and Hammond-Kosack, 2017). The genome of the reference strain of B. cinerea, B05.10, contains 4 % of TE, which are localized essentially in the telomeric and centromeric regions of the core chromosomes, or on the two dispensable chromosomes (Amselem et al., 2011;Porquier et al., 2016), i.e., chromosomes not essential to immediate survival and missing in some or most individuals (Soyer et al., 2018). The genome of F. graminearum contains very little TE identified to date (0.3 % (Cuomo et al., 2007;King et al., 2015;King, Urban and Hammond-Kosack, 2017)).
In this study, we compare nucleosome repartition patterns in relation with genome structure and gene expression level in these four phytopathogenic Ascomycota. To gain insight into the role of nucleosome positioning and occupancy in regulating fungal pathogen transcription, we applied micrococcal nuclease digestion of mono-nucleosomes couple with high-throughput sequencing (MAINE-seq or MNase-seq) with regards to mRNA abundance to concomitantly map nucleosome-rich regions and transcriptionally active regions during fungal growth.

Strains and culture conditions
The studied fungi were cultured independently in the media and conditions classically used for each of them. Leptosphaeria maculans 'brassicae' v23.1.3 and Leptosphaeria maculans 'lepidii' IBCN84 mycelia were inoculated into 100 mL of Fries liquid medium (1 g/L NH4NO3, 5 g/L C4H12N2O6, 1 g/L KH2PO4, 0.5 mg/L MgSO4 7H2O, 130 mg/L CaCl2, 100 mg/L NaCl, 30 g/L C12H22O11 and 5 g/L Yeast extract). Tissues were harvested after growing for seven days in the dark at 25°C. Botrytis cinerea strain B05.10 (10 6 spores/mL) was grown for two days on solid Malt Medium (MM, 20 g/L malt extract, 5 g/L yeast extract and 15 g/L agar) covered with a cellophane layer (Simon et al., 2013;Kelloniemi et al., 2015). The plates were incubated in a growth chamber (Sanyo MLR-350H) at 23°C with an alternation of 14 h of white light and 10 h of darkness. After two days of culture, mycelia were ground in liquid nitrogen and stored to -80°C until further processing. Fusarium graminearum strain CBS185.32 (Centraal Bureau voor Schimmelcultures, Utrecht, the Netherlands) was grown for three days in modified liquid MS (glucose was substituted with sucrose) as previously described (Boutigny et al., 2009). All cultures were done in three biological replicates.

Preparation of nucleosomal DNA
Fungal material was harvested and treated with microccocal nuclease (MNase, cat. #MS0247S, New England BioLabs). For Lmb and Lml, ∼300 mg of mycelium were digested with 5 µL of MNase for 10 min at 37°C , directly followed by DNA purification as previously described (Soyer et al., 2020). For F. graminearum, mycelia were harvested by filtering and immediately homogenized for 1 min at 30 Hz using a TissueLyzer (Qiagen). Then, 100 mg of ground mycelium was digested for 10 min at 37°C with 15 µL of MNase in 600 µL of digestion buffer (0.6% v/v IGEPAL, 50 mM NaCl, 2 mM Pefabloc, 50 mM Tris-HCl pH8, 10 mM CaCl2). The reactions were stopped with 10mM EDTA and the samples treated with RNAse followed by proteinase K prior DNA purification with phenol/chloroform and ethanol precipitation. For B. cinerea, we digested 100-200 mg of mycelium per sample with 1 µL MNase at 37°C  for 10 min. The reactions were stopped by adding 10mM EDTA and samples treated with RNAse A followed by proteinase K. DNA purification was realized with the "Nucleospin Gel and PCR clean up kit" (Macherey Nagel, cat #740609.250). For all samples, nucleic acid quantification was performed by UV spectrometry using a Nanodrop-ND 1000 apparatus, and digestion profiles were checked by 2% agarose gel electrophoresis. Nucleosomal DNA was stored at -20°C until DNA library preparation.

Extraction of total RNA
For Lmb and Lml, total RNA was extracted from mycelium grown for one week in Fries liquid medium as previously described (Fudal et al., 2007). For F. graminearum, mycelia were harvested by filtering, rinsed twice with sterile deionized water, and flash frozen in liquid nitrogen. One milliliter of TRIzol TM Reagent (Thermo Fischer Scientific) was added to 200 mg of mycelium before grinding for 1.5 min at 30 Hz using a TissueLyzer (Qiagen). Total RNA was then extracted using a previously published protocol (Hallen et al., 2007). For B. cinerea, total RNA was extracted from frozen ground mycelium using a previously published protocol (Kelloniemi et al., 2015). All total RNA samples were stored at -80°C until preparation of RNA library.

Preparation of sequencing libraries, high-throughput sequencing, and read pre-processing
MNase-seq libraries were prepared from purified nucleosomal DNA using the kit NEBNext Ultra DNA Library Prep Kit for Illumina (cat. # E7370L New England BioLabs) following the manufacturer's instructions. The NEBNext Ultra Directional RNA Library Prep Kit for Illumina (cat. # E7420L New England BioLabs) was used to prepare all RNA-seq libraries, following the manufacturer's instructions. Sequencing was performed by the GenomEast platform, a member of the 'France Génomique' consortium (ANR-10-INBS-0009). Samples were run in 9-plex on an Illumina HiSeq 4000 in paired mode, 2x50 bp. Initial read quality check was performed using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Raw reads were then pre-processed with Trimmomatic v0.32 (Bolger, Lohse and Usadel, 2014) to clip out any remaining sequencing adapter sequence and crop low quality nucleotides (minimum accepted Phred score of 30). Reads in pairs of 40 bp or more in length were used in the present analysis.

Transcriptome analyses
RNA-seq reads were mapped against their respective reference genomes (see Table 1) using STAR v2.5.1 (Dobin et al., 2013). TPM counts (Transcripts Per Million reads, (Li et al., 2010)) were computed using the count TPM tool provided with the MNase-Seq Tool Suite (MSTS; Supplementary Figure 1) that was developed in-house to analyse genome-wide nucleosome positioning data combined with RNA-seq data (Lapalu 2023).

MNase-seq analyses
MNase-seq paired-end reads were mapped using Bowtie2 software ran in very-sensitive mode (Langmead and Salzberg, 2012). MSTS (MNase-seq Tool Suite) was used to compute all phasograms, dinucleotide composition, as well as nucleosome density profiles of genomic compartments and/or gene lists (Supplementary Figure 1). Lists of near-universal single copy orthologs were obtained by running BUSCO3 for Fungi (Simão et al., 2015;Waterhouse et al., 2018) on each reference genome studied. Graphical visualisations were computed with MSTS and Matlab R2020b (MathWorks). Frequency distributions of read coverage per base, obtained with the Phasogram function of MSTS were scaled (zscore) and plotted with Matlab. For each replicate, phases, standard errors (se), R 2 (coefficient of determination), and p-values (F-test) were determined after linear regression fitting to the first four successive peak positions.

Genome-wide nucleosome spacing
We first explored nucleosome landscapes in the four fungal genomes by measuring the average distance between nucleosomes genome-wide; we computed phasograms, i.e., frequency distributions of coverage per base genome-wide for all four species (Figure 1 and Supplementary Table 2). Phasograms obtained in nucleosome mapping resemble oscillating sine wave signals, for which period is the length of DNA bound to the histone octamer plus the length of the DNA stretch to the next nucleosome, averaged genome wide (Valouev et al., 2011). Phasing signals were observed genome-wide over 1,200 bp sliding windows revealing six to seven nucleosome peaks in a wave signal decaying in intensity with increasing distance and significant linear regression on peak apex positions, as previously described (Valouev et al., 2011). We found that, in the fungi studied, nucleosomes are 161 to 172 bp distant from each other (centre to centre), also called nucleosome repeat length or NRL (i.e., the length of DNA wrapped around the histone octamer plus linker DNA), depending on the considered species and culture condition. In B. cinerea, average NRL is estimated at 169 bp ( Figure 1A). In F. graminearum, this distance reaches 172 bp ( Figure  1B). In Lmb and Lml ( Figures 1C and 1D), average NRL is 166 bp and 161 bp, respectively. Considering these values and the canonical length of nucleosomal DNA (147 bp), linker DNA length can be estimated to stand, in average, between 14 to 19 bp for respectively Lml and Lmb, 22 bp for B. cinerea, and 25 bp for F. graminearum. Nucleosome phasing genome-wide seems to be particularly tight in F. graminearum, with very little deviation in the measured phases ( Figure 1B and Supplementary Table 2). In contrast, higher deviations are observed for B. cinerea ( Figure 1A and Supplementary Table 2).
Colin Clairet et al. Figure 1: Nucleosome phasing in the four fungi studied. Main graphs display scaled (z-score) phase frequencies (y-axis) as a function of position (in base pair; x-axis). Graphs in inserts show peak positions (in base pairs; y-axis) as a function of peak order (x-axis). For each replicate, phases, standard errors (se), R 2 (coefficient of determination), and p-values (F-test) are determined after linear regression fitting to the first four successive peak positions (see Supplementary Table 2). Repl = replicate. A. Botrytis cinerea; B. Fusarium graminearum; here, phase value in replicate #2 was measured for four successive peaks excluding peak #1 for which an apex was not clearly visible at the beginning of the profile; C. Leptosphaeria maculans 'brassicae'; D. Leptosphaeria maculans 'lepidii'.
Nucleosome spacing influences the formation of the higher order chromatin fibre, often referred to as the 30-nm chromatin fibre (Szerlong and Hansen, 2011). Several structural models of the chromatin fibre have been proposed, all underlining nucleosome-nucleosome interactions including the length of linker DNA fragments as major driving factors (Zhu and Li, 2016). Notably, the chromatin fibre was found to be narrower (21-nm in diameter) for a short NRL of 167 bp (Robinson et al., 2006). Similarly, increasing NRLs lead to increasingly wider fibres, reaching a highly compact 30-nm solenoid structure for an NRL of 197 bp. Typically, NRLs are ~175-200 bp in plants, Caenorhabditis elegans, and humans (Valouev et al., 2008(Valouev et al., , 2011Locke et al., 2013;Zhang, Zhang and Jiang, 2015), and ~165 bp and 154 bp in the yeasts S. cerevisiae and Schizosaccharomyces pombe, respectively (Yuan et al., 2005;Lantermann et al., 2010). In the present study, we found NRL values remarkably constant between biological replicates, a phenomenon sometimes referred to as clamping that involves ATP-dependent chromatin remodelers in purified experimental systems (Lieleg et al., 2015). The species B. cinerea and F. graminearum show similar NRLs in the middle range (169 bp to 172 bp), possibly indicating intermediate levels of compaction of the higher order chromatin fibre. These values are similar to those obtained for the Pezizomycotina A. fumigatus, i.e., linker length ranging from 21 to 27 bp, using an MNase treatment similar to the one used in the present study (Nishida et al., 2009). Lml and Lmb distinguish themselves with shorter NRLs of only 161-166 bp, suggesting a narrower chromatin fibre structure. In the yeast S. cerevisiae, linker length was shown to be the result of the competition for binding between the chromatin remodeling factors ISW1a (Imitation SWItch) and CHD1 (Chromodomain Helicase DNA-binding), the latter mediating shorter length after the eviction of histone H1 (Ocampo et al., 2016). Considering that sequence polymorphisms in CHD1 has been previously associated with variations of linker length (Hughes and Rando, 2015), a seducing possibility is that such variation at the protein level may account for a portion for the inter-species differences observed here. Regarding the very closely related species Lmb and Lml, Lmb presents longer NRLs than Lml. We hypothesized this peculiarity may be explained by large AT-rich regions displayed by the Lmb genome, not encountered in the Lml genome (Rouxel et al., 2011;Grandaubert et al., 2014). Indeed, DNA sequence is a major determinant of nucleosome landscapes (Struhl and Segal, 2013), in particular AT stretches that confer more rigidity to the chromatin fiber.

Nucleosome distribution profiles
Read density was plotted genome-wide in one kb-long sliding non-overlapping windows along chromosomes for all four fungi (Figures 2 to 4, and Supplementary Figures 2 to 6). The density profiles obtained for B. cinerea show remarkable regularity of nucleosome density genome-wide (Figure 2A and Supplementary Figure 2). Nevertheless, we could observe that almost all occasional thin peaks of density were correlated with the positions of BOTY retro-transposons (Porquier et al., 2016). Out of 48 complete copies of BOTY in the genome, 31 show a peak of nucleosome density. Notably, they correspond to the BOTY elements with an equilibrated percentage of GC (43-45%) while the 17 copies that do not show such a peak are those with a lower percentage of GC (14-24%) probably because they have undergone Repeat-Induced Point mutation, or RIP (Porquier et al., 2016). Peaks of density were rarely observed for TE other than BOTY (Supplementary Figure 3). We also investigated nucleosome spacing in regions occupied by BOTY and non-BOTY TE. Phasograms were plotted as described above restricting our analysis to BOTY or other TE ( Figures 2B and 2C, respectively). Much larger phases can be observed in other TE regions (178.4-187.5 bp) when compared to BOTY regions (171.3-172.5 bp) or genome-wide (168.2-169.4 bp, Figure 1A), indicating larger nucleosome spacing in TE-occupied regions. BOTY-containing regions, which positions correlate with discrete peaks of nucleosome density, exhibit slightly larger phases than genome-wide. Thus, the observed peaks of read density may be the result of increased nucleosome occupancy, i.e., a measure of the stability of a nucleosome at a given position in a multiple cell sample, rather than a denser deposition of nucleosomes. BOTY is one of the largest TE identified in the B05.10 strain (6.4-6.6 kb), and that's definitively the TE with the largest genome coverage i.e., 0.96% (Porquier et al., 2016). Notably, the majority of B. cinerea small interfering RNA (siRNA) predicted to silence host plant genes are derived from the copies of BOTY and related elements that show an equilibrated percentage in GC (Weiberg et al., 2013;Porquier et al., 2021). As the production of these siRNA effectors is activated during the early phase of plant infection, we could speculate that the high nucleosome occupancy on the loci of production (i.e. un-RIPped BOTY TEs) is a mechanism to restrict their production during saprophytic growth. The observation of two distinct chromatin states characterizing TEs in Verticillium dahliae supports this proposition (Cook et al., 2020). This hypothesis remains to be tested by investigating nucleosome occupancy during in planta development.
In F. graminearum, regions equally packed with nucleosomes are interspaced with areas with lower density ( Figure 3A and Supplementary Figure 4). Strikingly, this profile mirrors the previously described SNP density profiles in F. graminearum (Laurent et al., 2017). We investigated whether or not this profile was the result of increased spacing between nucleosomes in regions found denser in SNPs. Phasograms were plotted restricting our analysis to the SNP-enriched polymorphic islands or the rest of the genome, as defined by Laurent et al. (Laurent et al., 2017). Wave signals similar to the ones observed genome-wide were obtained, with phases of 172.3-172.4 bp in polymorphic islands ( Figure 3B, Supplementary Table 3) and 171.6-171.9 bp outside these regions ( Figure 3C, Supplementary Table 3). These results indicate that nucleosomes appear well-arrayed genome-wide, with very similar phases in polymorphic islands vs. nonpolymorphic islands. Thus, the observed drops in read density profile cannot be explained by a depletion in nucleosomes but may rather be the result of reduced nucleosome occupancy. This observation suggests increased frequencies of transient nucleosome positioning events in F. graminearum fast evolving polymorphic islands (Laurent et al., 2017) and thus more relaxed chromatin structures. Here, nucleosome dynamics may enable fast evolution of particular genome segments while regions defined by higher occupancies may secure sequence conservation.  Figure 6), suggesting that AT-rich regions are particularly dense in nucleosomes. Considering the remarkable compartmentalized organization of the genome of Lmb (Rouxel et al., 2011), absent from Lml, differences of nucleosome phasing and occupancy in TE-and ATrich vs. GC-equilibrated and gene-rich regions were investigated. A region was considered AT-rich if it contained less than 40 % of GC. As described earlier, AT-rich regions represent one-third of the Lmb genome divided in 419 regions of 1 to 320 kb in length. Examination of unprocessed mapping outputs reveals that the number of fragments (read pairs) mapped in AT-rich and GC-equilibrated regions were very similar, with 23.8 million and 24.8 million fragments, respectively, which is far from the 1/3 vs. 2/3 ratio expected. In terms of coverage depth, mean coverage is 207 vs. 135 fragments for AT-rich and GCequilibrated regions, respectively, which could suggest higher nucleosome occupancy in the former. We explored this hypothesis and compared phasograms for AT-rich vs. GC-equilibrated regions ( Figure 4B and 4C, Supplementary Table 3). Average NRLs were found larger in AT-rich than GC-equilibrated compartments, measured at 169.2 bp and 164.2 bp respectively, suggesting lower nucleosome frequencies in AT-rich regions than in GC-equilibrated regions. Nonetheless, coverage density is higher in AT-rich regions ( Figure 4A and Supplementary Figure 5), consistently with our hypothesis of higher nucleosome occupancy in these regions and thus less accessible genome compartment, in heterochromatic state. This is in accordance with the recent genome-wide mapping of histone modifications performed by Soyer et al. (2020) on Lmb and Lml in which the Histone H3 Lysine9 tri-methylation heterochromatin mark was found associated with TE-and AT-rich regions of Lmb including on the dispensable chromosome (Soyer et al., 2020). Finally, signal intensity in phasograms appeared more stable on the long nucleotide range in GCequilibrated than in AT-rich regions, an observation in line with the well-known destabilizing effect of AT stretches on nucleosome positioning leading to fuzzier signals (Kaplan et al., 2009;Ponts et al., 2010;Tillo et al., 2010;Bunnik et al., 2014;Russell et al., 2014;Jin, Finnegan and Song, 2018). All together, these data support a heterochromatic state of Lmb AT-rich regions mediated by strong nucleosome occupancy during axenic growth. Since the AT-rich regions host many fungal effector genes expressed during primary infection of oilseed rape leaves, we may assume that these regions are decondensed during infection, allowing the action of specific transcription factors. We tried to perform MNase-seq experiments at an early stage of oilseed rape infection by Lmb but the number of fungal reads was too low to be able to reliably analyze fungal nucleosome positioning. To go further, techniques to enrich in fungal material prior to MNase treatment should be considered. A. Coverage density profiles were computed for nonoverlapping 1 kb-long bins along the four chromosomes of F. graminearum. In green are plotted SNP density profiles as previously described (Laurent et al., 2017). "nd" indicates the highly variable 3' end of chromosome 4 for which SNP were not called. In blue are plotted the z-scored average nucleosome density profile (see Supplementary Figure 2 for individual replicate plots). Black arrows indicate centromeres (King et al., 2015); B and C. Nucleosome phasing in polymorphic islands (B) or outside polymorphic islands (C) as previously defined (Laurent et al., 2017). Main graphs display phase frequencies (y-axis) as a function of position (in base pair; x-axis). Graphs in inserts show peak positions (in base pairs; y-axis) as a function of peak order (x-axis); Phases +/-standard errors (se), R 2 (coefficient of determination), and p-values (F-test) are determined after linear regression fitting to the first four successive peak positions (see Supplementary Table 3). Repl = replicate.

2011) (B) and GC-equilibrated regions (C). Main graphs display phase frequencies (y-axis) as a function of position (in base
pair; x-axis). Graphs in inserts show peak positions (in base pairs; y-axis) as a function of peak order (x-axis); for each replicate, phases, standard errors (se), R 2 (coefficient of determination), and p-values (F-test) are determined after linear regression fitting to the first four successive peak positions (see Supplementary Table 3). Repl = replicate.

Sequence composition and nucleosome positioning
Literature data report that distribution of bases in nucleosome core DNA is non-random and exhibits a ~10 bp AA/TT/AT/TA offset with GG/CC/GC/CG dinucleotide frequency (Satchwell, Drew and Travers, 1986;Segal et al., 2006;Segal and Widom, 2009b). Here, we investigated di-nucleotide frequencies -i.e., the incidence of a given neighbouring pair of nucleotides in a sequence -in nucleosomal DNA segments in all four fungi. Averaged di-nucleotide contents centred around all fragments (read pairs) were plotted ( Figure  5) and periodicities investigated by autocorrelation analyses (Supplementary Figure 7). Autocorrelation plots reveal the previously described ~10 bp-periodicities (Laurent et al., 2017) (Satchwell, Drew and Travers, 1986;Segal et al., 2006;Segal and Widom, 2009b) for all studied fungi while showing differences in instant autocorrelation coefficient profiles (Supplementary Table 4 and Supplementary Figure 7). Signal is indeed very regular in F. graminearum and, to a lesser extent, Lmb, whereas Lml and B. cinerea show more irregular autocorrelation profiles. These results are consistent with our previous observation that nucleosomes are tightly phased in F. graminearum whereas somehow fuzzier (higher deviation) in B. cinerea and Lml ( Figure 1A and 1D).
Di-nucleotide frequency graphs display A+T dinucleotides frequency waves oscillating out of phase with G+C ones. For all studied fungi, GC dinucleotides are centred on nucleosome dyads ( Figure 5). Considering that there are 16 possible combinations of di-nucleotides, equilibrated distribution of di-nucleotide contents should contain 25% of AT/TA/AA/TT and 25% of GG/CC/CG/GC. Our observations reveal a skewed distribution in favour of AT dimers marked for B. cinerea and Lmb ( Figure 5A and 5C. The presence of ATrich regions in Lmb and, to a lesser extent, in B. cinerea genomes may explain such a result (Rouxel et al., 2011;Porquier et al., 2016). Di-nucleotide frequencies of Lmb AT-rich and GC-equilibrated regions were thus inspected (Supplementary Table 5 and Supplementary Figure 8). As one would expect, A+T frequencies are particularly high in AT-rich regions, while maintaining alternance with G+C di-nucleotides and ~10 bp periodicity. Nucleosome positioning is believed to be particularly hard-wired to DNA sequence, and especially the largely documented anti-nucleosome effect of Poly(dA:dT) tracts (Segal and Widom, 2009a;Struhl and Segal, 2013). Here, while nucleosome phase was indeed found 5 bp longer in Lmb AT-rich regions than in GC-equilibrated regions, occupancy was nonetheless higher in the former leading to the formation of the previously suggested heterochromatic state of these regions, which has consequences on gene expression and recombination (Rouxel et al., 2011;Soyer et al., 2014Soyer et al., , 2020Gay et al., 2020). These observations suggest the mobilisation of trans-acting chromatin remodeling factors to maintain heterochromatin structures on such disfavouring sequences. Importantly, we found that GC periodicity at nucleosome dyads is preserved even within AT-rich regions, suggesting such pattern in an AT-rich environment is sufficient to permit efficient wrapping of DNA around nucleosomes and strong occupancy. This longer NRLs in addition with high nucleosome density in AT-rich regions should have an impact on gene expression in these regions enriched in effector gene (Rouxel et al., 2011;Soyer et al., 2020).

Nucleosome landscapes of fungal gene units
Nucleosome occupancy profiles of gene units were investigated in all fungi ( Figure 6). As previously reported in other eukaryotes, translation start sites (the ATG codon) are preceded by a nucleosomedepleted region (NDR) and immediately followed by a well-positioned +1 nucleosome ( Figure 6A). Variations of the exact position of these features relative to the ATG start codon, as well as variations in the intensity of the NDR valley, are nonetheless observed between fungal species. For example, the NDR and the centre of the +1 nucleosome (or +1Nucl) are found at -154 bp and +14 bp, respectively, in F. graminearum whereas they are located at -99 bp and + 26 bp, respectively, in B. cinerea. The fungus B. cinerea shows the deepest NDR valley among all observed profiles. In Lmb and Lml, NDRs are found at -129 and -144 bp from ATG, respectively, and +1Nucl at +19 bp and +55 bp. Finally, nucleosome profiles upstream of NDRs appear fuzzy for all fungi but Lmb, with varying degrees of fuzziness. This fuzziness is no longer visible when nucleosome profiles are centred on TSS for F. graminearum (NTSS = 6,212 genes) and B. cinerea (NTSS = 11,701 genes) ( Figure 6C). The NDR is more defined, and located immediately upstream of the TSS, with a minimum detected at -58 bp and -20 bp upstream of the TSS of F. graminearum and B. cinerea, respectively. These values are consistent with the binding of the RNA polymerase II ~50 bp upstream of TSS, observed in active promoters of mammalian and Drosophila cells (Core, Waterfall and Lis, 2008;Min et al., 2011;Kwak et al., 2013). Considering nucleosome environments at stop codons ( Figure 6B), strongly arrayed nucleosomes are particularly found upstream the stop codon, fewer signal variations being observed downstream. Here, the stop codon is a clear boundary for nucleosome arraying and occupancy in all fungi and all conditions investigated. A nucleosome seems remarkably well positioned exactly on stop codons in F. graminearum in particular. The analysis was repeated on TTS in F. graminearum (NTTS = 5,292 genes) and B. cinerea (NTTS = 11,701 genes) ( Figure 6D). In B. cinerea, signal appeared strong and well-arrayed, decaying downstream of the TTS. In F. graminearum, strong positioning of nucleosomes -2 and -1 (-178 bp and -16 bp relative to TTS, respectively) followed by a deep 3' end NDR (+100 bp downstream of TTS) can be observed. Nucleosome positioning on TSS and TTS could not be analyzed for Lmb and Lml since TSS and TTS annotations are not supported by experimental data for these species' genomes as it is the case for F. graminearum (King, Urban and Hammond-Kosack, 2017;Basenko et al., 2018), or by collaborative annotation as for B. cinerea (Pedro et al., 2019).
For comparison purposes, nucleosome profiling was repeated restricting our analysis to genes identified as BUSCO lineage-specific single-copy evolutionary conserved orthologs in Fungi (Simão et al., 2015;Waterhouse et al., 2018). Overall broad patterns remain similar to those obtained while investigating whole genomes, with notably somehow more regular oscillations patterns (Figure 7). In F. graminearum, whilst the distance NDR-to-TSS (-60 bp) remains very similar to the one measured earlier genome-wide ( Figure 7C), distance NDR-to-ATG increases by 29 bp (Figure 7A), whereas in B. cinerea the distance NDRto-TSS reduces to 0 bp and nearly no increase in the distance NDR-to-ATG is observed (+ 3 bp). Similarly, the distance NDR-to-ATG is only 5 bp longer than genome-wide in Lmb whereas it increases by 23 bp in Lml. Towards the 3' end of the gene unit, nucleosome signals around stop codons are similar to those obtained genome-wide for all fungi ( Figure 7B). However, the deep NDR found downstream of TTS of F. graminearum genome-wide can no longer be observed at "Fungi" TTS loci ( Figure 7D). Similar to genomewide profile, strongly positioned -1Nucl and -2Nucl are still visible at -5 bp and -209 bp, respectively, the latter being more intense and defined. The general profile of a fungal gene unit shares similarities with those previously described in various eukaryotes: the ATG codon is decorated by a well-positioned +1 nucleosome and preceded by an NDR. This "+1 nucleosome" is an extremely well-conserved feature among eukaryotes spread across the tree of life, and nucleotide sequence only is not sufficient to explain such consistency. Such stability requires the intervention of ATP-dependent chromatin remodelers, belonging to one of the families CHD, INO80, ISWI, or SWI/SNF (Reyes, Marcum and He, 2021). Nucleosome landscapes can thus be viewed as the final result of active positioning forces (the action of chromatin remodelers) combined with destabilizing nucleotide content, including poly(A+T) tracts (see above). Recently, this scenario was proposed to be species-specific (Barnes and Korber, 2021;Oberbeckmann, Krietenstein, et al., 2021), supporting that several combinatorial nucleosome arraying rules can form during the course of evolution. Indeed, when we restricted our analysis to conserved single-copy orthologous fungal genes, the overall profiles and the intensities of "+1 nucleosome" and NDRs were more homogenous between fungi. Similarly, while an NDR can be observed downstream of F. graminearum TTS, it is no longer visible when the analysis is restricted to conserved fungal genes, indicating again an evolutionary component. Promoters are typically found in NDRs upstream of +1 nucleosomes (Yuan et al., 2005;Jiang and Pugh, 2009). NDR sizes largely depend on the action of the SWI/SNF ATP-dependent remodeler RSC (Remodeling the Structure of Chromatin) complex (Krietenstein et al., 2016;Yague-Sanz et al., 2017;Wagner et al., 2020) that seem to facilitate initiation of transcription by preventing the filling of NDRs with nucleosomes (Ocampo et al., 2019).

Nucleosome landscapes of gene units according to gene expression
Same analyses were repeated for genes categorised according to their expression levels (expressed in TPM counts, see Materials and Methods). The general variations in nucleosome profiles around translation start sites are similar for all expression categories in all considered fungi and culture conditions: ATG codons are immediately followed by a well-positioned +1 nucleosome and preceded by a dip in nucleosome density (Figure 8). Remarkable variations are nonetheless observed with regard to positions of +1 nucleosomes and NDRs, as well as the amplitude of the nucleosome signal difference between them (here defined as ∆nucl = |signal+1nucl -signalNDR|), depending on gene expression. ATG-centred nucleosome profiles for genes not expressed in our conditions (TPM = 0) show remarkably reduced ∆nucl when compared to those measured for genes more expressed, and a distance to the ATG reduced (Figure 8 and Supplementary Table  6). Conversely, highly expressed genes (TPM50) display the deepest NDRs located at the furthest upstream the ATG. Similar trends are observed when profiles are centred on B. cinerea and F. graminearum TSS ( Figures 8E and 8F, respectively). Moreover, NDRs were usually found further from the ATG site than those in genes expressed at lower levels or not expressed. Finally, the nucleosome wave signal decay phenomenon was observed at distances from the ATG codons shorter in poorly expressed genes than in highly expressed genes, although the +1 nucleosome remained fairly well conserved. Nucleosome depletion can be the side-effect of active transcription with the binding of pre-initiation complex resulting in nucleosome eviction as previously shown in yeast (Venters and Pugh, 2009). Indeed, we evidenced valleys of nucleosome signals upstream of ATG codons found deeper in highly expressed genes. On the contrary, genes not expressed showed little or no NDR, depending on the considered species or culture condition. In our conditions, the amplitude between the NDR and the +1 nucleosome seemed to be an informative measure of gene expression level: the higher this value is, the more genes are expressed. This feature was less strict when TSS were considered, raising the question of different mechanisms of transcription regulation depending on gene unit structures. A possible scenario is that other factors, such as the general regulator Reb1 (RNA polymerase Enhancing Binding protein), that modulate the action of remodeling factors (Ghassabi Kondalaji and Bowman, 2021) may contribute to the profiles observed here. Strikingly, a recent study evidenced the role of such factors as barriers to fine-tune the action of remodelers, with the consequence of modulating nucleosome spacing and phasing distances (Oberbeckmann, Niebauer, et al., 2021). Oberbeckmann and Colleagues proposed a model consisting in promoter NDRs (maintained by the RSC complex, see above) insulated upstream by the barrier factor Reb1 and downstream by the +1 nucleosome, the spacing between the two landmarks being controlled by the remodelers INO80 or ISW2, maintaining longer vs. shorter distances respectively (Oberbeckmann, Niebauer, et al., 2021). In this model and consistently with our observations, gene bodies are characterized by high nucleosome density with shorter NRLs. Here, the marked oscillations of the signal obtained for highly expressed genes suggests that the presence of well-arrayed nucleosomes around ATG combined with high amplitudes and longer distances of nucleosome signal between the NDR valley and the +1 nucleosome peak could be a hallmark of robust active gene expression. On the contrary, an absent or virtually absent NDR upstream of ATG may mark genes with variable expression levels, often exhibiting a TATA-box in their promoters (Tirosh and Barkai, 2008). Such promoters displayed enhanced sensitivity to mutations in yeast, an observation arguing for a link between chromatin structure and the evolution of gene expression (Hornung, Oren and Barkai, 2012). A further confounding element is the observation that introducing Reb1 binding sites in such promoters reduced sensitivity to mutation (Hornung, Oren and Barkai, 2012), which was interpreted in 2012 (when the work was published) as blocking nucleosome formation and introducing an NDR. In the light of the model recently proposed by see above), this analysis may now be re-visited as providing the necessary binding factor that constraints NDR formation and maintenance by RSC and INO80/ISW2.
At the end of gene units, wave signal decay was also observed downstream of stop codons in all studied fungi (Supplementary Figure 9). In Lmb and Lml, strong nucleosome positioning was found immediately upstream the stop codon, directly followed by a strong nucleosome valley in the case of Lml or a general decrease in signal in Lmb (Supplementary Figure 9A and 9B). In B. cinerea and F. graminearum, nucleosome signal was more defined with clear oscillations decaying past the stop codon ( Supplementary Figures 9C  and 9D). In F. graminearum, a lesser signal intensity seemed to characterise highly to moderately expressed genes (Supplementary Figure 9D), a feature visible only for highly expressed genes in B. cinerea (Supplementary Figure 9C). Altogether, our observations highlight the conservation of a nucleosome immediately before the stop codon, followed by a decrease in signal as a mark of gene expression in the four studied fungi. Markedly, when TTS of B. cinerea were considered, not expressed genes revealed a remarkably regular signal of weak amplitude. On the contrary, nucleosome signal around F. graminearum TTS of genes not expressed were characterized by a strong wave signal of well-arrayed nucleosomes (Supplementary Figure 9F). Recently, RSC was found to also play a positive effect on transcription termination, in addition to initiation (Ocampo et al., 2019). Considering that residence time of RNA polymerase II at the NDR downstream the TTS facilitates fast re-initiation of transcription (by recycling the RNA pol II; (Shandilya and Roberts, 2012;Cole et al., 2014), a scenario could be that RSC may be a transcription tuning knob, its presence stimulating transcription initiation while promoting RNA pol II dissociation to terminate transcription. Accordingly, the signal fuzziness we observed around the stop codons and TTS may reflect the highly dynamic nature of nucleosome positioning/removal catalyzed by RSC to regulate transcription termination. Overall, our results showed strong association of nucleosome landscapes at gene unit boundaries with expression levels in the four studied ascomycetes.

Conclusion
The present study explored nucleosome landscapes of four phytopathogenic filamentous fungi with the aim of unravelling common features as well as potential specificities. Our general observation of nucleosome positioning genome-wide revealed shorter nucleosome-repeat lengths in L. maculans 'brassicae' and L. maculans 'lepidii' compared to B. cinerea and F. graminearum, suggesting a more compact chromatin fibre. High nucleosome occupancy was further observed in AT-rich regions of L. maculans 'brassicae', a feature a priori unexpected considering the well-described destabilising effect of AT-stretches but in line with the heterochromatic nature of these peculiar regions. High nucleosome occupancy was also observed at the loci of BOTY retrotransposons in the genome of B. cinerea. On the contrary, regions with reduced occupancy were observed in F. graminearum and co-localised with highly polymorphic regions described as prone to genetic evolution. As a whole, our results plead in favour of evolution of not only the positions of nucleosomes but also their occupancy, both likely hard-wired to genome sequence evolution, with regions defined by higher occupancies possibly securing sequence conservation. Evolution of genome sequences marked with peculiar chromatin signature profiles in relation with host adaptation has been previously described in the fungal pathogen V. dahliae (Cook et al., 2020). Considering how gene expression may relate on nucleosome patterning, an element of fungal specialization may rely on how chromatin remodeling proteins as well as promoter and other underlying genomic sequences have diversified. 172.0 0.0 1 3.94e-16 † Fitting was performed using the positions of the first four successive peaks of phase frequencies (peaks #1 to #4); ‡ fitting in replicate #2 was performed using the positions of four successive peaks excluding peak #1 for which an apex was not clearly visible at the beginning of the profile (i.e., peaks #2 to #5) 163.6 0.6 1 1.57e-5 † Fitting was performed using the positions of the first four successive peaks of phase frequencies (peaks #1 to #4); ‡ fitting was performed using the positions of three successive peaks excluding peak #1 for which an apex was not clearly visible at the beginning of the profile (i.e., peaks #2 to #4); *as previously defined (Diolez et al., 1995);**as defined by Laurent et al. (2017) (Laurent et al., 2017).