The histone code of the fungal genus Aspergillus uncovered by evolutionary and proteomic analyses

Chemical modifications of DNA and histone proteins impact the organization of chromatin within the nucleus. Changes in these modifications, catalysed by different chromatin-modifying enzymes, influence chromatin organization, which in turn is thought to impact the spatial and temporal regulation of gene expression. While combinations of different histone modifications, the histone code, have been studied in several model species, we know very little about histone modifications in the fungal genus Aspergillus, whose members are generally well studied due to their importance as models in cell and molecular biology as well as their medical and biotechnological relevance. Here, we used phylogenetic analyses in 94 Aspergilli as well as other fungi to uncover the occurrence and evolutionary trajectories of enzymes and protein complexes with roles in chromatin modifications or regulation. We found that these enzymes and complexes are highly conserved in Aspergilli, pointing towards a complex repertoire of chromatin modifications. Nevertheless, we also observed few recent gene duplications or losses, highlighting Aspergillus species to further study the roles of specific chromatin modifications. SET7 (KMT6) and other components of PRC2 (Polycomb Repressive Complex 2), which is responsible for methylation on histone H3 at lysine 27 in many eukaryotes including fungi, are absent in Aspergilli as well as in closely related Penicillium species, suggesting that these lost the capacity for this histone modification. We corroborated our computational predictions by performing untargeted MS analysis of histone post-translational modifications in Aspergillus nidulans. This systematic analysis will pave the way for future research into the complexity of the histone code and its functional implications on genome architecture and gene regulation in fungi.


InTRoDuCTIon
Inside eukaryotic cells, the membrane-bound nucleus facilitates the storage and replication of the genetic material as well as the regulation of gene expression [2]. Within the nucleus, DNA is arranged as a chromatin fibre formed by repeated nucleosome units consisting of histone protein octamers (two copies each of H2A, H2B, H3 and H4 histones) that are wrapped by 145-147 bp of DNA [3,4]. The chromatin structure is regulated by chemical modifications such as methylation of DNA on cytosines (5mC) and adenines (6mA), and different post-translational modifications (PTMs) on the histone proteins, including methylation and acetylation [5,6]. The combination of these different modifications, the histone code, influences the recruitment of diverse nucleosome-associated proteins (e.g. DNA binding proteins), which can cause dynamic transitions between the condensed and transcriptionally silent heterochromatin and the accessible and transcriptionally active euchromatin [7][8][9]. Chromatin dynamics is therefore key for organisms to adapt gene expression patterns in response to a variety of developmental and environmental clues.
In recent decades, several enzymes involved in orchestrating chromatin dynamics have been identified and functionally characterized in diverse eukaryotes [5,9,10]. These enzymes can be broadly classified into three distinct types: (i) writers that place modifications on nucleosomes or linker DNA; (ii) erasers that remove these modifications; and (iii) readers that first recognize specific modifications and subsequently recruit writers, erasers or other proteins to the correct genomic location [11]. Chromatin modifiers that catalyse the same type of chemical modification often harbour functionally and evolutionarily conserved catalytic domains. For instance, DNA methylation (5mC) is catalysed by different DNA methyltransferases (DNMTs) which contain a DNA methyltransferase domain, while histone mono-, di-and tri-methylation are nearly exclusively catalysed by proteins containing a SET [Su (var) 3-9, E(z) and Trithorax] domain [10]. Next to these catalytic domains, the enzymes often contain accessory domains, which can contribute to the overall stability of the enzyme or establish interactions with other proteins [12]. Furthermore, these multi-domain proteins, together with other accessory protein subunits, can form multi-protein complexes to find the correct target region within the genome and modify the exact amino acid residue, collectively contributing to the specific functions of these complexes [13,14]. For example, the histone methyltransferase complex PRC2 (Polycomb Repressive Complex 2) is responsible for the mono-, di-and tri-methylation of histone protein H3 on the lysine at position 27 (abbreviated as H3K27me1/2/3; the same abbreviation principle will be applied below for other histone PTMs) [15][16][17][18]. Evolutionary analyses have suggested that the key enzymes responsible for chromatin modifications have evolved and siginificantly expanded early in eukaryotic evolution [19]. An analysis of 94 gene families with roles in chromatin modifications reveals that 87 out of 94 gene families emerged more than one billion years ago, and 48 were already present in the last eukaryotic common ancestor [19].
Among eukaryotes, the fungal kingdom is remarkably large and diverse, comprising an estimated 1.5-5 million species with various morphologies, life cycles and ecological niches [20]. The most well-studied model fungal species in respect of chromatin modifications is the unicellular yeast Saccharomyces cerevisiae, in which 75 histone PTMs on 43 distinct amino acid positions have been reported, with histone methylation and acetylation being most abundant [11]. In addition to ancient chromatin modifiers, lineage-specific innovations such as duplications and losses were also found in fungi. For example, Dim-2 and RID [Repeat-Induced Point Mutation (RIP)-Defective] are considered to be fungal-specific DNMTs, while DNMT3 that is typically found in animals and plants seems to be absent in fungi [21][22][23]. However, S. cerevisiae is not representative of the complete fungal kingdom, and detailed studies on PTMs are needed for other, less studied fungal species.
The fungal genus Aspergillus comprises approximately 350 species, including pathogens of humans, animals and crops, food contaminants, as well as important cell factories for industrial and medical applications [24]. Chromatin dynamics influences gene expression in several Aspergilli and regulates growth, sexual development, secondary metabolite biosynthesis and virulence [11,[25][26][27]. Thus far, very few well-known chromatin modifications and their catalytic protein complexes have been identified and studied in just a few Aspergillus species [11]. For example, gene knockout of the H3K9 methyltransferase Clr4

Impact Statement
Uncovering how fungi regulate gene expression and genome organization is an important step to understand how they control biological processes such as growth or biosynthesis of bioactive molecules. Despite the known importance of chromatin modifications in controlling a plethora of biological processes across eukaryotes, they remain particularly poorly understood in filamentous fungi, even in model organisms such as Aspergilli. Our systematic evolutionary analysis provides a robust framework for the functional analysis of chromatin modifications in Aspergillus species and other fungi. Our results not only implicate candidate enzymes as playing a role in new chromatin modifications, but they also point at species that have experienced duplications or losses of genes encoding enzymes for specific chromatin modifications. We therefore expect that this work will set the stage for future research into the complexity of the histone code and its functional implications on gene regulation and genome organization in fungi.
homologue in Aspergillus fumigatus reduced histone H3K9me1/3 results in reduction of radial growth and conidia production, as well as delayed conidiation [26]. Similarly, CclA, a subunit of the COMPASS (COMplex of Proteins ASsociated with SET1) complex, is essential for H3K4me2/3 in A. fumigatus, and the deletion strain showed increased secondary metabolite production and decreased growth [28]. The SAGA (Spt-Ada-Gcn5 Acetyltransferase) complex, which is responsible for histone acetylation, has been functionally characterized in Aspergillus nidulans [29]. Notably, this complex in A. nidulans lacks the Sgf11 and Sus1 subunits that have been described in S. cerevisiae, suggesting that this complex evolves in distinct trajectories in different fungi [29]. Besides these few examples, we so far lack a complete view of the occurrence and evolution of chromatin modifiers in the genus Aspergillus. Here, we report the phylogenetic analyses of core enzymes involved in the deposition and removal of histone PTMs, and assess the conservation of all the subunits belonging to chromatin regulator complexes in 94 diverse Aspergilli. We also performed proteomics analyses to corroborate our computational predictions and reveal a broader spectrum of histone modifications.

Phylogenetic analyses of chromatin modifiers
Phylogenetic trees for six major chromatin modifiers were generated using their conserved catalytic domains. Catalytic domain HMMs (Hidden Markov Model) of DNA methyltransferase domain (PF00145), histone methylation SET domain (PF00856), acetyltransferase domain (PF00583), histone acetyltransferase MOZ/SAS FAMILY domain (PF01853), histone deacetylase domain (PF00850) and histone demethylase Jumonji C (JmjC) domain (PF02373) were downloaded from the Pfam database (https://pfam.xfam.org). Hmmsearch, which is part of the HMMER package (v3.1b2) [30], with the Gathering Cut-Off threshold (-cut_ga) was used to identify occurrences of these domains in the predicted proteomes of all 109 fungal species considered here. To filter out short and incomplete fragments, we applied a cutoff calculated by 50 % length coverage of the query sequence (domain) or the hit. Subsequently, the identified domain sequences were aligned using MAFFT v7.271 [31]. TrimAl v1.2 (-gt 0.1) was used to remove positions in the alignments with gaps in more than 90 % of the sequences [32]. Maximum-likelihood phylogenetic trees were reconstructed using IQ-TREE v1.6.10 [33,34] using ultrafast bootstrap and Shimodaira-Hasegawa approximate likelihood ratio test (SH-aLRT) [35] with 1000 pseudo-replicates, and standard model selection to automatically determine the best-fit model. We predicted additional domains in the analysed proteins using Hmmsearch with all Pfam domain profiles that are present in the Pfam database (v28.0) and visualized them along the phylogenetic trees using iTOL [36]. Lastly, additional TBLASTN searches were performed on the NCBI website (https:// blast.ncbi.nlm.nih.gov/Blast.cgi) and checked manually to confirm dubious absences in the domain trees.

Species tree construction
To construct a species phylogeny of the fungi analysed fungi here, we used BUSCO v4.0.1 (Benchmarking Universal Single-Copy Ortholog assessment tool) [37] to identify the occurrence of 758 fungal single-copy genes (Table S2) within the 109 predicted proteomes (Table S1). The identified homologues of single-copy BUSCO genes were subsequently concatenated into a single super-alignment. Then, the maximum-likelihood species phylogeny was estimated using IQ-TREE v1.6.10 [33,34] by employing partitioned model selection, which allows selecting a substitution model for each BUSCO gene. Branch supports were obtained using ultrafast bootstrap and SH-aLRT [35] with 1000 pseudo-replicates separately.

Chromatin modifier subunits homology search and gene tree construction
We obtained previously studied chromatin modifier complexes from the literature, and the associated protein sequences were downloaded from NCBI (https://www.ncbi.nlm.nih.gov) (Table S4). For each subunit of the chromatin modifier complexes, homologues in the analysed 109 predicted proteomes were detected based on sequence similarity searches. First, BLASTP (v2.2.31+) [38] searches were performed, and an e-value threshold of 10 −2 was applied. Multiple blast matches (High-Scoring Pairs: HSPs) of the same query were concatenated if they were aligned with different non-overlapping regions. Only the non-redundant region was retained when HSPs of the same query were overlapping. For these two conditions, the highest e-value of the different matches was used to represent the merged match. Only matches with the e-value lower than 10 −5 and covering more than 50 % of the query or subject were retained for the following analyses. If we observed fewer than ten hits for the BLASTP search after filtering in Aspergilli, we used PSI-blast (Position-Specific Initiated blast v2.2.31+) [38] with five iterations (-num_iterations 5) using the best Aspergilli hit from the BLASTP as the query; if no Aspergilli match was retained, we used the same query as for the BLASTP searches. The same filter criterion as for the BLASTP selection was applied to the PSI-blast hits, and the results of these two homology search strategies were subsequently joined.
To account for apparent protein absences due to erroneous gene annotation, we tested these potential losses by using TBLASTN (blast v2.2.31+) [38] to query protein sequences against the genome assemblies. Matches with e-value lower than 10 −5 together with 4500 nt upstream and downstream were extracted, and Exonerate v2.2.0 [39] was used to predict the gene structure [options: protein2genome, extensive search, showtargetgff, showcigar and showquerygff options were set as true, and '>%ti(%tcb -%tce)\ n%tcs\n' option was used to output the complete aligned coding sequence]. Transeq v6.6.0.0 [40] was used to translate the coding sequence into protein sequences. Together with the protein sequences identified from BLASTP and PSI-blast, all protein sequences from each subunit were retrieved to build a phylogenetic tree (MAFFT, TrimAl and IQ-TREE settings were the same as for the tree reconstruction using the catalytic domains above).
For each sequence included in the phylogenic trees, the domain organization was identified and visualized as described above. When the domain organization of an Aspergillius protein differed from its orthologues, the corresponding gene model was manually checked in the MycoCosm genome browser [1] and corrected. The presence and absence patterns of 83 subunits from 15 different chromatin modifier complexes were visualized using the Seaborn library in Python [41,42]. Lastly, additional TBLASTN searches for dubious absences in the outgroup species were performed on the NCBI website as absences could be caused by older and more fragmented assemblies in the JGI database compared to NCBI.

Histone enrichment and mS
Spores of 4-or 5-day-old sporulating A. nidulans (AnWT pabaA1 veA1, kindly provided by Prof. Joseph Strauss) growing on ME (malt extract) agar medium were harvested with 2 ml sterile water and cultured in 50 ml ME liquid medium at 28°C and 200 rpm for 16 h. The sample was centrifuged at 8000 g for 10 min to discard supernatant, and germinating spores were resuspended in 25 ml 0.8 M NaCl and centrifuged again. Ten millilitres of protoplasting solution (20 mg ml −1 Trichoderma lysing enzyme and 5 mg ml −1 driselease in 0.8 M NaCl) were used to resuspend spores, which were then incubated at 30°C, 100 rpm for 2-3 h. Protoplasts were retrieved by filtering through sterile Miracloth and were concentrated by centrifugation at 3000 rpm and 4 °C for 5 min. Histone proteins were extracted according to an adapted histone protein extraction protocol [43]. Briefly, 10 8 protoplasts were resuspended in 1 ml Nuclei Isolation Buffer (1×PBS, 0.5 mM PMSF, 5 µM leupeptin, 5 µM aprotinin and 5 mM Na-butyrate) with 1 % Triton X-100, followed by vigorous pipetting through a 200 µl pipette several times. This was then centrifuged at 2300 g for 15 min at 4°C, the supernatant was discarded, and the nuclear pellet resuspended in 200 µl Nuclei Isolation Buffer with 0.1 % SDS and 250 U Pierce Universal Nuclease and incubated at 37°C for 30 min to digest nucleic acids. The protein concentration was measured using the Bradford Protein Assay [44].
The samples were separated on a 17 % SDS-PAGE gel and histone bands were excised, chemically acylated with propionic anhydride and in-gel digested with trypsin, followed by peptide N-terminal derivatization with phenyl isocyanate (PIC) [45]. Chemical acylation of lysines (which occurs on unmodified or mono-methylated residues) impairs trypsin cleavage, resulting in proteolityc cleavage at arginine residues only. This treatment allows us to obtain histone peptides of proper length for MS analysis. The samples were then desalted on handmade StageTips columns [46]. Peptide mixtures were separated by reversed-phase chromatography on an EASY-Spray column (Thermo Fisher Scientic), 25 cm long (inner diameter 75 µm, PepMap C18, 2 µm particles), which was connected online to a Q Exactive Plus instrument (Thermo Fisher Scientific) through an EASY-Spray Ion Source (Thermo Fisher Scientific), as previously described [45].
The acquired RAW data were analysed using the integrated MaxQuant software v.1.6.10. The Uniprot UP000000560 database was used for identification of A. nidulans histone peptides. Enzyme specificity was set to Arg-C. The estimated false discovery rate (FDR) was set at a maximum of 1 %. The mass tolerance was set to 6 ppm for precursor and fragment ions. Two missed cleavages were allowed, and the minimum peptide length was set to four amino acids. Min. score for modified peptides and min. delta score for modified peptides were set to 1. Variable modifications included lysine propionylation, monomethylation +propionylation, dimethylation, trimethylation and acetylation. N-terminal PIC labelling was set as a fixed modification [45]. To reduce the search time and the rate of false positives, with increasing the number of variable modifications included in the database search, the raw data were analysed through multiple parallel MaxQuant jobs [47], setting different combinations of variable modifications. Peptides identified by MaxQuant with Andromeda score higher than 50 and localization probability score higher than 0.75 were quantified, either manually or by using a version of the EpiProfile 2.0 software [48] adapted to the analysis of histones from A. nidulans. Identifications and retention times were used to guide the manual quantification of each modified peptide using QualBrowser version 2.0.7 (Thermo Fisher Scientific). Site assignment was evaluated from MS2 spectra using QualBrowser and MaxQuant Viewer. Extracted ion chromatograms (XICs) were constructed for each doubly charged precursor, based on its m/z value, using a mass tolerance of 10 ppm. For each histone modified peptide, the percenatge relative abundance (%RA) was estimated by dividing the area under the curve (AUC) of each modified peptide for the sum of the areas corresponding to all the observed forms of that peptide [49]. The AUC values are reported (Table  S6) and visualized using GraphPad Prism.

RESuLTS AnD DISCuSSIon
DNA methylation, histone methylation and histone acetylation are among the most studied chromatin modifications in animals, plants and fungi [50][51][52]. Protein complexes involved in writing, reading and erasing these modifications have been reported in a few model fungal species to date [16,22], but the evolution and conservation of these complexes in diverse fungal species have not been largely systematically explored. Besides some research on DNA methyltransferases [21,23] and SET domain-containing proteins [17] that have been performed throughout the fungal kingdom, we still lack the evolutionary overview of other chromatin modifiers and in-depth analyses for important fungal groups. To start uncovering the occurrence and evolutionary histories of chromatin modifiers involved in these major modifications in the fungal genus Aspergillus, we first focused on six conserved catalytic domains reported to be involved in histone modifications: DNA methyltransferase domain for DNMTs, SET domain for histone methyltransferases (HMTs), acetyltransferase domain and MOZ/SAS FAMILY domain for histone acetyltransferases (HATs), histone deacetylase domain for histone deacetylases (HDACs) and JmjC domain for histone demethylases (HDMTs). We identified homologues harbouring these conserved domains in 109 predicted proteomes; the completeness of the predicted proteomes was generally high as 97 had more than 97 % complete BUSCO predicted proteins (Table S1). The 94 Aspergilli species cover 18 recognized sections: Polypaecilum, Cremei, Restricti, Aspergillus, Cervini, Clavati, Fumigati, Candidi, Circumdati, Janorum, Terrei, Flavi, Ochraceorosei, Bispori, Cavernicolarum, Usti, Nidulantes and Nigri [53][54][55][56][57][58][59], and 15 other species were used as outgroups (Fig. S1). Two of the outgroup species are Basidiomycota and the remaining species cover each class in the phylum Ascomycota, which provides a robust framework for our evolutionary analyses within Aspergilli. The outgroup species include well-studied model species, some of which are important for crop production, food fermentation and antibiotics production [60,61]. We used the identified homologues for each catalytic domain to determine their evolutionary history using maximumlikelihood phylogenies (Figs 1-4, S2 and S3).

Aspergilli harbour conserved DnA methyltransferases for genome defence and methylation maintenance
The conserved DNA methyltransferase domain is the only domain known to catalyse the addition of methyl groups to DNA specifically at cytosines (5mC) [21]. Based on previous domain architecture analyses, phylogeny reconstructions and functional experiments, five distinct protein families containing a DNA methyltransferase domain have been recognized in fungi: DNMT1, DNMT2 and DNMT5 are commonly found in eukaryotes, while Dim-2 and RID are considered to be fungal-specific [21]. Among these five families, DNMT2 seems to function as a tRNA methyltransferase rather than DNA methyltransferase [62]. In line with previous observation [21,63], DNMT2 only occurs in our study in two Basidiomycota and one Ascomycota (Schizosaccharomyces pombe), suggesting it has been lost in most Ascomycota.
Three (DNMT1, DNMT5 and Dim-2) of the four remaining families catalyse the maintenance of DNA methylation pattern from the parental DNA strand onto synthesized daughter DNA strands, while Dim-2 is also able to catalyse de novo methylation [64,65]. Based on our species selection and phylogenetic analyses, we here found DNMT1 and Dim-2 to be absent in all Aspergilli and only present in four or seven out of 15 outgroup species, respectively, while DNMT5 is found in most Aspergilli and outgroup species (Fig. 1). Only four outgroup species harbour DNMT1, the two Basidiomycota species, Botrytis cinerea and Plectania melastoma, which corroborates previous findings [21,66]. Dim-2 was identified only in seven Ascomycota outgroup species, and matches in B. cinerea and Zymoseptoria tritici were included after an additional TBLASTN search in NCBI (Table S6). The original Dim-2 absence in B. cinerea by local Hmmsearch is caused by the older version and fragmented assembly at JGI compared to NCBI. The Z. tritici IPO323 strain used in this study carries a non-functional Dim-2 as previously reported [67] and no significant hit was found with Hmmsearch or BLASTP searches, but several fragmented hits with TBLASTN searches in NCBI. This loss of Dim-2 function in Z. tritici IPO323 is very recent because other isolates harbour a functional Dim-2 gene [67]. Its loss leads to increased activity of transposable elements (TEs) and decreased mutation rate, thus influencing genome defence and genome evolution [67]. In contrast to the DNMT1 and Dim-2 losses in Aspergilli, DNMT5 is identified in 88 out of 94 Aspergilli and 10 out of 15 outgroup species (Fig. 1), with a conserved domain architecture consisting of a DNA methyltransferase domain, an SNF2 family domain (PF00176), and a Helicase conserved C-terminal domain (PF00271), the latter two domains being involved in transcription regulation [68]. The five outgroup species without DNMT5 are Saccharomyces cerevisiae, Schizosaccharomyces pombe, Neurospora crassa, Fusarium graminearum and Cladonia grayi, among which the two yeast species also lack the other two DNA methylation-maintaining enzymes, and thus accordingly are considered to be incapable of DNA methylation [69][70][71]. Five out of the six Aspergillus species that lost DNMT5 belong to a single clade within the Nidulantes section (A. undulatus, A. olivicola, A. filifera, A. venezuelensis and A. stella-maris), and the sixth one is A. robustus which is assigned to the Circumdati section (Fig. S1). In the Basidiomycota fungus Cryptococcus neoformans, DNMT5 maintains DNA methylation using the same cofactors as DNMT1 (70,72). Thus, it is most likely that DNMT5 is similarly involved in DNA methylation maintenance in Aspergilli as in Cryptococcus neoformans. Moreover, DNA methylation maintenance does not appear to be required for the survival of fungi, as shown with the absence of all three of DNMT1, DNMT5 and Dim-2 not only within the Nidulantes section but also in Saccharomyces cerevisiae and Schizosaccharomyces pombe (Fig. 1). However, the absence of these genes could have an impact on the regulation of certain genes such as secondary metabolite biosynthetic gene clusters and it would be interesting to compare species within the Nidulantes section in more detail to identify how such loss may have affected gene expression and genome organization. RID is a more recently acquired single-domain DNMT that specifically occurs in fungi and is involved in a fungal-specific genome defence mechanism called RIP (repeat-induced point mutation) [21]. RID is a putative de novo methyltransferase for TEs in both CG and non-CG contexts relating to RIP during the sexual stage in N. crassa; however, the capacity of RID to catalyse 5mC has not yet been demonstrated [72]. Here, we identified RID in all analysed Aspergilli species as well as in 10 outgroup species. A RID-like homologue knock-out mutant in the model filamentous Ascomycota Podospora anserina shows blocked sexual development, as gamete formation and fertilization can happen but dikaryotic cells cannot develop normally [73]. Interestingly, the knock-out mutant of DmtA (homologue of RID in A. flavus) shows altered conidiation, sclerotial production, aflatoxin biosynthesis and virulence [27], suggesting that RID might have additional functions next to genome defence. Additionally, fungi lost the gene for  de novo DNA methyltransferase, DNMT3, around 50-150 million years ago [74,75], and only DNMT5 and RID were consistently found in Aspergilli (Fig. 1). Thus, we hypothesize that Aspergilli can only maintain DNA methylation via DNMT5 but do not spur genome-wide de novo DNA methylation. This hypothesis seems to be corroborated as bisulphite sequencing in A. flavus has uncovered very little genome-wide DNA methylation [76]. However, systematic analyses of DNA methylation patterns and experimental evidence for the roles of DNMTs in Aspergilli are still largely missing.

Aspergilli lack SET7 (KMT6), but SET-domain proteins are abundant and other histone methyltransferases are conserved
Histone methylation is typically catalysed by SET domain-containing proteins [17,77]. We identified a total of 1548 SET domaincontaining proteins and classified them into 15 groups based on their domain architecture and phylogenetic support that contain species both from Aspergilli as well as from the outgroup. Among these groups, 107, 111, 109, 102 and 108 sequences represent homologues of the known HMTs SET9 (KMT5), Dim-5 (KMT1), SET1 (KMT2), Ash1 (KMT2H) and SET2 (KMT3), respectively   S1) without Fusarium oxysporum and Aspergillus kawachii since they were excluded in the latter analysis. Each row represents a subunit in the histone modification complexes; the order agrees with the summary list (Table S5). Each row shows the presence/absence of a specific complex subunit. For each complex, the first row represents the catalytic subunit. Each box in the matrix indicates the presence (dark colour) or absence (light colour) of a subunit (Table S4). The brown boxes indicate additional hits found with TBLASTN searches on NCBI (Table S6). (Fig. 2). The remaining homologues were divided into nine groups that we broadly classified as SET domain-containing proteins. Among them, eight groups contain a SET domain but lack experimentally determined functions, and 117 sequences could be grouped as ribosomal lysine N-methyltransferases (Fig. 2). Enzymes in this group contain a SET domain together with a Rubisco LSMT substrate-binding domain (PF09273). While the homologue in plants transfers a methyl group to the large subunit of the Rubisco holoenzyme complex [78,79], the homologue in Saccharomyces cerevisiae has been characterized as a ribosomal lysine N-methyltransferase [80].
The well-described HMTs, SET1, Dim-5, SET2, Ash1 and SET9, all contain distinct accessory domains and catalyse methylation on lysine residues at different positions of H3 or H4 tails [17]. SET1 is composed of the SET domain, the N-SET domain (PF11764) as well as the SET_assoc domain (PF11767), and is responsible for H3K4me1/2/3 in Saccharomyces cerevisiae [81]. SET1 occurs in a single copy in all species (Fig. 2), indicating that this enzyme and its catalysed PTM play fundamental roles in fungal biology. Dim-5 contains a SET domain and Pre-SET domain (PF05033) and acts as an H3K9 methyltransferase inducing heterochromatin formation notably in N. crassa [82]. Dim-5 was found in all species included in this study except Saccharomyces cerevisiae and A. calidoustus, among which Saccharomyces cerevisiae is reported to lack H3K9me [82,83]. How such a loss has  Aspergilli. (a). Schematic presentation of the in-gel digestion method used to process histones isolated in A. nidulans. After SDS-PAGE separation, gel bands corresponding to histone H3 and H4 are excised and chemically acylated with propionic anhydride (which derivatizes unmodified and mono-methylated lysines). Thanks to the derivatization, trypsin cutting at lysines is impaired, and histones are cleaved only at the C-terminus of arginine. Finally, digested peptides are N-terminally derivatized with phenyl isocyanate (PIC). (b). A. nidulans histone H3 and H4 sequences. Cleavage sites after PRO-PIC digestion are shown. Peptides identified in this study are highlighted in light red and labeled using their starting-ending amino acids. (c). Relative abundance of post-translational modifications detected on peptides of histone H3 and H4 from A. nidulans using unbiased mass spectrometry analysis. unmod: unmodified peptide; me1: monomethylation; me2: dimethylation; me3: trimethylation; ac: acetylation. The term "combinations" refers to multiple modifications occurring on the same peptide (Table S6). (d). Histone modifications are predicted by genomic and phylogenetic analyses, and/or proved by mass spectrometry in A. nidulans. Histone modifications verified by mass spectrometry are labeled as solid dots, while the ones that are only predicted by phylogenic analysis are labeled as solid dots with a red outline.
impacted gene expression in A. calidoustus in comparison to its closest relative A. carlsbadensis deserves further investigation. Both SET2 and Ash1 catalyse H3K36me in N. crassa and share the SET domain and the AWS domain (PF17907) (Fig. 2), yet SET2 harbours one additional SRI domain (PF08236) that interacts with RNA polymerase II and links H3K36me with transcript elongation [84][85][86]. In N. crassa, Ash1 was shown to contribute to H3K36me at sub-telomeric regions, while SET2 acts primarily at euchromatic regions [87,88]. Considering the common origin of SET2 and Ash1 as well as their similar activities, an ancient duplication followed by the loss or gain of the SRI domain is likely to have contributed to their functional diversification (Fig. 2). SET2 was found in all outgroup species, while Ash1 was lost in both Basidiomycota and two yeasts. False absences of Ash1 in B. cinerea and F. graminearum were corrected as they could be retrieved by TBLASTN search at NCBI (Table S6). Both enzymes are conserved in all Aspergillus species except in A. carbonarius from the Nigri section, suggesting that either other SET domaincontaining proteins exhibit functional redundancy or that H3K36 methylation is not fully required in this species. In Fusarium fujikuroi, deletion of SET2 and Ash1 strongly affected vegetative growth and conidiation [88]. SET9 homologues harbour only a SET domain, and were shown to catalyse H4K20me in Schizosaccharomyces pombe [89]. This group is well conserved in most of the species except in Saccharomyces cerevisiae [90] and in the two closely related Aspergilli, A. novoparasiticus and A. arachidicola that both belong to section Flavi (Figs 2 and S1).
SET7 is the catalytic subunit of PRC2, and is highly conserved in mammals, plants and N. crassa, in which it has been well characterized to catalyse methylation on H3K27, a PTM that is deposited at facultative heterochromatin [16]. Even though PRC2 is considered to be a key chromatin modifier complex, its catalytic subunit SET7 was lost in many fungal species [17]. SET7 has previously been reported to be absent in A. nidulans and A. fumigatus [15][16][17], and we here demonstrate that SET7 is absent from all analysed Aspergilli as well as from Penicillium (Fig. 2), suggesting that this enzyme was probably lost at or before the divergence of the last common ancestor of these two fungal genera. Thus, in line with previous reports [17,25,91,92], our results indicate that methylation on H3K4, H3K9, H3K36 and H4K20 are probably catalysed by conserved enzymes in most Aspergilli, while H3K27 methylation is absent. The loss of this methylation is usually reported to trigger serious phenotypes in humans [93] and mice [94]. Loss of H3K27me3 in fungi typically leads to gene expression changes without severe growth defects [16,95], and thus further investigating how Aspergilli and other fungi are able to cope with such a loss will probably provide new insights on the histone code and its evolution in fungi.
Importantly, the occurrence of a high number of uncharacterized SET domain-containing proteins opens the possibility that these contribute to not yet characterized histone modifications. In particular, four uncharacterized SET domain-containing groups with 144,195,113, and 113 members are conserved in Aspergilli and thus probably play important roles in different cellular processes.

Histone acetyltransferases are conserved in Aspergilli
HATs have been classified into five families: MYST (MOZ, Ybf2/Sas3, Sas2, and Tip60)-related HATs, Gcn5-related acetyltransferases (GNATs), p300/CBP HATs, general transcription factor HATs and nuclear hormone-related HATs [96][97][98]. We focused here only on the first two HAT families because these are known to be involved in the regulation of histone proteins and appear largely conserved in fungi [96]. The MOZ/SAS FAMILY domain is the catalytic domain of MYST-related HATs [99]. As previously reported [100], the three characterized MYST-related HATs, Sas3, Sas2 and Esa1, belong to distinct monophyletic clades, and these paralogues are nearly fully conserved in the fungi analysed here (Fig. 3). MOZ/SAS FAMILY domain and zf-MYST domain (PF17772) are present in Sas3, which can on its own acetylate histone H3 and H4, and H2A weakly in vitro [101]. Within the NuA3 (Nucleosomal Acetyltransferase of histone H3) complex, Sas3 histone acetylation activity is in vitro restricted on H3K14 and to a less extent on H3K23, while in Saccharomyces cerevisiae, only H3K14ac can be detected [102]. Sas3 occurs in most of the species but is lost in Schizosaccharomyces pombe, as previously observed [103,104], and is duplicated in A. carbonarius (Fig. 3). In addition to the two domains found in Sas3, Sas2 and Esa1 carry a PHD_4 domain (PF16866) or Tudor-knot domain (PF11717), respectively (Fig. 3). Esa1 is present as a single copy in all species included in this study, and its functions are known to contribute to the acetylation of K5, K8, K12 and K16 on histone H4, as well as on H2AZK14 in Saccharomyces cerevisiae [105][106][107][108]. Sas2 is similarly found in most of the species except for two Basidiomycota as well as for A. carbonarius. Sas2 is involved in chromatin-mediated gene regulation and H4K16 acetylation in both Saccharomyces cerevisiae and Candida albicans [105,109,110]. Although Esa1 and Sas2 catalyse H4K16ac, it was shown in Candida albicans that they are differentially recruited at different stages of development [105]. Because both Esa1 and Sas2 are well conserved in the genus Aspergillus, we hypothesize that these two HATs may also have distinct roles during development.
GNATs have an acetyltransferase domain as their catalytic core, and 2604 acetyltransferase domain-containing proteins were found in our search (Fig. S2). The genes for two well-characterized GNATs, Gcn5 and Elp3, occur in a single copy in every species included in our study (Fig. S2). Gcn5 harbours a Bromodomain (PF00439) in addition to the Acetyltransferase domain and is a crucial subunit of distinct HAT complexes such as the HAT-A2 complex, the SAGA complex and the ADA complex in Saccharomyces cerevisiae [96,98]. All these complexes contain two coactivator proteins, Ada2 and Ada3, that directly bind to Gcn5 [111] and the recruitment of different other subunits determines the location where they catalyse acetylation [96]. Elp3 is characterized by the Acetyltransferase domain, Radical SAM (PF04055) domain and Radical SAM C-terminal domain (PF16199), and has been shown to catalyse H3K14ac and H4K8ac in Saccharomyces cerevisiae [112][113][114]. Elp3 also exhibits a dual enzymatic activity as it catalyses demethylation in mouse zygotes [115], which is due to the Radical SAM domain that catalyses the demethylation of methylated lysyl residues on histones [116].
Next to Gcn5 and Elp3, we have identified 2386 additional acetyltransferase domain-containing proteins. Because their functions remain unknown due to the lack of homology to sequences with characterized functions, they are valuable candidates that may acetylate histones or other proteins. Interestingly, it has been shown that Gcn5 also exhibits succinyltransferase activity in humans [117] and can also function as a histone crotonyltransferase to regulate gene transcription in vitro [118]. It is expected that other less common histone modifications, i.e. histone lysine butyrylation, propionylation or malonylation, may also be catalysed by histone acetyltransferases because these histone modifications all use acyl-Coenzyme A to complete the reactions [119][120][121]. Further research is needed to determine whether these less common histone modifications are catalysed by the known histone acetyltransferases Gcn5 and Elp3, or as of yet unstudied HATs. In conclusion, all five characterized fungal HATs described above are highly conserved in the genus Aspergillus, showing their importance in fungal biology. Many other acetyl-transferase domain-containing proteins are encoded in Aspergillus genomes and may be responsible for histone acetylation or additional less common histone modifications such as butyrylation, succinylation and malonylation.

Duplication of RpdA histone deacetylase gene in the Flavi and Janorum sections
Five families of HDACs can be grouped into two classes based on sequence similarity: class I (HOS1, HOS2 and Rpd3/RpdA) and class II (HdaA and HOS3) [122,123]. They all harbour the Histone deacetylase catalytic domain, and our phylogenetic tree corroborates the HDAC classification into two distinct classes (Fig. 4). HdaA is found as a single copy in each species included in this study (Fig. 4) . It harbours the Histone deacetylase domain and C-terminal Arb2 domain (PF09757) that serves as an anchor to target the centromeric heterochromatin region [124,125]. HdaA deacetylates lysine in histones H3 and H2B, but not H4 or H2A in Saccharomyces cerevisiae [126,127], and the ΔhdaA knock-out strains show increased or reduced production of different secondary metabolites in A. nidulans [128], A. niger [128] and A. fumigatus [129]. By contrast, HOS3 is found in only 80 Aspergillus species, as it is lost in nine species (A. nomiae, A. pseudocaelatus, A. tamarii, A. pseudotamarii, A. caelatus, A. sergii, A. novoparasiticus, A. minisclerotigenes, A. arachidicola) in the Flavi section, three species (A. desertorum, A. stercorarius and A. multicolor) in the Nidulantes section, and two species (A. heteromorphus and A. carbonarius) in the Nigri section. HOS3 is well conserved in all outgroup species except Schizosaccharomyces pombe (Fig. 4) [130]. HOS3 has been shown to catalyse deacetylation with distinct specificity for histone H2AK7, H2BK11, H3K14 and H3K23, as well as H4K5 and H4K8 in Saccharomyces cerevisiae cell extracts [131]. Notably, HOS1 was only detected in three outgroup species, Saccharomyces cerevisiae and both Basidiomycota (Fig. 4). The specificity of Hos1 in catalysing deacetylation on histone proteins is not yet known [132], but it can catalyse the deacetylation of the Smc3 subunit of cohesin and influence the regulation of chromosome segregation during mitosis in S. cerevisiae [133]. HOS2 is found in 15 outgroup species and most Aspergilli except A. oryzae (Fig. 4). It is known to deacetylate lysine residues in H3 and H4 histone tails [134]. The conservation of HOS2 in Aspergilli suggests that it plays a key role, which is supported by the deletion of HOS2 in A. niger where the deletion strains displays severe reduction in growth, sporulation, SM biosynthesis and stress resistance [128].
Within class I HDACs, the Rpd3/RpdA clade is divided into two sub-clades (Fig. 4). One is nearly fully conserved (except in A. oryzae) and comprises the well-characterized Rpd3/RpdA HDACs, while the other one only comprises sequences from the Flavi and Janorum sections (Fig. 4). Our phylogenetic reconstruction suggests a duplication event that has occurred in the ancestor of all Aspergilli giving rise to a paralogue that was retained in two sections only and could explain the A. oryzae loss in the conserved sub-clade. The two paralogues from the same species vary in length as they share a conserved N-terminal catalytic domain but have a different C terminus (Fig. 4) . Follow-up studies are needed to investigate the function of both paralogues in Aspergillus species. Rpd3 is a crucial deacetylase in eukaryotes, as it is required for deacetylation in most locations on core histone proteins, except H4K16 in vitro [135] that is catalysed by HOS2 [134]. Saccharomyces cerevisiae lacking Rpd3 is sensitive to high osmolarity and shows compromised expression of osmotic stress genes [136]. Deletion of RpdA in A. nidulans is lethal [137], while only heterokaryon disruptants could be obtained in A. oryzae [138], indicating its crucial role in the fungal genus Aspergillus. In summary, class I and class II deacetylases are well conserved in Aspergilli, and duplication events like the one observed for RpdA may confer new activities. While only one paralogue was retained in A. oryzae, the presence of two paralogues in sections Flavi and Janorum may provide functional redundancy or diversification, which remains to be determined.
Lastly, we reconstructed the phylogeny of histone demethylases harbouring the Jumonji C (JmjC) domain (PF02373) (Fig. S3). Six groups are distinguished by their domain distribution: KDM4A/JHDM3/JMJD2, KDM5/JARID, DMM-1, JMJD6/PKDM11, KDM2/ Epe1 and KDM6B/JMJD3, which also agrees with a previous study [91]. The first five groups occur in most of the outgroup and Aspergillus species, but KDM6B/JMJD3 is only present in several outgroup species but not in any of the Aspergillus species. Certain JmjC domain-containing proteins were shown to catalyse the demethylation of histone proteins, e.g. KDM5/JARID for H3K4me2/3 and KDM4A/JHDM3/JMJD2 for H3K9me2/3 and H3K36me2/3 [139]. Histone demethylase activity has not been detected for DMM-1 [140] or Epe1 [141] in N. crassa, while KDM2 can demethylate H3K36me1/2 in mouse [142] and JMJD6/PKDM11 is an arginine demethylase in Arabidopsis [143]. Interestingly, the human KDM6B/JMJD3 is responsible for the demethylation of H3K27me2/3 [144]. As we did not find SET7 and thus conclude that H3K27me3 is absent, it would be tempting to speculate that H3K27 demethylase should be absent in Aspergillus too. However, based on the fact that other enzymes (e.g. UTX) can also catalyse the demethylation of H3K27 [144,145], we cannot unambiguously correlate the absence of KDMB6B/JMJD3 to SET7 loss. More studies are needed to investigate the regulation mechanism and interplay of proteins for H3K27 methylation.

Presence-absence patterns of subunits that compose histone modifier complexes
Histone PTMs are catalysed by the conserved domains found in HMTs, HATs and HDACs (Figs 2-4 and S2), but these chromatinmodifying enzymes often do not act alone. Instead, they are part of protein complexes that comprise accessory subunits, which are involved in the recognition and induction of specific histone modifications at the correct locations [13,14]. To further investigate the presence-absence patterns (PAPs) of known chromatin modifier complexes in Aspergilli, we focused on 15 chromatin modifier complexes with 83 subunits previously identified and studied in Saccharomyces cerevisiae or N. crassa (Table S4). Forty-four subunits are present in all Aspergillus and outgroup species. Most of the accessory subunits are either fully conserved or absent in Aspergilli (Fig. 5 and Table S4). Ten subunits show PAPs in both Aspergilli and outgroup species, while nine subunits appear specific to Saccharomyces cerevisiae as they are not found in any other species and are likely to play a specific role in this species. Seventeen subunits show PAPs in the outgroup species, while eight are absent and nine are present in the Aspergilli. Lastly, three subunits are present in all outgroup species but show PAPs in Aspergilli.
Five histone methyltransferase complexes, COMPASS, DCDC (Dim-5/-7/-9/CUL4/DDB1 Complex), SET2 Complex, PRC2 and Dot1 (disruptor of telomeric silencing 1) (KMT4), were included in our analyses (Fig. 5). Five COMPASS subunits (SET1, CclA, Swd1, Swd2 and Swd3) are fully conserved, suggesting an important function for H3K4 methylation. Consistent with a PAP in most outgroup species and the loss in Aspergilli, Spp1p is not essential for H3K4 methylation [146]. As Sdc1 shows a PAP in Aspergilli and outgroup species, it is likely that this subunit also does not play a key role in H3K4 methylation. In the DCDC complex, while Dim-7 and Dim-9 seem to be dispensable for H3K9 methylation, Dim-5, Cul4 and Dim-8 are highly conserved, and we thus expect these enzymes to play key roles in the activity of this complex ( Fig. 5 and Table S4). However, it was shown that Dim-7 and Dim-9 are responsible for heterochromatin recognition, which leads to subsequent recruitment of the complete DCDC to induce H3K9 methylation in N. crassa [147]. Thus, it is very likely that other accessory subunits are involved in recruiting DCDC at specific chromatin locations in Aspergilli. Furthermore, the unique absence of Dim-8 in A. indicus may suggest altered H3K9 methylation in this species. Consistent with the absence of SET7 (Fig. 2), all components of the PRC2 complex are lost in Aspergilli (Fig. 5). Lastly, the interaction and crosstalk of multiple chromatin modifications in Aspergillus species should be further explored. For example, heterochromatin formation and regulation are done by H3K9me3, DNA methylation and H3K27me3 in N. crassa. In this species, H3K9me3 is recognized by an adapter protein, HP1, to recruit DNA methyltransferase Dim-2 for DNA methylation [72], and cause the redistribution of H3K27me3 [15]. In Aspergillus species, although HP1 and the writer for H3K9me3, DCDC complex, are well conserved (Fig. 5), the inability for genome-wide de novo DNA methylation and H3K27me3 suggest different mechanisms of heterochromatin formation and regulation in Aspergilli, which remain to be determined.
We included eight HAT complexes that differ by their distinct catalytic subunits (Fig. 5). The HAT-A2, ADA, SAGA and Elongator complexes use an acetyltransferase domain-containing protein as catalytic subunit [29,148]. The SAS (Something About Silencing), NuA4 (Nucleosomal Acetyltransferase of histone H4), and NuA3 complexes use a MOZ/SAS FAMILY domain-containing protein [149]. The NuB4 (Nuclear Hat1p-containing type B histone acetyltransferase) complex uses the histone acetyltransferase HAT1 domain-containing protein [150]. Gcn5, Ada2 and Ada3 are fully conserved in all fungi and constitute the core enzymatic module of the HAT-A2, ADA and SAGA complexes. This catalytic module contributes to catalysing histone acetylation, and other subunits in these complexes guide them to specific locations and link them to the transcriptional machinery [151]. Besides the catalytic module, the ADA complex also harbours the Ahc1 subunit that is found in Saccharomyces cerevisiae only [148,152] . Twelve additional subunits of the SAGA complex are well conserved in Aspergilli and outgroup species (Fig. 5). The SAGA complex also contains a deubiquitylase (DUB) module, for which the composition differs between Aspergilli and other species. Sgf11 is unique to Saccharomyces cerevisiae and Sus1 was probably lost in the ancestor of Pezizomycotina because we can only identify Sus1 orthologues in both Basidiomycota and the yeasts (Fig. 5). The Elongator complex appears more variable as four subunits are mostly conserved, but Elp6 is specific to Saccharomyces cerevisiae, and Hpa2 and Hpa3 show PAPs in both Aspergilli and outgroup species.
The NuA4, SAS and NuB4 complexes are well conserved in Aspergilli as well as in the outgroup species (Fig. 5). These complexes appear slightly different in Saccharomyces cerevisiae with the presence of specific subunits (Eaf5 in the NuA4 complex; Sas4 and Sas5 in the SAS complex). Both SAS and NuA4 complexes can catalyse H4K16ac with opposite effects during cell growth and morphogenesis in Candida albicans [105]. It would be interesting to assess whether these two complexes exhibit similar activities in Aspergilli. For the NuA3 complex, two accessory subunits are well conserved among Aspergilli while two others, Pdp3 and Nto1, show different PAPs, and Yng1 is specific to Saccharomyces cerevisiae. Pdp3 and Nto1 seem to be lost in a complementary manner, suggesting a possible functional redundancy (Fig. 5). Finally, two histone deacetylase complexes, Sin3-Rpd3 and HCHC, are fully conserved (Fig. 5). The Sin3-Rpd3 complex is slightly different in Saccharomyces cerevisiae with the specific Ume1 subunit.
Overall, our analysis showed that the core catalytic modules of complexes involved in histone PTMs are mostly conserved in all fungi, including Aspergilli, but variability in accessory subunits could contribute to differences in the localization of histone PTMs along the genome.

unbiased mS captures the diverse histone modifications in A. nidulans and verifies the absence of H3K27 methylation
Our genomic and phylogenetic analyses have identified the presence of multiple chromatin modifiers in Aspergilli, suggesting that these species have the ability to establish a wide variety of histone modifications. To assess the presence of histone PTMs predicted in our evolutionary analyses (Figs 2-5), we extracted histone proteins from A. nidulans and performed unbiased MS analyses (Fig. 6a, b). H3 and H4 peptides with PTMs corresponding to H3K4me, H3K9me, H3K36me and H3K79me were detected (Fig. 6c), which we anticipated being catalysed by SET1, Dim-5, Ash1 and SET2 and Dot1, respectively. H4K20me, which is catalysed by SET9 in Schizosaccharomyces pombe [89], was not detected (Fig. 6c), probably due to technical reasons because the unmodified form of histone H4 peptide 20-24 containing the K20 residue was not detected. For the acetylation on H3K9, H3K14, H3K18, H3K23, H4K5, H4K8, H4K12 and H4K16, they were confidently detected (Fig. 6c). We did not observe any peptide that would indicate the presence of H3K27me in A nidulans, while we did detect H3K27ac (Fig. 6c and Table S6). These results corroborate our computational analysis that the PRC2 complex is absent in A. nidulans, as well as previous reports on the absence of H3K27me3 in several Aspergilli [17,153,154]. Thus, the MS results verify the occurrence of histone modifications predicted to be present based on our comparative analyses in the fungal genus Aspergillus. Moreover, the differences between the relative abundance of distinct histone modifications provide interesting information to investigate these abundances and identify specific conditions under which these modifications increase or decrease in abundance.

ConCLuSIon
Changes in chromatin organization play an important role in regulating the expression of genes [130]. The fungal genus Aspergillus is well known for its ability to thrive in a broad range of different environmental niches and the huge potential for secondary metabolite biosynthesis [11,155]. However, a detailed and complete overview of occurrence and evolution of chromatin modifiers in Aspergilli is a prerequisite to be able to better understand how these fungi regulate gene expression. In the future, more research is needed to uncover how the absence of de novo DNA methylation impacts genome evolution in Aspergilli. We demonstrate that PRC2 is lost throughout Aspergilli, a finding consistent with the lack of H3K27me detection by MS. Although histone acetyltransferases and deacetylases are overall conserved, our study also revealed specific duplication and loss events, which require further investigation to understand their impact on chromatin and gene regulation in certain Aspergillus species. Moreover, other less common histone modifications such as histone lysine succinylation, crotonylation, butyrylation, propionylation and malonylation, are all probably catalysed by HATs and HDACs as they all use Coenzyme A (CoA) to complete the reactions [119][120][121], but may require specific, yet to be discovered subunits. In conclusion, our study provides a detailed overview of the evolutionary routes of the chromatin modifier complex in the fungal genus Aspergillus, and therefore generates the necessary framework to identify targets for functional studies to understand how chromatin regulates gene expression in Aspergilli.

Funding information
Xin Zhang is funded by the Chinese Scholarship Council (CSC) (201907720028). This work was supported by EPIC-XS, project number 823 839, funded by the Horizon 2020 programme of the European Union.

Acknowledgement
The sequence data were produced by the US Department of Energy Joint Genome Institute (http://www.jgi.doe.gov/) in collaboration with the user community. In particular, we thank Scott Baker, Mikael Andersen and Joseph Spatafora for their permission to use the sequence data. We thank Jos Houbraken for validating the Aspergillus species tree and Joseph Strauss for providing the A. nidulans strain. We thank Alessandro Vai for support with analysis of MS data for histone PTM identification on Aspergillus core histone sequences.

Conflicts of interest
The authors declare that there are no conflicts of interest.