Skip to main content

Loss of EZH2-like or SU(VAR)3–9-like proteins causes simultaneous perturbations in H3K27 and H3K9 tri-methylation and associated developmental defects in the fungus Podospora anserina

Abstract

Background

Selective gene silencing is key to development. It is generally accepted that H3K27me3-enriched heterochromatin maintains transcriptional repression established during early development and regulates cell fate. Conversely, H3K9me3-enriched heterochromatin prevents differentiation but constitutes protection against transposable elements. We exploited the fungus Podospora anserina, a valuable alternative to higher eukaryote models, to question the biological relevance and functional interplay of these two distinct heterochromatin conformations.

Results

We established genome-wide patterns of H3K27me3 and H3K9me3 modifications, and found these marks mutually exclusive within gene-rich regions but not within repeats. We generated the corresponding histone methyltransferase null mutants and showed an interdependence of H3K9me3 and H3K27me3 marks. Indeed, removal of the PaKmt6 EZH2-like enzyme resulted not only in loss of H3K27me3 but also in significant H3K9me3 reduction. Similarly, removal of PaKmt1 SU(VAR)3–9-like enzyme caused loss of H3K9me3 and substantial decrease of H3K27me3. Removal of the H3K9me binding protein PaHP1 provided further support to the notion that each type of heterochromatin requires the presence of the other. We also established that P. anserina developmental programs require H3K27me3-mediated silencing, since loss of the PaKmt6 EZH2-like enzyme caused severe defects in most aspects of the life cycle including growth, differentiation processes and sexual reproduction, whereas loss of the PaKmt1 SU(VAR)3–9-like enzyme resulted only in marginal defects, similar to loss of PaHP1.

Conclusions

Our findings support a conserved function of the PRC2 complex in fungal development. However, we uncovered an intriguing evolutionary fluidity in the repressive histone deposition machinery, which challenges canonical definitions of constitutive and facultative heterochromatin.

Background

Histones are subjected to a variety of post-translational covalent modifications [1] that may impact the overall degree of packing of the genome. In most organisms, the opened euchromatin is enriched in tri-methylation of lysine 4 and lysine 36 (H3K4me3 and H3K36me3), two concomitant modifications associated with active transcription [2]. In contrast, histones present in the compacted heterochromatin are enriched in either tri-methylated H3K27 (H3K27me3) or tri-methylated H3K9 (H3K9me3) [3]. The compact architecture of heterochromatin limits the accessibility of the transcription machinery to the embedded DNA, thereby silencing gene expression. Repeat-rich genomic regions enriched in H3K9me3 are referred to as ‘constitutive’ heterochromatin because the subsequent silencing may be constant across development [4]. In contrast, the ‘facultative’ heterochromatin corresponds to the deposition of H3K27me3 on gene-rich regions, whose silencing is transient and dynamic across developmental processes, allowing cell type-specific differentiation and rapid adaptation of gene expression [5].

H3K9me3 is catalyzed by the SET-domain SU(VAR)3–9 enzymes [6,7,8,9] and bound by the chromodomain of Heterochromatin Protein 1 (HP1) [10,11,12]. Subsequent HP1 oligomerization, which results in nucleosome binding [13] and phase separation [14,15,16], further enhances chromatin compaction and spreads this structure over large genomic compartments. When present, high levels of DNA cytosine methylation are found in H3K9me3-enriched regions. This modification is associated with genome stability by preventing either expression of transposable elements or mitotic recombination. As demonstrated in mouse and drosophila models, when H3K9me3 enzymatic activity is absent or reduced, the embryos die because of diverse developmental defects [17]. Accordingly, H3K9me3-dependent heterochromatin is considered as a barrier to cell fate changes, by preventing some transcription factors to bind DNA [18].

First identified in Drosophila, the conserved Polycomb group (Pc-G) protein complexes were shown to be both writers (Polycomb Repressive Complex 2, PRC2) and readers (Polycomb Repressive Complex 1, PRC1) of H3K27me3 [19,20,21]. The catalytic subunit of PRC2 is made of the SET-domain protein enhancer of zest homolog 2 (EZH2), which tri-methylates H3K27. The H3K27me3-polycomb modification is instrumental for maintaining transcription repression established during early development and, therefore, has been linked to biological functions related to cell fate, as diverse as X-inactivation [22] and hematopoiesis [23]. Conversely, when PCR2 functions are impaired, development of plants [24], insects [25] and mammals [26, 27] is disrupted, presumably because lineage-specific programs are no longer accurate. Mutations affecting EZH2 lead to cell proliferation in cancer [13], while others are associated with congenital disorders such as the Weaver syndrome [28] or Ataxia-telangiectasia disease [29].

In animals and flowering plants, H3K9me3 and H3K27me3 marks are almost always mutually exclusive [30], but the molecular basis for this exclusion is unknown. To further investigate the genomic distribution and functional interplay of these two marks, filamentous fungi of the Pezizomycotina clade are a valuable alternative to higher eukaryotic models. Unlike ascomycete yeasts [8], most of them display activity of both H3K9 and H3K27 histone methyltransferase enzymes [31]. The role of the H3K27 histone methyltransferase was investigated in the model fungus Neurospora crassa [32,33,34,35,36,37,38,39,40], the human pathogen Cryptococcus neoformans [41], the plant endophyte Epichloë festucae [42] and the plant pathogens Fusarium graminearum [43], Fusarium fujikuroi [44], Zymoseptoria tritici [45] and Magnaporthe oryzae [46]. As for metazoans and plants, H3K27me3 is essential for proper development of F. graminearum [43], F. fujikuroi [44] and M. oryzae [46] (including growth, sporulation, fertility, pathogenicity, and expression of mycotoxins, pigments, and other secondary metabolites). Paradoxically, deletion of the PRC2 enzymatic subunit of N. crassa, C. neoformans or Z. tritici does not result in obvious development defects. However, H3K9me3 methylation is required for N. crassa normal growth and full fertility [47, 48] and in the wheat pathogen Z. tritici, lack of H3K9me3 induces genomic rearrangements. H3K9me3 and H3K27me3 modifications do not overlap in these fungal genomes, a pattern also found in animals and plants. Interestingly, depletion of H3K9me3 in either N. crassa or Z. tritici leads to a massive redistribution of H3K27me3 in genomic compartments that were previously embedded in constitutive H3K9me3 heterochromatin [35, 45]. Altogether, these reports present a contrasted picture, suggesting that connections between heterochromatin features and regulation of gene expression need to be further investigated in fungi.

In this study, we used the model system P. anserina to establish the genome-wide distribution of H3K4me3, H3K9me3 and H3K27me3. Most of the H3K9me3 is found in repeat-rich regions also enriched in H3K27me3, while both the H3K9me3 and H3K27me3 are co-exclusive with the H3K4me3 marks. We demonstrated that absence of PaKmt1, the H3K9me3 methyltransferase, resulted in loss of H3K9me3 and in significant reduction of H3K27me3, but had a limited impact on P. anserina growth, differentiation and sexual reproduction. Likewise, removal of PaHP1 diminished H3K9me3 and H3K27me3 modifications, with impact on P. anserina physiology similar to that of the PaKmt1 mutation. Conversely, we observed that absence of PaKmt6, the H3K27me3 methyltransferase, resulted in drastic depletion of both H3K27me3 and H3K9me3 and caused severe defects in most aspects of the life cycle including growth, differentiation processes and sexual reproduction. These findings point towards an interdependence of the two types of heterochromatin within the repeated regions of the genome as well as a conserved function of Polycomb silencing in P. anserina development.

Results

Chromatin-modifying enzymes are present in the P. anserina genome and expressed throughout its life cycle

The P. anserina genome contains a single gene (Pa_6_990) encoding a putative H3K9me3 methyltransferase, named PaKmt1 (342 aa) and a single gene (Pa_1_6940) encoding a putative H3K27me3 methyltransferase, named PaKmt6 (1090 aa). Both of these highly conserved proteins show a canonical SET domain typical of bona fide histone methyltransferases (Fig. 1a and Additional file 1: Fig. S1A–C). PaKmt1 displays pre- and post-SET domains (PS50867 and PS50868) (Additional file 1: Fig. S1D). PaKmt6 contains a CXC domain (PS51633), a cysteine-rich motif frequently found in polycomb group (Pc-G) proteins, but no chromodomain. PaKmt1 and PaKmt6 cluster in two different clades of histone methyltransferase (Fig. 1b). Their evolutionary history is canonical since the topology of the two branches of the tree is consistent with the species tree.

Fig. 1
figure1

Structure and functions of the histone methyltransferases PaKmt1 and PaKmt6. a Domain structure of histone methyltransferases PaKmt1 and PaKmt6. Sizes in amino acid (aa) are given (right). PRE-SET (red, IPR007728), SET (orange, IPR001214) and POST-SET (green, IPR003616) conserved domains are required for H3K9 methyltransferase activity of KMT1 homolog proteins. CXC (yellow, IPR026489) is a cysteine-rich conserved domain located in the H3K27 methyltransferase catalytic domain of KMT6 homologs. b Phylogenetic analysis of KMT1 and KMT6 histone methyltransferase homologs. This tree regroups H3K9 (blue) and H3K27 (red) methyltransferase proteins from fungi, plants and metazoans. They cluster in two different clades of histone methyltransferases. Their evolutionary history is canonical since the topology of the two branches of the tree is consistent with the species tree. Bootstraps are given. Filamentous fungi: Podospora anserina (Pa), Neurospora crassa (Nc), Fusarium graminearum (Fg), Trichoderma reesei (Tr), Epichloë festucae (Ef) Botrytis cinerea (Bc), Magnaporthe oryzae (Mo), Zymoseptoria tritici (Zt), Leptosphaeria maculans (Lm), Aspergillus nidulans (An), Aspergillus fumigatus (Af), Penicillium oxalicum (Po), Ascobolus immersus (Ai), Puccinia graminis (Pg), Pneumocystis jirovecii (Pj), Conidiobolus coronatus (Cc), Mucor circinelloides (Mc); yeasts: Schizosaccharomyces pombe (Sp) and Cryptococcus neoformans (Cn); the worm Caenorhabditis elegans (Ce), the fruit-fly Drosophila melanogaster (Dm), the mouse Mus musculus (Mm) and the model plant Arabidopsis thaliana (At). Accession numbers for proteins used in alignments are listed in Additional file 24: Table S6. c Detection of histone post-translational modifications Western-immunoblot analysis of H3K9me3, H3K27me3 and H3K4me3 modifications in wild-type (WT); single mutants ΔPaKmt1, ΔPaKmt6; double mutant ΔPaKmt1ΔPaKmt6; complemented strains ΔPaKmt1-PaKmt1+ (PaKmt1, PaKmt1-GFP-HA), ΔPaKmt6-PaKmt6+ (PaKmt6, PaKmt6-HA). Nuclear protein samples were extracted from protoplasts. A total of 35 µg of nuclear protein was loaded in each lane. Antibodies were directed against native H3 histone (H3), H3K9me3, H3K27me3 and H3K4me3 modifications. We detected some marginal cross-reactivity between the anti-H3K9me3 antibody and the H3K27me3 modification (evaluated at 3%, see Additional file 7: Fig. S7). Expected size of histone 3 is 17 kDa (arrows indicated 15 kDa)

RNA-seq experiments detected PaKmt6 transcripts of predicted structure in both vegetative mycelium and fruiting bodies [49]. However, while annotation of the PaKmt1 gene predicted that its ORF was made of two exons, RNA-seq experiments showed ambiguous profiles for PaKmt1 transcripts [49]. Dedicated reverse-transcription polymerase chain reaction (RT-PCR) experiments revealed two distinct transcripts (Additional file 2: Fig. S2). From 1-day-old mycelium cDNA, we amplified a single 1108-bp fragment corresponding to the PaKmt1 spliced transcript. However, from either 4-day-old mycelium RNA extracts or fruiting body RNA extracts (2-day and 4-day post-fertilization fruiting bodies), we amplified an additional 1170-bp fragment, which corresponds to a PaKmt1 unspliced transcript. The translation of the 1170-bp fragment results in a truncated protein of 76 amino acids lacking the conserved functional motifs (pre-SET and SET domains). This observation suggests that, in addition to the functional full length PaKmt1 H3K9 methyltransferase, P. anserina may produce a truncated version of this protein at specific developmental stages of its life cycle.

Genome-wide distribution of H3K4me3, H3K9me3 and H3K27me3 histone modifications in the P. anserina genome

Using western-blotting, we first showed that H3K4me3, H3K9me3 and H3K27me3 marks were present in the P. anserina genome (Fig. 1c, WT). Then, to further characterize the genomic patterns of these three histone modifications, we performed ChIP-seq experiments on chromatin samples collected from vegetative mycelia grown 60 h at 27 °C. For mapping genome-wide patterns of histone modifications, we used two distinct methodologies (i) domainogram analysis [50] based on enrichment within multiscale windows and (ii) peak calling analysis using the MACS2 algorithm (Fig. 2a and Additional file 19: Table S1). Clustering analysis of histone marks in the wild-type background established that the biological replicates clustered together with high correlation coefficients (Additional file 3: Fig. S3A).

Fig. 2
figure2

Distribution of histone marks in wild-type P. anserina. a Wild-type distribution of H3K4me3, H3K9me3 and H3K27me3 on chromosomes one and five. Domainograms (top) show significance of enrichment of H3K4me3, H3K9me3, H3K27me3 marks in windows of varying size. Color-coding of p-value is indicated (right). ChIP-seq patterns (bottom) display histone modification coverage and MACS2 detected peaks. Both patterns are represented with the same scale on chromosomes one and five of P. anserina. Both centromeres are depicted (pink). The “Mat region” (orange) corresponds to the 800 kbp region devoid of recombination and containing the mating-type locus [93]. Repeated sequences are depicted in light blue. b Peak size distribution in the wild-type P. anserina. Box plots showing sizes (in base pairs) of the MACS2-predicted peak of H3K4me3 (green), H3K9me3 (red) and H3K27me3 (blue). The box is delimited, bottom to top, by the first quartile, the median and the third quartile. Numbers above each box represent the number of peaks detected. Outliers are not shown. c Repartition of histone marks on annotated genome features in wild-type P. anserina. The pie charts represent the proportion of peaks that overlap the genomic features indicated in the legend panel

The H3K4me3 modification (Fig. 2a, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 19: Table S1, Additional file 20: Table S2 and Additional file 21: Table S3) covered nearly 11% of the P. anserina genome, with 3050 detected peaks. This mark appeared to be present in blocks of approximately 1.3 kb (Fig. 2b) in gene-rich regions, mostly located in coding sequences (CDS, Figs. 2c and 3a and Additional file 20: Table S2) and in 5′ UTR (Figs. 2c and 3a). As in other eukaryotes, H3K4me3 showed a positive correlation with active transcriptional activity (Fig. 3b). On rare instances, when H3K4me3 peaks appeared to be localized in repeats (Figs. 2c and 3a, Additional file 5: Fig. S5A, B and Additional file 21: Table S3), it mainly corresponded to transposable element (TE) relics or transposase genes inserted into coding sequences (49 out of the 63 H3K4me3 peaks found in repeats).

Fig. 3
figure3

Presence of histone marks in promoters, genes and TEs in the wild-type strain. a Metaplots. Top panel: plots of averaged ChIP-seq signal. Bottom panel: heatmaps divided in K-means built clusters representing the association versus non-association of the indicated histone modifications with the specific genomic regions. Promoter regions were defined as 1 kbp 5′ of the start codon ± 0.5 kbp of surrounding sequence (N = 10,835 Additional file 20: Table S2); coding sequences or CDS were aligned by their two ends (indicated by START and STOP) ± 1 kbp of surrounding sequence (N = 10,839 Additional file 20: Table S2); repeats were defined as TE body and duplications and the rDNA array ± 0.2 kbp surrounding each region (N = 1680 Additional file 21: Table S3). Histone modification levels in the heatmaps were calculated for non-overlapping 10 bp windows within the specific genomic regions and sorted by average value of each row. b Violin plots of gene expression according to the histone mark. Gene expression was inferred from the Transcripts Per Kilobase Million (TPM) values calculated in [49]. Log2 of the TPM values was plotted according to the marks detected for the CDS. The three lines of the violin plot represent the first quartile, the median and the third quartile. The white dot and the number above represent the mean. The numbers on top represent the number of genes in each category. K4: H3K4me3, K9: H3K9me3, K27: H3K27me3. C. Euler diagram showing overlap between the different histone marks. H3K4me3, H3K9me3 and H3K27me3 are represented by the green, red and blue areas, respectively. The numbers indicate the number of base pairs covered by each mark, according to MACS2-predicted peaks. Overlap between areas shows the number of base pairs covered by two marks

Conversely, H3K9me3 (Fig. 2a, Additional file 3: Fig. S3B–D and Additional file 19: Table S1, Additional file 20: Table S2 and Additional file 21: Table S3) was the least abundant mark found in P. anserina genome (1090 detected peaks, 2.94% of the genome), forming approximately 0.9-kb-long blocks (Fig. 2b). This mark is mostly found at repeat-rich regions, i.e., TE (728 peaks, Additional file 21: Table S3), centromeric regions and telomeric regions (Figs. 2c and 3a), so that the abundance of H3K9me3 modification in the P. anserina genome reflected its repeat content, omitting the rDNA cluster [51]. In accordance with its genomic distribution, the H3K9me3 modification showed a negative correlation with active transcriptional activity (Fig. 3b). Annotated TE families (i.e., Copia, Gypsy TEs, Tc1 mariner TEs and soloLTR relics) were marked with H3K9me3, and transcriptionally silenced, with the exception of MITEs (Additional file 5: Fig. S5A–C, Additional file 6: Fig. S6 and Additional file 21: Table S3). In the P. anserina genome, nearly all copies of TEs differ by polymorphisms generated by the repeat-induced point mutation (RIP) mechanism [51]. Originally described in N. crassa [52], RIP is a mutagenic process that occurs during the sexual dikaryotic stage of many filamentous fungi (for review see [53]). In addition of C–G to T–A transition accumulations, the RIPed repeats are embedded in H3K9me3-rich heterochromatin, which leads to transcriptional gene silencing. MITEs are an exception to this pattern because their short length (< 500 pb, Additional file 21: Table S3) makes them less prone to RIP.

The second most common mark found in the P. anserina genome in terms of peak numbers (1941 detected peaks, Additional file 19: Table S1), H3K27me3 nonetheless covered 18% of the genome (Fig. 2a, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 21: Table S3), showing large blocks approximately 3.5 kb long (Fig. 2b). This mark was found enriched at sub-telomeric regions (Figs. 2c and 3a, Additional file 3: Fig. S3C, D and Additional file 4: Fig. S4), but also on about one third of the annotated coding sequences (2992 detected peaks, Additional file 20: Table S2). Most of these H3K27me3-marked coding sequences were likely not expressed during the vegetative growth phase (Fig. 3b). Importantly, H3K27me3 modification was also found on all TE families, including members of MITEs, soloLTR relics and unclassified TEs (Additional file 5: Fig. S5A, B, Additional file 6: Fig. S6 and Additional file 21: Table S3).

Zooming in on the P. anserina genomic patterns, it appeared that distribution of H3K4me3 was found to be mutually exclusive with the distributions of both H3K9me3 and H3K27me3 (Fig. 3c, Additional file 3: Fig. S3C, D, and Additional file 4: Fig. S4). This observation was further confirmed by both the domainograms, where domains with most significant enrichment of the H3K9me3 and H3K27me3 marks, but not the H3K4me3 mark, showed a large degree of overlap (Fig. 2a and Additional file 20: Table S2) and the Spearman correlation test of individual bins, where the H3K9me3 and H3K27me3 IPs were grouped in the same cluster and were negatively correlated with the H3K4me3 cluster (Additional file 3: Fig. S3A). On the same note, we noticed that chromosome five is significantly enriched in H3K9me3 and depleted in H3K4me3 (chi2 test, p-value < 2.2e−16), which is in line with its high TE and repeat content. However, we did find a few genes that harbored both H3K4me3 and H3H27me3 marks in these experimental conditions (Additional file 20: Table S2). Remarkably, this set was significantly enriched (p-value < 0.001) in genes encoding proteins involved in non-allelic heterokaryon incompatibility (HI), transcription and chromatin remodeling (Additional file 20: Table S2). More importantly, patterns of H3K27me3 and H3K9me3 modifications, determined either as peaks or by significant genomic coverage, were often found overlapping in gene-poor compartments scattered across the genome (Fig. 3c, Additional file 3: Fig. S3C, D, Additional file 4: Fig. S4 and Additional file 6: Fig. S6), especially on repeats, which included most of the TE families (Fig. 3a: cluster_1, and Additional file 21: Table S3). Indeed, we identified 954 H3K27me3 and H3K9me3 common peaks, representing 90% of the H3K9me3 peaks (Additional file 19: Table S1 and Additional file 21: Table S3). This observation was consistent with the Spearman clustering analysis showing that H3K9me3 and H3K27me3 IP samples belong to the same cluster (Additional file 3: Fig. S3A).

PaKmt1 and PaKmt6 are histone methyltransferases and deletion of their genes impacts both H3K9me3 and H3K27me3 distribution

To determine the histone methyltransferase(s) required for the deposition of the heterochromatin marks, we constructed ΔPaKmt1 and the ΔPaKmt6 deletion mutants (Additional file 8: Fig. S8). Double ΔPaKmt1ΔPaKmt6 mutant strains were obtained by crossing the corresponding single mutants. All three mutants were viable, demonstrating that neither PaKmt1 nor PaKmt6 are essential genes. Western-blotting showed that H3K9me3 was almost absent in proteins extracted from the ΔPaKmt1 mutant (Fig. 1c), whereas the H3K27me3 was absent from the ΔPaKmt6 sample. The faint remaining signal seen in the ΔPaKmt1 sample could be due to the cross-reactivity between the anti-H3K9me3 antibody and the H3K27me3 modification, which was evaluated at 3% of the specific anti-H3K9me3/H3K9me3 immunogenic reaction (Additional file 7: Fig. S7). Importantly, neither mark was present in the double ΔPaKmt1ΔPaKmt6 mutant, showing that no additional H3K9 or H3K27 methylation exists in P. anserina. As expected, H3K4me3 was detected in all samples. We concluded that PaKmt1 is required for H3K9me3, while PaKmt6 is required for H3K27me3.

ChIP-seq experiments performed with either ΔPaKmt1 mutants or ΔPaKmt6 mutants showed that the H3K4me3 patterns remained unchanged in all the mutant genetic contexts, suggesting that loss of H3K27me3 or H3K9me3 has no impact on the H3K4 methyltransferase activity (Fig. 4a, b, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 19: Table S1 and Additional file 20: Table S2).

Fig. 4
figure4

Distribution of histone marks in ΔPaKmt6 and ΔPaKmt1 strains. a Visualization of H3K4me3, H3K9me3 and H3K27me3 localization across the genome in wild-type and mutant backgrounds. Domainograms (top) show significance of enrichment of H3K4me3, H3K9me3 and H3K27me3 marks in windows of varying size, according to the indicated p-value (right). ChIP-seq patterns (bottom) display histone modification coverage and MACS2 detected peaks. Both patterns are represented with the same scale on chromosomes one and five of P. anserina. Both centromeres are depicted (pink). The “Mat region” (orange) corresponds to the 800 kbp region devoid of recombination and containing the mating-type locus [93]. Repeated sequences are depicted in light blue. b Peak size distribution in the ΔPaKmt1 and ΔPaKmt6 mutant strains. Box plots showing sizes in base pairs of the MACS2-predicted peaks of H3K4me3, H3K9me3 and H3K27me3, shown in green, red and blue, respectively. The box is delimited, bottom to top, by the first quartile, the median and the third quartile. Numbers on top of each box represent the number of peaks detected. Outliers are not shown

According to domainogram analysis, ChIP-seq experiments performed in the ΔPaKmt1 background revealed a nearly comparable span of H3K9me3 domains, as compared to WT (Fig. 4a, b, Additional file 3: Fig. S3B–D, Additional file 20: Table S2), while MACS2 analysis did not detect any significant peak enrichment (Fig. 4a, b, Additional file 19: Table S1, Additional file 20: Table S2 and Additional file 21: Table S3 and Additional file 4: Fig. S4). Because of differences in significance calling between domainograms and MACS peak calling, combined with a low degree of cross-reactivity of the anti-H3K9me3 antibody (Additional file 7: Fig. S7), this remaining signal could be due to H3K27me3 detection in this mutant genomic context devoid of H3K9me3. This would mean that PaKmt1 (i) is responsible for the wild-type distribution of H3K9me3 modifications at the specific locations described above and (ii) is the main P. anserina H3K9me3 methyltransferase. Unexpectedly, our data revealed variations of the H3K27me3 patterns in the ΔPaKmt1 mutant background (Fig. 4a, b, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 19: Table S1, Additional file 20: Table S2 and Additional file 21: Table S3). First, the average length of the H3K27me3 blocks was reduced by two thirds (Fig. 4b). Second, about 30% of the H3K27me3 peaks were missing, especially those located in the gene-rich blocks (Fig. 4a, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 20: Table S2). Third, some the remaining H3K27me3 peaks were re-localized within TEs, where they replaced the missing H3K9me3 marks (Additional file 21: Table S3).

ChIP-seq experiments performed using the ΔPaKmt6 mutant strains revealed that this deletion drastically reduced the overall H3K27me3 content (significant genome coverage dropped from 18.95 to 0.04%, (Fig. 4a, b, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 19: Table S1, Additional file 20: Table S2 and Additional file 21: Table S3). To our surprise, we noticed that the loss of PaKmt6 severely impacted the amount of H3K9me3 marks, since in the ΔPaKmt6 background, H3K9me3 genome coverage dropped from 2.94 to 0.18% (Fig. 4a, b, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 19: Table S1 and Additional file 21: Table S3). This means that the diminution in H3K9me3 peaks was as strong as the diminution in H3K27me3 peaks in this mutant background. Nevertheless, the remaining H3K9me3 peaks in the ΔPaKmt6 background were not displaced, nor were the missing H3K9me3 blocks replaced by H3K27me3 ones (Additional file 3: Fig. S3C, D, Additional file 9: Fig. S9 and Additional file 21: Table S3).

Loss of H3K27me3 in ΔPaKmt6 background modifies gene expression

To assay the impact of H3K27me3 and/or H3K9me3 loss on gene expression, we performed RT-qPCR experiments on vegetative mycelia grown 60 h at 27 °C, for either the wild-type genotype (5 independent samples) or mutant genotypes (ΔPaKmt1 and ΔPaKmt6, 5 independent samples each). We first focused on six genes that harbor H3K27me3 marks in wild-type background (Fig. 5a, Additional file 10: Fig. S10 and Additional file 20: Table S2). In the ΔPaKmt6 background (Fig. 5a, Additional file 20: Table S2), all of them lost their H3K27me3 marks and were over-expressed (fold change ≥ 2). This result is consistent with H3K27me3 being involved in regulation of gene expression during development programs since Pa_1_6263, which showed the most dramatic up-regulation effect (fold change > 25, p-value = 0.001 Fig. 5a), is normally expressed only during the sexual phase [49]. Expression of genes encoding secondary metabolites (SM) is known to be sensitive to chromatin compaction [43, 44, 54, 55]. Indeed, when we assayed the expression of two such genes (Pa_6_7270 and Pa_6_7370) in the ΔPaKmt6 background, we found them both up-regulated (Fig. 5a, Additional file 20: Table S2). Notably, loss of PaKmt1 had no significant impact (fold change ≤ 2) on the expression of this set of genes, even if some of the H3K27me3 marks were lost in this mutant background (Additional file 10: Fig. S10 and Additional file 20: Table S2).

Fig. 5
figure5

Gene expression in the ΔPaKmt6 strain and growth features of the ΔPaKmt6 and ΔPaKmt1 strains. a Relative expression of selected genes in ΔPaKmt6 mutant strains. Upper panel: histone marks in the wild-type strain. For each selected gene, a black dot shows the presence of the mark. Conversely a white dot shows its absence. The error bars represent the 95% confidence interval. Pa_4_1170 and Pa_1_16300 encode two putative Glycoside Hydrolases from families 7 and 61, respectively [51]. Pa_6_7270 and Pa_6_7370 belong to a secondary metabolite gene cluster located on a chromosome 7 sub-telomeric region. Pa_5_10 encodes the meiotic drive element Spok2 located on Chromosome 5 [56]. Pa_1_1880 encodes the putative DNA repair RAD50 protein. Pa_7_9210 encodes a putative exonuclease. They are both readily expressed at all stages of the life cycle and their expression remained unchanged in both of the two mutant backgrounds. Middle panel: normalized expression ratio of the selected genes relative to the wild-type strain. Overexpression (fold change ≥ 2, dark black line): Pa_1_6263: expression ratio = 25,987, p-value = 0.001; Pa_4_1170: expression ratio = 6375, p-value = 0.004; Pa_1_16300: expression ratio = 7588, p-value = 0.003; Pa_6_7270: expression ratio = 5280, p-value < 0001; Pa_6_7370: expression ratio = 6016, p-value = 0.002; Pa_5_10: expression ratio = 2212, p-value < 0.001; Pa_1_1880: expression ratio = 1.216, p-value = 0.192; Pa_7_9210: expression ratio = 0.756, p-value = 0.061. Three normalization genes (AS1, GPD and PDF2) were selected with geNorm [105] among eight housekeeping genes. Average expression stability of these three normalization genes is M = 0.230 and their pairwise variation is V3/4 = 0.071. See “Methods” for details and Additional file 22: Table S4 for Cq and selection of normalization genes and NRT-qPCR controls. Lower panel: histone marks in the ΔPaKmt6 mutant strain. They are indicated as in the upper panel (see above). b Quantification of expression of two Tc1_mariner-like members in ΔPaKmt1 and ΔPaKmt6 mutant strains. ΔPaKmt1 genetic context: no overexpression (fold change ≥ 2, dark black line). ΔPaKmt6 genetic context: Tc1_mariner-like_pelobates family: expression ratio = 7.4, p-value = 0.009 in ΔPaKmt6 genetic context and of Tc1_mariner-like_rainette family: expression ratio = 2.005, p-value = 0.004 in ΔPaKmt6 genetic context. Normalization as described in a. c Phenotypic characterization of ΔPaKmt1 and ΔPaKmt6 single mutant strains. The mat+ strains of the indicated genotypes were grown on M2 minimal medium for 5 days at 27 °C. d Vegetative growth kinetics of ΔPaKmt1 and ΔPaKmt6 single mutant strains on M2 minimal medium. See “Methods” section for details

We also tested the expression of three TE families (Fig. 5b) in the wild-type background. The Copia_Ty1_nephelobates and the Tc1_mariner-like_pelobates are heavily RIPed and harbor both H3K9me3 and H3K27me3 marks (Additional file 21: Table S3 and Additional file 5: Fig. S5A, B). In consequence, their respective expression was silenced (Additional file 5: Fig. S5C and Additional file 22: Table S4). In contrast, the Tc1_mariner-like_rainette family shows fewer RIP mutations, along with H3K4me3 modifications but neither H3K9me3 nor H3K27me3 marks (Additional file 21: Table S3) and therefore was expressed (Additional file 22: Table S4). In the ΔPaKmt6 background, expression of the Tc1_mariner-like_pelobates was significantly increased (fold change = 7.4, p-value = 0.009 Fig. 5b), while RT-qPCR was inconclusive for Copia_Ty1_nephelobates (i.e., RT-qPCR not different from NRT-qPCR, see “Methods” section). As expected, expression of the Tc1_mariner-like_rainette family was even further increased (fold change = 2.005 p-value = 0.004 Fig. 5b).

Loss of PaKmt6 has stronger impact on vegetative growth than loss of PaKmt1 and results in male gamete over-production

We then assessed the impact of loss of either H3K27me3 or H3K9me3, or both, on the P. anserina life cycle.

Pigmentation, branching and the overall aspect of ΔPaKmt1 mycelia were similar to those of wild-type strains, except for reduced production of aerial hyphae (Fig. 5c) and slight growth defects (Fig. 5d and Table 1), including an unusual failure to resume growth after it has been stopped for a few days (Additional file 11: Fig. S11A). We ruled out the possibility that these phenotypes result from activation of Crippled Growth (CG), an epigenetically triggered cell degeneration [57], by showing that the ΔPaKmt1 mutation did not promote CG (Additional file 11: Fig. S11B) but behaved like the wild-type strain in that respect. This suggests that the reduced capability to resume growth may be indicative of less plasticity to perform an adaptive developmental program under changing environmental conditions.

Table 1 Mycelium growth at different temperatures

Pinkish pigmentation, crooked branching pattern and the overall aspects of ΔPaKmt6 mycelia were clearly different from those of wild-type strains (Fig. 5c). First of all, ΔPaKmt6 mutant growth rate was 43% lower than that of wild-type strains during the first 60 h of growth on M2 at 27 °C (Fig. 5d). After this time, progression of the ΔPaKmt6 colony margins became so erratic and irregular that the standard assay consisting in thallus radius measurements was no longer accurate. Indeed, ΔPaKmt6 mycelia showed “stop-restart” growth resulting in patchy mycelium sectors, which literally sprouted out of the colony margins like cauliflower inflorescences (Fig. 5c). These outgrowths are manifest in the increasingly large error bars along the corresponding growth curve (Fig. 5d). However, we were able to grow the ΔPaKmt6 mutants in race tubes to determine that this mutation had no effect on senescence (Additional file 12: Fig. S12).

In addition to growth defects, the ΔPaKmt6 mutant mycelium showed a powdery phenotype. At the microscopic scale, we observed that this corresponds to an over-production of male gametes (spermatia) (app. 130 times more than in the wild-type strains, Fig. 6a). We tested the ability of these ΔPaKmt6 male gametes to fertilize wild-type female gametes (ascogonia). To do so, we sprayed a defined quantity of either mutant or wild-type male gametes onto wild-type female gametes of compatible mating type. Given that one fruiting body originates from one fertilization event, we established that the ΔPaKmt6 male gametes were as efficient as the wild-type ones, since a similar number of fruiting bodies was obtained in each condition (4 independent experiments). This functional test indicates that the ΔPaKmt6 mutant male gamete over-production had no obvious impact on their biological properties.

Fig. 6
figure6

Sexual features of the ΔPaKmt6 strain. a Spermatium production in wild-type and mutant backgrounds. See “Methods” section for details. b ΔPaKmt6 homozygote crosses result in non-homogeneous distribution of late appearing fruiting bodies showing multiple morphological defects. These crippled structures could be due to lack of neck, lack of ostiols, lack of setae, or non-canonical pear-shape envelopes, supernumerary necks, misdirected orientation, upside down growth in agar, fusions of fruiting bodies, etc.

Finally, when we constructed the ΔPaKmt1ΔPaKmt6 double mutant, the phenotypes generated by the ΔPaKmt6 null allele were epistatic upon those generated by the ΔPaKmt1 allele.

The ΔPaKmt6 mutant formed crippled fruiting bodies yielding a reduced number of ascospores

We then investigated the ability of the histone methyltransferase mutant strains to perform sexual reproduction. Sexual development and ascospore production of crosses involving two ΔPaKmt1 parental strains were indistinguishable from those of wild-type strains. By contrast, when we performed homozygous ΔPaKmt6 crosses, fruiting bodies were formed, which means that fertilization occurred, but their developmental time frame was delayed. They also displayed a large palette of morphological defects (e.g., lack of necks or extra necks, lack of ostioles, upward/downward mis-orientation, development into agar, fusion) and their ascospore production was significantly reduced (Fig. 6b).

Heterozygous crosses showed that when the ΔPaKmt6 null allele was present in the female gamete genome and the wild-type allele was present in the male gamete genome, the fruiting bodies were crippled and the ascospore production reduced, as in the homozygous ΔPaKmt6 crosses (Additional file 13: Fig. S13A). Conversely, heterozygous crosses involving female wild-type strains and male mutant strains produced well-formed and fully fertile fruiting bodies, indicating a ΔPaKmt6 maternal effect.

To determine whether the ΔPaKmt6 mutant’s reduced fertility resulted from either a fruiting body envelope defect or a zygotic tissue defect, we set up tricaryon crosses involving the Δmat strain [58, 59]. The Δmat strain lacks the genes required for fertilization, so it does not participate either as male or female in the sexual reproduction. However, the Δmat mycelium is able to provide maternal hyphae required to build fruiting bodies. Consequently, the Δmat strain can complement mutants defective for fruiting body tissue formation but not for zygotic tissue defect. In this study, we observed that the Δmat; ΔPaKmt6 mat+; ΔPaKmt6 mat− tricaryon yielded few but fully matured and well-formed fruiting bodies that produced normal ascospores (Additional file 13: Fig. S13B). This experiment demonstrates that ΔPaKmt6 mutation compromised the fruiting body envelope formation but had no impact on zygotic development. This result is in line with the observed ΔPaKmt6 maternal effect (see the above paragraph).

Finally, crosses involving the double ΔPaKmt1ΔPaKmt6 mutant strains displayed a typical ΔPaKmt6 phenotype with respect to the perithecium number, morphology and ascospore production, demonstrating that the ΔPaKmt6 null allele was epistatic. In addition, one third of the asci recovered from homozygous ΔPaKmt1ΔPaKmt6 crosses showed fewer pigmented ascospores than asci from wild-type crosses (with a color panel spanning from white to green). In general, deficit in pigmentation is indicative of incomplete maturation during the ascospore formation process [60].

The PaKmt6 null allele reduces the ascospore germination efficiency

Homokaryotic ascospores originating from homozygous ΔPaKmt1 crosses germinated like those of wild-type crosses (> 96% germination rate, N = 100, Table 2). However, about two thirds of normal looking ascospores produced by the homozygous ΔPaKmt6 crosses did not germinate (3 independent experiments, N = 100 for each experiment, Table 2). The reduced germination efficiency was further enhanced in ΔPaKmt1ΔPaKmt6 double mutants (Table 2). On germination medium, wild-type ascospore germination occurs at a germ pore located at the anterior side of the ascospore, by the extrusion of a spherical structure. This transient structure, called a germination peg, gives rise to a polarized hypha that develops into mycelium. In the ΔPaKmt6 mutants, the ascospore germination process was stopped at an early stage, since germination pegs were not formed. However, the ΔPaKmt6 defect was ascospore autonomous since ascospores from wild-type progeny obtained in WT × ΔPaKmt6 crosses germinated normally, whereas those from ΔPaKmt6 progeny did not. In addition, PaKmt6+PaKmt6 dikaryotic ascospores had a wild-type germination rate showing that the ΔPaKmt6 defect was recessive (2 independent experiments, N = 80 dikaryotic ascospores).

Table 2 Ascospore germination efficiency

It was previously shown that unmelanized ascospores carrying the pks1-193 mutant allele germinate spontaneously at high frequency (over 95%), even on non-inducing media [60, 61]. Crosses involving pks1-193 ΔPaKmt6 mat- and pks1+ ΔPaKmt6 mat+ strains yielded asci made of two unmelanized (pks1-193 ΔPaKmt6) and two melanized (pks1+ ΔPaKmt6) ascospores. As expected, about one third of the ΔPaKmt6 melanized ascospores germinated on inducing medium. By contrast, almost all of the ΔPaKmt6 unmelanized spores germinated, even on non-inducing medium. This result shows that the absence of melanin suppressed the ΔPaKmt6 ascospore germination defect, as found for ΔPaPls1 and ΔPaNox2 mutants [60].

Loss of PaHP1 leads to overall reduction and partial relocation of H3K9me3 modifications

We searched the P. anserina genome for a homolog of heterochromatin protein 1 (HP1, [62]). This protein is known not only to bind the H3K9me3 marks of the constitutive heterochromatin, but also to perform a wide variety of functions ranging from cohesion of sister chromatids for the maintenance of telomeres, DNA repair and replication, cell cycle control and even RNA splicing [63]. We identified one gene (Pa_4_7200) encoding the protein PaHP1 (252 aa, 57% identity, e-value = 4e−78) (Additional file 14: Fig. S14A). It displays conserved chromo (PS00598) and chromo-shadow (PS50013) domains (Additional file 14: Fig. S14A, B and Additional file 15: Fig. S15A). By self-aggregation, chromo-shadow domain-containing proteins can bring together nucleosomes and, thus, condense the corresponding chromatin region [13]. RNA-seq data showed that the PaHP1 gene is expressed in both vegetative mycelium and fruiting body tissues [49]. We also explored PaHP1-GFP cellular localization in 2-day-old growing mycelium. In wild-type PaHP1-GFP-expressing strains, the fluorescence signal was nuclear and punctate (Fig. 7a). However, when the PaHP1-GFP allele was expressed in a ΔPaKmt1 background, the GFP fluorescence was still nuclear but no longer punctate (Fig. 7a), as already observed in N. crassa [48]. On the contrary, when expressed in a ΔPaKmt6 background, the PaHP1-GFP signal appeared punctate as in wild-type nuclei.

Fig. 7
figure7

Chromatin modifications in the ΔPaHP1 strain. a Subcellular localization of PaHP1. Fluorescence microscopy pictures show PaHP1-GFP chimeric protein during vegetative growth in mycelium. Nuclei are marked using the PaH1-mCherry allele that encodes histone H1 protein tagged with mCherry fluorophore (red). Indicated scale is the same for all the pictures (5 µm). PaHP1-GFP displays nuclear localization and forms foci that likely correspond to heterochromatin domains (arrows). This specific localization is disturbed only in ΔPaKmt1 mutants. b Visualization of histone mark localization across the genome in the ΔPaHP1 strain. Domainograms (top) show significance of enrichment of H3K4me3, H3K9me3 and H3K27me3 marks in windows of varying size, according to the indicated p-value (right). ChIP-seq patterns (bottom) display histone modification coverage and MACS2 detected peaks. Both patterns are represented with the same scale on chromosomes one and five of P. anserina. Both centromeres are depicted (pink). The “Mat region” (orange) corresponds to the 800 kbp region devoid of recombination and containing the mating-type locus [93]. Repeated sequences are depicted in light blue. c Peak size distribution in the ΔPaHP1 mutant strain. Box plots showing sizes in base pair of the MACS2-predicted peak of H3K4, H3K9 and H3K27 trimethylation are shown in green, red and blue, respectively. The box is delimited, bottom to top, by the first quartile, the median and the third quartile. Numbers on top of each box represents the number of peaks detected. Outliers are not shown

To explore potential role(s) of PaHP1, we deleted the corresponding Pa_4_7200 gene (Additional file 15: Fig. S15B, C). ChIP-seq experiments performed in the ΔPaHP1 background showed that the H3K4me3 pattern was similar to that of wild-type strains (Fig. 7b, c, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 19: Table S1 and Additional file 20: Table S2), while both H3K9me3 and H3K27me3 marks were reduced. The H3K9me3 genome coverage was reduced by half (Fig. 7c, Additional file 3: Fig. S3B–D, Additional file 4: Fig. S4 and Additional file 19: Table S1, Additional file 20: Table S2 and Additional file 21: Table S3), while the H3K27me3 loss of coverage was similar to that of ΔPaKmt1 strains. Interestingly, on TEs, we detected replacement of some of the missing H3K27me3 by H3K9me3 modifications. This was especially noticeable for the soloLTR relics, although this was also the case for all the annotated TE families, albeit to a lesser extent (Additional file 21: Table S3). In addition, we observed both focalized spreads of H3K9me3 marks on TEs (Additional file 3: Fig. S3C, D, Additional file 21: Table S3) and massive H3K9 methylation on the rDNA cluster (Additional file 16: Fig. S16). These observations may account for the H3K9me3 block widening in the ΔPaKmt1 genetic context. Altogether these results suggest that, in the absence of the recruitment platform mediated by the binding of PaHP1 to the constitutive heterochromatin, the compaction may be partly relaxed and the histone modifications lost, but also that the boundaries between heterochromatin and euchromatin may be less clear-cut. However, in this mutant context, the H3K9me3 reduction has no impact on gene expression (Additional file 17: Fig. S17A).

Loss of PaHP1 has no more impact than that of PaKmt1 on the P. anserina life cycle

Throughout all the P. anserina life cycle, the ΔPaHP1 mutants behave similarly to the ΔPaKmt1 ones (Table 1, Additional file 17: Fig. S17B–D). Importantly, the double ΔPaKmt1ΔPaHP1 mutants did not show any additional or more severe phenotypes (Table 1, Additional file 17: Fig. S17A). When the ΔPaHP1 allele was associated with the ΔPaKmt6 one, all the phenotypes displayed by the single ΔPaKmt6 mutants were epistatic upon those of the ΔPaHP1 mutants. The only exception is related to the germination efficiency of the triple ΔPaKmt1ΔPaHP1ΔPaKmt6, which was similar to that of the double ΔPaHP1ΔPaKmt6 rather than the double ΔPaKmt1ΔPaKmt6 mutant (Table 2). This may be due to a slightly reduced pigmentation of the ascospores produced by the homozygous ΔPaKmt1ΔPaHP1ΔPaKmt6 crosses. Indeed, even if only the darker wild-type looking ones were used in this test, slightly reduced pigmentation facilitates germination.

Discussion

In eukaryotes, the deposition of histone modifications is essential to chromatin features, which impact cell fate through regulation of gene expression. In this study, we aimed at understanding how heterochromatic territories are assembled within the P. anserina genome and how this genome organization influences the P. anserina life cycle.

H3K27me3 and H3K9me3 modifications are not mutually exclusive

To get the first hints of P. anserina’s genome organization, we took advantage of ChIP-seq experiments to uncover genome-wide chromatin landscapes in a wild-type background. Our results highlighted two uncommon features of both the constitutive and the facultative heterochromatin. Firstly, in most animals and plants, the H3K9me3 modification covers fairly large domains, while the H3K27me3 modification forms shorter ones. But for the P. anserina genome, we observed the opposite, since H3K27me3 blocks were on average three times longer than H3K9me3 blocks. Second, if as expected the H3K9me3 modification was restricted to repeat-rich regions, mostly located near telomeres and centromeres [64], the H3K27me3 modification encompasses both annotated CDS and H3K9me3 marked repeats. Co-occurrence of H3K9me3 and H3K27me3 at same genomic loci is of particular interest because these two heterochromatic marks were long considered mutually exclusive in flowering plants, mammals, and fungi [30, 65]. Yet recently, this dogma has been challenged. Mass spectrometry experiments demonstrated that H3K27me3 and H3K9me3 can coexist within the same histone octamer in embryonic stem cells [35, 66]. Co-occurrence of H3K9me3 and H3K27me3 was found at TE loci of the ciliate Paramecium tetraurelia [67], the nematode Caenorhabditis elegans [68] and the bryophyte Marchantia polymorpha [69]. In mammals, these two marks were found together in rare instances, along with DNA methylation, either at a subset of developmentally regulated genes in mouse extra-embryonic early embryo lineages [70] or in some human cancer cells [71,72,73]. Likewise, such H3K9me3 and H3K27me3 overlap, exclusively associated with DNA methylation, has also been described for A. thaliana TE loci [74].

In the basidiomycete yeast C. neoformans, removal of the Enhancer-of-zeste-like protein EZH eliminates H3K27me3, while loss of its reader, the chromodomain Ccc1 protein, redistributes H3K27me3 to new genomic locations that coincide with H3K9me2 heterochromatin domains [41]. Further disruption of C. neoformans H3K9 methyltransferase Clr4 causes loss of both H3K9me2 and redistributed H3K27me3 modifications [41]. These findings suggest that co-factors, such as readers, are responsible for anchoring the H3K27 methyltransferase to its specific facultative heterochromatin compartments, that otherwise can be recruited to H3K9me constitutive heterochromatin, resulting in H3K9me/H3K27me co-occurrence. P. anserina has neither Ccc1 homologs nor a PRC1 complex. If such recognition modules were lost, it would explain why PaKmt6 is recruited by both facultative and constitutive heterochromatin. Interestingly, H3K9me3 and H3K27me3 modifications also overlap in the sub-telomeric regions of some other filamentous fungi, such as Z. tritici, Fusarium oxysporum and Leptosphaeria maculans ‘brassicae’ [32, 33, 38, 75,76,77].

Are bivalent-like genes present in fungi?

In this study, we showed that H3K4me3 modification does not overlap with either of the H3 repressive marks at a genome-wide scale. However, we identified 50 genes (Additional file 20: Table S2) that make an exception to this rule by harboring both H3K4me3 and H3H27me3 marks. In mammals, this pattern is typical of bivalent genes found in totipotent embryonic stem cells [78], which play a central role in developmental and lineage specification processes. Unlike those of animals, fungal nuclei remain totipotent throughout the entire life cycle [79, 80]. Yet the molecular basis of this peculiar feature is not understood and it will be of interest to explore the function of the H3K4me3 and H3K27me3 dual-marked genes identified in the P. anserina genome (Additional file 20: Table S2). So far, we noted that a sizeable portion of them belongs to the het gene family. In filamentous fungi, the HET proteins trigger a cell death reaction as a consequence of a non-self-recognition process known as heterokaryon incompatibility [81]. It has been proposed that these het genes could be components of a fungal innate immune system. Since immune response must be quick and massive upon non-self-contacts, marking these genes with both repressive and activator modifications may be a way to fine tune their expression. Recently, a large number of genes associated with both H3K27me3 and H3K4me2 were identified in L. maculans ‘brassicae’ [77], indicating that bivalent domains occur in other fungi.

PaKmt1 and PaKmt6 are not independent chromatin-modifying enzymes

ChIP-seq experiments carried out on mutants allowed us to test for potential cross-talk between constitutive and facultative heterochromatin. Our data established that P. anserina’s constitutive and facultative heterochromatin are inter-dependent since removal of either PaKmt1 or PaKmt6 impacted both marks although not to the same extent. This is in contrast with the corresponding mutant phenotypes reported so far in N. crassa [33, 35], Z. tritici [45] and C. neoformans [41] where H3K9me3 content is not affected by loss of H3K27me3 and conversely (Table 3). In addition, H3K27me3 reduction in ΔPaKmt1 was not associated with any substantial genomic redistribution in contrast to what has been described in the corresponding N. crassa and Z. tritici mutants [33, 35, 45] (Table 3). Lack of PaHP1 further supported the observation that either type of heterochromatin requires the presence of the other, since in this mutant background, we saw fewer H3K9me3 and H3K27me3 marks but also reduced size of facultative heterochromatin domains. These concomitant observations suggest that initiation and propagation of the deposition of chromatin marks by the PRC2 complex are inter-dependent processes in P. anserina. Alternatively, if, as shown in A. thaliana for the FLC locus [89], that initiation and propagation of the deposition of chromatin marks by the PRC2 complex were independent in P. anserina, both of them would be impaired when the H3K9me3–heterochromatin pathway is damaged.

Table 3 Phenotypes associated with either loss of KMT1, KMT6 or HP1 in fungal species

In general, two distinct enzymes catalyze methylation of either H3K9 or H3K27 residues. However, the ciliate P. tetraurelia is endowed with only one class of histone methlytransferase, the Enhancer-of-zeste-like protein Ezl1, which methylates both H3K9 and H3K27 residues [67]. This observation suggests that eukaryotes initially acquired a single, versatile, methyltransferase that duplicated and became sub-functionalized in the course of evolution. Thus many taxa have two enzymes, with distinct substrate specificities. This ancestral versatility might be partially restored (or unmasked) by the lack of competition when one of these enzymes is missing as in the ΔPaKmt1 and ΔPaKmt6 mutant genetic contexts. Under this hypothesis, we can propose that (i) the P. anserina SU(VAR)3-9-like homolog (PaKmt1) may have a residual H3K27 histone-methyltransferase activity, that would account for the remaining H3K27me3 marks in the ΔPaKmt6 mutants, whereas (ii) SET-domain Enhancer-of-zeste-like PaKmt6 cannot be a catalytic substitute for PaKmt1, since H3K9me3 peaks were no longer detected in the ΔPaKmt1 mutants.

Constitutive chromatin is dispensable whereas facultative chromatin is required for most P. anserina developmental programs

Assessing the phenotype of ΔPaKmt1, ΔPaHP1, and ΔPaKmt6 single mutants, we established that the loss of both PaKmt1 and PaHP1 has almost no effect on the P. anserina life cycle, whereas the absence of PaKmt6 drastically impairs every developmental step. To recapitulate, the ΔPaKmt6 inactivation resulted in (i) slow, irregular, and thermo-sensitive mycelium growth, (ii) absence of aerial hyphae, (iii) striking over-production of male gametes, (iv) reduced fertility of the female gametes and/or of fertilization efficiency, (v) formation of crippled and mis-orientated fruiting bodies, (vi) reduced yields of ascospores and (vii) reduced germination efficiency of ascospores (Table 3). This means that almost every aspect of the life cycle is impacted by the loss of most of H3K27me3 marks and consequently of the facultative heterochromatin (see below). By contrast, the single ΔPaKmt1 and ΔPaHP1 mutants showed only marginal defects on aerial hyphae production and resumption of growth after a period of arrest (Table 3). More importantly, the single ΔPaKmt1 and ΔPaHP1 and double ΔPaKmt1 ΔPaHP1 mutants are fully fertile. These phenotypes contrast with those of N. crassa, as well as most filamentous fungi studied to date, i.e., Aspergillus fumigatus, E. festucae, L. maculans, Penicillium oxalicum, and Z. tritici (Table 3). Conversely, for F. graminearum and F. fujikuroi, Enhancer-of-zeste-like enzymes are key for their developmental integrity, as in P. anserina (Table 3). F. graminearum kmt6 mutants show aberrant germination patterns, over-pigmented mycelium, restricted growth, reduced pathogenicity, over-production of several SM clusters, and sterility, while loss of Kmt1 or HP1 orthologues does not lead to any noticeable phenotypes [43]. Similarly, F. fujikuroi kmt6 mutants form stunted mycelia that produce a reduced yield of conidia [44].

Conclusions

Recent studies (Table 3) suggest that the canonical definitions of facultative and constitutive heterochromatin may be challenged in fungi. Further contributing to this complex picture, in this study we showed that besides the H3K9me3 modifications that are likely installed by RIP [90], P. anserina repeats also harbor H3K27me3 modifications. It has already been postulated that the ancestral role of Enhancer-of-zeste enzymes, along with the SET-domain SU(VAR)3–9 enzymes, may be to silence TEs, through the concerted action of both H3K27me3 and H3K9me3 [67, 69, 91]. The H3K9me3/H3K27me3 overlap we have observed could therefore be a relic of the concerted action of these two classes of histone methyltransferases. Recently, a BAH-PHD protein mediating Polycomb silencing [40] and an ISWI chromatin remodeling factor were identified in N. crassa [92]. Interestingly, the phylogenetic distribution of this BAH-PHD protein suggests that it belongs to an ancestral Polycomb repression system which still exists in Fungi and early emerging eukaryotes but no longer in animals. The search for P. anserina PRC2 associated proteins will help us to understand Polycomb silencing in the absence of canonical PRC1, as well as functional cross-talk between facultative and constitutive heterochromatin.

Methods

Strains and culture conditions

The strains used in this study derived from the “S” wild-type strain of P. anserina that was used for sequencing [51, 93]. Standard culture conditions, media and genetic methods for P. anserina have been described [94, 95]. Most recent protocols can be accessed at either [96] or http://podospora.i2bc.paris-saclay.fr. Mycelium growth is performed on M2 minimal medium, in which carbon is supplied as dextrin and nitrogen as urea. Ascospores germinate on a specific germination medium (G medium) containing ammonium acetate. The methods used for nucleic acid extraction and manipulation have been described [97, 98]. Transformations of P. anserina protoplasts were carried out as described previously [99].

Identification and deletions of the PaKmt1, PaKmt6 and PaHP1 genes

PaKmt1, PaKmt6 and PaHP1 genes were identified by searching the complete genome of P. anserina with tblastn [100], using the N. crassa proteins DIM-5 (NCU04402), SET-7 (NCU07496) and HP1 (NCU04017). We identified three genes: Pa_6_990, Pa_1_6940 and Pa_4_7200 renamed PaKmt1, PaKmt6 and PaHP1, respectively. To confirm gene annotation, PaKmt1 transcripts were amplified by RT-PCR experiments performed on total RNA extracted from growing mycelia (1-day- and 4-day-old mycelium) and developing fruiting bodies (2 and 4 days post-fertilization) using primers 5s-DIM5 and 3s-DIM5 (Additional file 23: Table S5). Functional annotation was performed using InterProScan analysis (http://www.ebi.ac.uk/interpro/search/sequence-search), Panther Classification System (http://www.pantherdb.org/panther/), PFAM (http://pfam.xfam.org/) and Prosite (http://prosite.expasy.org/).

Deletions were performed on a Δmus51::bleoR strain lacking the mus-51 subunit of the complex involved in end-joining of DNA double strand breaks as described in [60]. In this strain, DNA integration mainly proceeds through homologous recombination. Replacement of the PaKmt1 wild-type allele with hygromycin-B resistance marker generated viable mutants carrying the ΔPaKmt1 null allele. Replacement of the PaKmt6 wild-type allele with a nourseothricin resistance marker generated viable but severely impaired mutants. Replacement of the PaHP1 wild-type allele with a hygromycin-B resistance marker also generated viable mutants carrying the ΔPaHP1 null allele. All deletions were verified by Southern blot (Additional file 8: Fig. S8 and Additional file 12: Fig. S12A).

Construction of fluorescent-HA-tagged chimeric proteins PaHP1-GFP-HA, PaKmt1-mCherry-HA, PaKmt6-GFP-HA

The pAKS106 and pAKS120 plasmids were described in [101]. All the cloned inserts were sequenced before transformation.

To gain insight into PaHP1 localization, we constructed a PaHP1-GFP-HA fusion chimeric protein. To this end, the PaHP1 gene was PCR amplified from the S strain genomic DNA using primers FC3-BamH1 and FC4-GFP (Additional file 23: Table S5). The resulting amplicon harbors 520 bp of promoter sequence followed by PaHP1 CDS cleared from its stop codon. In parallel, eGFP CDS was amplified from peGFP plasmid using FC5-HP1 and FC6-HindIII primers. The two DNA segments were fused by PCR using primers FC3-BamH1 and FC6-BamH1. This chimeric DNA fragment was digested with BamHI and HindIII and cloned into the pAKS106 plasmid linearized with the same enzymes. When transformed into ΔPaHP1 mutant strain, the PaHP1-GFP-HA allele is able to restore a wild-type phenotype (growth rate and aerial mycelium), showing that the tagged version of PaHP1 protein is functional.

PaKmt1-mCherry-HA allele construction followed the same experimental design used for PaHP1-GFP-HA construction. First, PaKmt1 was PCR amplified using Dim5mChFBamH1 and Dim5mChR and resulted in a DNA amplicon composed of 454 bp of promoter followed by the PaKmt1 CDS. The mCherry CDS was amplified from the pmCherry plasmid using mChDim5F and mChDim5RHindIII primers. Fragments were fused by PCR, digested and cloned in pAKS106. The Kmt1-mCherry-HA allele was able to restore wild-type growth phenotype when transformed in ΔPaKmt1 mutant strains, indicating that this chimeric protein is functional.

PaKmt6 (Pa_1_6940) was PCR amplified from S genomic DNA using FC73-GA1KMT6 and FC74-GA1KMT6 primers, resulting in 4232 bp DNA amplicons composed of 912 bp of promoter followed by PaKmt6 CDS cleared from its stop codon. In parallel, the GFP allele was amplified from peGFP plasmid using FC75-GA1GFP and FC76-GA1GFP. Both fragments were cloned together in one step in the pAKS-106 plasmid previously digested with NotI and BamHI restriction enzymes. This cloning step was performed using Gibson Assembly kit (New England Biolabs). This yielded the chimeric PaKmt6-GFP-HA allele. An additional PaKmt6-HA allele was generated by PCR amplification from the wild-type strain DNA of PaKmt6 CDS along with its own promoter, using FC65-NotI and FC72-BamHI primers. These amplicons were digested and cloned in pAKS106 plasmid using NotI and BamHI restriction enzymes. When transformed into ΔPaKmt6 mutant strain, both PaKmt6-GFP-HA and PaKmt6-HA alleles were able to restore wild-type growth and sexual development, suggesting that this chimeric version of the PaKmt6 protein is functional. However, no GFP fluorescence signal was detected in the ΔPaKmt6, PaKmt6-GFP-HA complemented strains.

To overexpress PaHP1, the endogenous promoter were replaced by the AS4 promoter that drives high and constitutive expression throughout the life cycle [102]. PaHP1-GFP allele was PCR amplified from pAKS106-HP1-GFP (this study) using primers FC21-BamH1 and FC6-HindIII. The resulting 1877 bp amplicons harbor the chimeric construction cleared from its stop codon and promoter. The PCR product was digested with BamHI and HindIII and then cloned in pAKS120 previously digested with the same enzymes. The pAS4-PaHP1-GFP-HA resulting allele was introduced in ΔPaHP1 strains by transformation.

pAS4-PaKmt1-mCherry-HA over-expressed allele construction followed the same experimental procedures as the one used for pAS4-HP1-GFP-HA construction. PaKmt1-mCherry CDS was PCR amplified using FC22-BamHI and mChDim5R-HindIII using pAKS106-PaKmt1-mCherry-HA plasmid as DNA template. Amplified DNA fragments were then digested and cloned in pAKS120 plasmid using HindIII and BamHI restriction enzymes. This allele was transformed in ΔPaKmt1 strain.

Complementation experiments with GFP-tagged proteins

To perform complementation experiments, the PaKmt1-mCherry allele was transformed into the ΔPaKmt1 strain. Because the PaKmt1-mCherry allele was linked to a phleomycin marker, 55 phleomycin-resistant transformants were recovered. Among those, five wild-type looking independent transgenic strains were selected and crossed to a wild-type strain. In the progeny, all the phleomycin- and hygromycin-resistant haploid strains behaved similarly to the wild-type strain. Two of these independent transformants were selected for further phenotypic characterization. To ask whether a constitutively over-expressed PaKmt1 allele could have an impact on the physiology of P. anserina, we introduced the AS4-PaKmt1-mCherry allele into the ΔPaKmt1 strain. All the phleomycin-resistant transformants looked like a wild-type strain (N = 30 phleomycin resistant). Two of them were crossed to a wild-type strain. In their progeny, all the phleomycin-resistant strains carrying either the PaKmt1+ allele or the ΔPaKmt1 were indistinguishable from wild-type strains.

Similarly, we transformed a PaHP1-GFP-HA allele, linked to a phleomycin marker, into a mutant ΔPaHP1 strain. Among the 20 phleomycin-resistant transformants that were recovered, five independent transformants displaying a wild-type phenotype were selected and crossed to the wild-type strain. In the progeny, all the phleomycin- and hygromycin-resistant haploid strains behaved similarly to the wild-type strain. Transformation of the AS4-PaHP1-GFP-HA allele into the ΔPaHP1 strain generated phleomycin-resistant transformants (N = 19) showing intermediate to full complementation but no extra phenotypes. Two transformants presenting full complementation were crossed to a wild-type strain. In their progeny, all the phleomycin-resistant strains carrying either the PaHP1+ allele or the ΔPaHP1 behaved like the wild-type strains.

Finally, the PaKmt6-HA allele and the PaKmt6-GFP-HA allele were independently transformed into a ΔPaKmt6 strain. For the two transformation experiments, among the phleomycin-resistant transformants that were recovered (N = 20 and N = 32, respectively), two independent transgenic strains displaying a wild-type phenotype were selected and crossed to a wild-type strain. In the progeny of these four crosses, all the phleomycin- and nourseothricin-resistant haploid strains behaved similarly to wild-type strain.

Phylogenetic analysis

Orthologous genes were identified using the MycoCosm portal [103] and manually verified by reciprocal Best Hit BLAST analysis. Sequences were aligned using Muscle (http://www.ebi.ac.uk/Tools/msa/muscle/) and trimmed using Jalview (version 2.9.0b2) to remove non-informative regions (i.e., poorly aligned regions and/or gap containing regions). Trees were constructed with PhyML 3.0 software with default parameters and 100 bootstrapped data sets [104]. The trees were visualized with the iTOL server (http://itol.embl.de/).

Phenotypic analyses

Spermatium counting was performed as follows: each strain was grown on M2 medium with two cellophane layers at 18 °C for 14 days (P. Silar’s personal communication). To collect spermatia, cultures were washed with 1.5 mL of sterile water. A Malassez counting chamber was used for enumeration.

Cytological analysis

Light microscopy was performed on 2-day-old mycelium. Explants taken from complemented strains expressing PaHP1-GFP-HA were analyzed by fluorescence microscopy. Pictures were taken with a Leica DMIRE 2 microscope coupled with a 10-MHz Cool SNAPHQ charge-coupled device camera (Roper Instruments), with a z-optical spacing of 0.5 mm. The GFP filter was the GFP-3035B from Semrock (Ex: 472 nm/30, dichroïc: 495 nm, Em: 520 nm/35). The Metamorph software (Universal Imaging Corp.) was used to acquire z-series. Images were processed using the ImageJ software (NIH, Bethesda).

Western blot analysis

Western blots were performed on nuclear extracts made from 200 mg of ground mycelia resuspended in 3 mL of New Cell Lysis Buffer [Tris–HCL pH 7.5 50 mM, NaCl 150 mM, EDTA 5 mM, NP-40 0.5%, Triton 1% and protease inhibitor (Roche-04693132001)]. The nuclear extracts were collected after centrifugation (300×g, 5 min) and resuspended in 300 µL of Nuclei Lysis Buffer (Tris–HCl pH 8 50 mM, EDTA 10 mM, SDS 0.5% and protease inhibitor (Roche-04693132001). Nuclei were then sonicated (Diagenode bioruptor plus sonicator, 20 min at high level). Protein contents of nuclear extracts were quantified using the Qubit Protein assay kit (Invitrogen). Then, 35 µg of nuclear extract were diluted in Laemmli sample buffer and loaded on a 10% SDS polyacrylamide gel. Membranes were probed using the following high affinity ChIP-grade monoclonal antibodies: anti-Histone 3 (Abcam ab12079), anti-H3K9me3 (Abcam ab8898), anti-H3K27me3 (Millipore 17-622) and anti-H3K4me3 (Active Motif 39159).

Chromatin immuno-precipitation

Each ChIP-seq experiment was performed on two independent biological replicates for each genotype (mat+ strains only). Chromatin extractions were made from mycelium grown for 60 h at 27 °C in liquid culture (KH2PO4 0.25 g/L, K2HPO4 0.3 g/L, MgSO4 0.25 g/L, urea 0.5 g/L, biotin 0.05 mg/L, thiamine 0.05 mg/L, oligo-element 0.1%, yeast extract 5 g/L and dextrin 5.5 g/L). Mycelium was filtered and quickly washed with phosphate-buffered saline (PBS), resuspended and then cross-linked in 100 mL of PBS with 1% paraformaldehyde for 15 min at 27 °C. Reaction was stopped by addition of glycine (60 mM). Cross-linked mycelium was filtered, washed with cold PBS buffer and dried on Whatman paper. Dried mycelium was frozen in liquid nitrogen and ground to powder in a mortar. Aliquots of mycelium powder (200 mg) were resuspended in 1 mL of Lysis Buffer [Hepes pH 7.5 50 mM, NaCl 0.5 M, EDTA pH 8 1 mM, Triton X-100 1%, Sodium deoxycholate 0.1%, CaCl2 2 mM and protease inhibitor (Roche)]. Chromatin was digested to 0.15–0.6-kb fragments using Micrococcal nuclease (NEB #M0247) for 30 min at 37 °C (10,000 gels units/mL). The enzymatic digestion was stopped by addition of EGTA (15 mM). After centrifugation (15,000×g, 15 min, twice) of the digested samples, pelleted genomic DNA fragments were discarded. Concentration of the soluble chromatin fractions was assayed using the fluorescence system Qubit dsDNA HS assay kit (Invitrogen). Prior to immune-precipitation, 5 µg of soluble chromatin was washed in 1.1 mL of Lysis buffer containing 30 µL of magnetic beads (Invitrogen 10001D, 4 h at 4 °C). From the 1.1-mL washed soluble chromatin cleared from magnetic beads, 100 µL was kept for non-IP input control. Specific antibodies were added individually to the remaining 1 mL of the washed soluble chromatin (IP samples) for an overnight incubation at 4 °C on a rotating wheel. Antibodies used in this study are anti-H3K9me3 (Abcam ab8898), anti-H3K27me3 (Millipore 17-622) and anti-H3K4me3 (Active Motif 39159). One additional mock sample was done with an anti-GFP antibody (Abcam ab290) in the wild-type strain only. Chromatin was immuno-precipitated by adding 20 µL of magnetic beads, incubated for 4 h at 4 °C on a rotating wheel. Beads were successively washed in 1 mL for 10 min in of the following buffers: Lysis buffer (twice), Lysis buffer plus 0.5 M NaCl, LiCl Wash Buffer (Tris–HCl pH 8 10 mM, LiCl 250 mM, NP40 0.5%, Sodium deoxycholate 0.5%, EDTA pH 8 1 mM; twice for H3K9me3 immuno-precipitation) and finally Tris–EDTA Buffer (Tris–HCl pH 7.5 10 mM, EDTA pH 8 1 mM). Elution of IP samples was done with TES buffer (Tris–HCl pH 8 50 mM, EDTA pH 8 10 mM with SDS 1%) at 65 °C. TES buffer was also added to non-immuno-precipitated input sample. Both IP and input samples were decrosslinked at 65 °C overnight and treated successively with RNAse A (Thermo Fisher Scientific) (0.2 mg/mL, 2 h at 50 °C) and Proteinase K (Thermo Fisher Scientific) (0.8 mg/mL, 2 h at 50 °C). DNA was then subjected to phenol–chloroform purification and ethanol precipitation. DNA pellets were resuspended in 30 µL of Tris–EDTA buffer.

RT-qPCR experiments

Cultures for RNA extractions were performed for 60 h at 27 °C in liquid culture (KH2PO4 0.25 g/L, K2HPO4 0.3 g/L, MgSO4 0.25 g/L, urea 0.5 g/L, biotin 0.05 mg/L, thiamine 0.05 mg/L, oligo-element 0.1%, yeast extract 5 g/L and dextrin 5.5 g/L). For each replicate, 100 mg of mycelium was harvested, flash-frozen in liquid nitrogen and ground with a Mikro-Dismembrator S (Sartorius, Göttingen, Germany) for one minute at 2600 rpm in a Nalgene Cryogenic vial (ref # 5000-0012, ThermoFischer Scientific, Waltham, USA) with a chromium steel grinding ball (ref # BBI-8546916, Sartorius, Göttingen, Germany). Total RNAs were extracted with the RNeasy Plant Mini Kit (ref # 74904, Qiagen, Hilden, Germany), according to the protocol described for Plants and Fungi with buffer RLT. Contaminating DNA was digested in RNA solutions with RNase-free DNase (ref # 79254, Qiagen, Hilden) and RNAs were purified once more with the RNeasy Plant Mini Kit (ref # 74904, Qiagen, Hilden, Germany), according to the protocol described for RNA cleanup. RNAs were quantified with a DeNovix DS-11 spectrophotometer (Willmington, USA), checked for correct 260/280 and 260/230 ratios, and RNA quality was checked by gel electrophoresis. Reverse transcription was performed with the SuperScript III Reverse Transcriptase (ref # 18080093, Invitrogen, Carlsbad, USA) and oligo-dT. A non-reverse transcribed (NRT) control was systematically performed on a pool of NRT controls to ensure that the Cq was above the Cq obtained from corresponding reverse transcribed RNAs. For genes without introns, NRT control was performed for each replicate sample to ensure that the Cq was above the Cq obtained from corresponding reverse transcribed RNAs. Experiments were run with five biological replicates for each strain, except for some TE cDNA detections (see the section below) (see Additional file 22: Table S4 for Cq). Each biological replicate was run in technical triplicates (see Additional file 22: Table S4 for Cq), except NRT control, which were in technical duplicates (see Additional file 22: Table S4 for Cq). When possible, primers were designed on two consecutive exons (Additional file 23: Table S5). Normalization genes were selected among a set of 8 housekeeping genes using geNorm (see Additional file 21: Table S3 for details of analysis) [105]. Selected normalization genes are AS1 (Pa_1_16650), GPD (Pa_3__5110) and PDF2 (Pa_7_6690), with an average expression stability M = 0.23 indicating a homogeneous set of samples and a V3/4 < 0.15 (see Additional file 22: Table S4 for details of analysis). RT-qPCR normalization according to the ΔΔCt method [106], standard error and 95% confidence interval calculations, and statistical analyses were performed using REST 2009 software (Qiagen, Hilden, Germany). Genes were defined as down-regulated in the mutant strain if the ratio of their transcript level in the mutant strain compared with that in the wild-type strain showed a-fold change > 0 and < 1, with a p-value < 0.05. On the other hand, genes were defined as up-regulated in the mutant strain if the ratio of their transcript level in the mutant strain compared with that in the wild-type strain showed a > onefold change, with a p-value < 0.05. Ratios with a 95% confidence interval, including 1, were not considered significant [107]. RT-qPCR experiments were MIQE compliant [108].

RT-qPCR experiments for detection of TE

Primers for detection of TE were designed to detect as many members as possible for each family (Additional file 23: Table S5). These primer couples detect genomic DNA as well as cDNA. As described above, NRT-qPCR controls were performed with these primer couples on each biological replicate and revealed for some of them that NRT-qPCR signals were too close to the RT-qPCR signals to allow a reliable analysis of the data. Biological replicates that displayed NRT-qPCR signals less than 2.9 cycles above the RT-qPCR signal were discarded and analyses were performed if three biological replicates at least were available (see Additional file 22: Table S4 for details of analysis).

qPCR analysis to test ChIP efficency

Quantitative PCR experiments were performed with dsDNA fluorescent detection method using the FastStart Universal SYBR Green Master from Roche (04913850001).

IP enrichments were assayed with the following protocol: 5 µL of IPed DNA samples and non-IPed input were diluted 10 times in nuclease free water. qPCR experiments were done on 2 µL of the diluted samples using the pairs of primers (0.5 µM) FC77-Actin/FC78-Actin and FC125-580/FC126-580 (see Additional file 23: Table S5). For each IP sample, relative concentrations were determined as an “Input percentage” with the following equation:

$$\% {\text{Input}} = \left( {{\text{E}}^{{\left( { {\text{CqInput}} {-}\left( {\frac{{\ln \left( {10} \right)}}{{{\text{lnE}}}}} \right)} \right) - {\text{CqSample}}}} } \right) \times 100.$$

Sequencing

ChIP-seq libraries were built using the NEB Next UltraII DNA library Prep kit for Illumina (New England Biolabs #E7645S/L) and Agencourt Ampure XP beads (Beckman Coulter, #A63880) according to the manufacturer’s protocols. PCR amplification was made of 12 cycles. ChIP-sequencing was performed by the I2BC High-throughput sequencing facility [NextSeq 500/550 High Output Kit v2 (75 cycles), SR 75 bp]. Reads were trimmed with Cutadapt 1.15 and filtered for control quality by FastQC v0.11.5.

ChIP-seq sequencing data analysis and visualization

Sequencing generated between 8 and 45 million reads depending on samples. The ChIP-seq analysis was performed on sets of data that contained 8 to 10 million reads at most. If needed, the number of reads was down-sized randomly. Reads were mapped on the S mat+ P. anserina genome [51] using Bowtie 2 software (version 2.3.0, see annex for results). For visualization, data were normalized using Deeptools 2.0 [109] bamcoverage with the BPM method. Spearman correlation factors were calculated on normalized data using Deeptools 2.0. For each condition, we made sure that the two biological replicates correlate (Additional file 14: Fig. S14 and Additional file 25: Table S7 and Additional file 26: Table S8) and then merge them for peak calling. Peak calling was performed using MACS2 software (version 2.1.1) using the mock sample as control. Further peak annotations and comparisons were done using Bedtools (version 2.26.0). Genome-wide visualization of peaks was also generated with Circos software [110]. All graphs have been drawn using ggplot2 [111]. Heatmap and metaplot figures have been generated using Deeptools2. Genomic regions have been extracted from the current genome annotation. Promoters have been arbitrarily defined as 1 kb 5′ of the start codon. Scores of the heatmap represent the mapped reads after normalization (see above), hence not only the MACS2-predicted peaks. The Euler diagram has been generated by the eulerr R package (v 6.0.0, https://cran.r-project.org/package=eulerr). The domainogram analysis [50] has been performed using scripts described in [112] and downloaded from the BBCF Github (EPFL, Lausanne, https://github.com/bbcf/). The domainograms were calculated on binned data sets (100 bp binning), using windows sizes ranging from 1 to 500 bins. Individual BRICKs (Blocks of Regulators In Chromosomal Kontext), which represent regions with significant clustering of ChIP-seq signal, were limited to 50 bins each.

To assess whether any of the biological functions observed in H3K4me3 + H3K27me3 gene lists (Additional file 20: Table S2) were present at a frequency greater than that expected by chance, p-values were calculated using hypergeometric distribution as described in [113].

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the NCBI Sequence Read Archive (SRA) (BioProject: PRJNA574032, https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA574032).

References

  1. 1.

    Strahl BD, Allis CD. The language of covalent histone modifications. Nature. 2000;403(6765):41–5.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  2. 2.

    Kouzarides T. Chromatin modifications and their function. Cell. 2007;128(4):693–705.

    CAS  Article  Google Scholar 

  3. 3.

    Bhaumik SR, Smith E, Shilatifard A. Covalent modifications of histones during development and disease pathogenesis. Nat Struct Mol Biol. 2007;14(11):1008–16.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  4. 4.

    Saksouk N, Simboeck E, Déjardin J. Constitutive heterochromatin formation and transcription in mammals. Epigenet Chromatin. 2015;8:3.

    CAS  Article  Google Scholar 

  5. 5.

    Trojer P, Reinberg D. Facultative heterochromatin: is there a distinctive molecular signature? Mol Cell. 2007;28(1):1–13.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  6. 6.

    Martens JHA, O’Sullivan RJ, Braunschweig U, Opravil S, Radolf M, Steinlein P, et al. The profile of repeat-associated histone lysine methylation states in the mouse epigenome. EMBO J. 2005;24(4):800–12.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. 7.

    Nakayama J, Rice JC, Strahl BD, Allis CD, Grewal SI. Role of histone H3 lysine 9 methylation in epigenetic control of heterochromatin assembly. Science. 2001;292(5514):110–3.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  8. 8.

    Rea S, Eisenhaber F, O’Carroll D, Strahl BD, Sun ZW, Schmid M, et al. Regulation of chromatin structure by site-specific histone H3 methyltransferases. Nature. 2000;406(6796):593–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  9. 9.

    Rice JC, Briggs SD, Ueberheide B, Barber CM, Shabanowitz J, Hunt DF, et al. Histone methyltransferases direct different degrees of methylation to define distinct chromatin domains. Mol Cell. 2003;12(6):1591–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  10. 10.

    James TC, Elgin SC. Identification of a nonhistone chromosomal protein associated with heterochromatin in Drosophila melanogaster and its gene. Mol Cell Biol. 1986;6(11):3862–72.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. 11.

    Eissenberg JC, James TC, Foster-Hartnett DM, Hartnett T, Ngan V, Elgin SC. Mutation in a heterochromatin-specific chromosomal protein is associated with suppression of position-effect variegation in Drosophila melanogaster. Proc Natl Acad Sci USA. 1990;87(24):9923–7.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  12. 12.

    Singh PB, Miller JR, Pearce J, Kothary R, Burton RD, Paro R, et al. A sequence motif found in a Drosophila heterochromatin protein is conserved in animals and plants. Nucleic Acids Res. 1991;19(4):789–94.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  13. 13.

    Canzio D, Chang EY, Shankar S, Kuchenbecker KM, Simon MD, Madhani HD, et al. Chromodomain-mediated oligomerization of HP1 suggests a nucleosome-bridging mechanism for heterochromatin assembly. Mol Cell. 2011;41(1):67–81.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Larson AG, Elnatan D, Keenen MM, Trnka MJ, Johnston JB, Burlingame AL, et al. Liquid droplet formation by HP1α suggests a role for phase separation in heterochromatin. Nature. 2017;547(7662):236–40.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  15. 15.

    Sanulli S, Trnka MJ, Dharmarajan V, Tibble RW, Pascal BD, Burlingame AL, et al. HP1 reshapes nucleosome core to promote phase separation of heterochromatin. Nature. 2019;575(7782):390–4.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Strom AR, Emelyanov AV, Mir M, Fyodorov DV, Darzacq X, Karpen GH. Phase separation drives heterochromatin domain formation. Nature. 2017;547(7662):241–5.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. 17.

    Liu S, Brind’Amour J, Karimi MM, Shirane K, Bogutz A, Lefebvre L, et al. Setdb1 is required for germline development and silencing of H3K9me3-marked endogenous retroviruses in primordial germ cells. Genes Dev. 2014;28(18):2041–55.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. 18.

    Becker JS, Nicetto D, Zaret KS. H3K9me3-dependent heterochromatin: barrier to cell fate changes. Trends Genet TIG. 2016;32(1):29–41.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  19. 19.

    Lewis EB. A gene complex controlling segmentation in Drosophila. Nature. 1978;276(5688):565–70.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  20. 20.

    Ringrose L, Paro R. Epigenetic regulation of cellular memory by the Polycomb and Trithorax group proteins. Annu Rev Genet. 2004;38:413–43.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  21. 21.

    Schuettengruber B, Cavalli G. Recruitment of polycomb group complexes and their role in the dynamic regulation of cell fate choice. Dev Camb Engl. 2009;136(21):3531–42.

    CAS  Google Scholar 

  22. 22.

    Wang J, Mager J, Chen Y, Schneider E, Cross JC, Nagy A, et al. Imprinted X inactivation maintained by a mouse Polycomb group gene. Nat Genet. 2001;28(4):371–5.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  23. 23.

    Majewski IJ, Blewitt ME, de Graaf CA, McManus EJ, Bahlo M, Hilton AA, et al. Polycomb repressive complex 2 (PRC2) restricts hematopoietic stem cell activity. PLoS Biol. 2008;6(4):e93.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  24. 24.

    Goodrich J, Puangsomlee P, Martin M, Long D, Meyerowitz EM, Coupland G. A polycomb-group gene regulates homeotic gene expression in Arabidopsis. Nature. 1997;386(6620):44–51.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  25. 25.

    Birve A, Sengupta AK, Beuchle D, Larsson J, Kennison JA, Rasmuson-Lestander A, et al. Su(z)12, a novel Drosophila Polycomb group gene that is conserved in vertebrates and plants. Dev Camb Engl. 2001;128(17):3371–9.

    CAS  Google Scholar 

  26. 26.

    O’Carroll D, Erhardt S, Pagani M, Barton SC, Surani MA, Jenuwein T. The polycomb-group gene Ezh2 is required for early mouse development. Mol Cell Biol. 2001;21(13):4330–6.

    PubMed  PubMed Central  Article  Google Scholar 

  27. 27.

    Pasini D, Bracken AP, Jensen MR, Lazzerini Denchi E, Helin K. Suz12 is essential for mouse development and for EZH2 histone methyltransferase activity. EMBO J. 2004;23(20):4061–71.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. 28.

    Gibson WT, Hood RL, Zhan SH, Bulman DE, Fejes AP, Moore R, et al. Mutations in EZH2 cause Weaver syndrome. Am J Hum Genet. 2012;90(1):110–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  29. 29.

    Li J, Hart RP, Mallimo EM, Swerdel MR, Kusnecov AW, Herrup K. EZH2-mediated H3K27 trimethylation mediates neurodegeneration in ataxia-telangiectasia. Nat Neurosci. 2013;16(12):1745–53.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. 30.

    Wiles ET, Selker EU. H3K27 methylation: a promiscuous repressive chromatin mark. Curr Opin Genet Dev. 2017;43:31–7.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  31. 31.

    Shaver S, Casas-Mollano JA, Cerny RL, Cerutti H. Origin of the polycomb repressive complex 2 and gene silencing by an E(z) homolog in the unicellular alga Chlamydomonas. Epigenetics. 2010;5(4):301–12.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  32. 32.

    Smith KM, Kothe GO, Matsen CB, Khlafallah TK, Adhvaryu KK, Hemphill M, et al. The fungus Neurospora crassa displays telomeric silencing mediated by multiple sirtuins and by methylation of histone H3 lysine 9. Epigenet Chromatin. 2008;1(1):5.

    Article  CAS  Google Scholar 

  33. 33.

    Jamieson K, Rountree MR, Lewis ZA, Stajich JE, Selker EU. Regional control of histone H3 lysine 27 methylation in Neurospora. Proc Natl Acad Sci USA. 2013;110(15):6027–32.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  34. 34.

    Basenko EY, Sasaki T, Ji L, Prybol CJ, Burckhardt RM, Schmitz RJ, et al. Genome-wide redistribution of H3K27me3 is linked to genotoxic stress and defective growth. Proc Natl Acad Sci USA. 2015;112(46):E6339–48.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  35. 35.

    Jamieson K, Wiles ET, McNaught KJ, Sidoli S, Leggett N, Shao Y, et al. Loss of HP1 causes depletion of H3K27me3 from facultative heterochromatin and gain of H3K27me2 at constitutive heterochromatin. Genome Res. 2016;26(1):97–107.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  36. 36.

    Klocko AD, Ormsby T, Galazka JM, Leggett NA, Uesaka M, Honda S, et al. Normal chromosome conformation depends on subtelomeric facultative heterochromatin in Neurospora crassa. Proc Natl Acad Sci USA. 2016;113(52):15048–53.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  37. 37.

    Bicocca VT, Ormsby T, Adhvaryu KK, Honda S, Selker EU. ASH1-catalyzed H3K36 methylation drives gene repression and marks H3K27me2/3-competent chromatin. Elife. 2018;7:e41497.

    PubMed  PubMed Central  Article  Google Scholar 

  38. 38.

    Jamieson K, McNaught KJ, Ormsby T, Leggett NA, Honda S, Selker EU. Telomere repeats induce domains of H3K27 methylation in Neurospora. Elife. 2018;7:e31216.

    PubMed  PubMed Central  Article  Google Scholar 

  39. 39.

    McNaught KJ, Wiles ET, Selker EU. Identification of a PRC2 accessory subunit required for subtelomeric H3K27 methylation in Neurospora. Mol Cell Biol. 2020. https://doi.org/10.1128/MCB.00003-20.

    Article  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Wiles ET, McNaught KJ, Kaur G, Selker JML, Ormsby T, Aravind L, et al. Evolutionarily ancient BAH-PHD protein mediates Polycomb silencing. Proc Natl Acad Sci USA. 2020;117(21):11614–23.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  41. 41.

    Dumesic PA, Homer CM, Moresco JJ, Pack LR, Shanle EK, Coyle SM, et al. Product binding enforces the genomic specificity of a yeast polycomb repressive complex. Cell. 2015;160(1–2):204–18.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  42. 42.

    Chujo T, Scott B. Histone H3K9 and H3K27 methylation regulates fungal alkaloid biosynthesis in a fungal endophyte-plant symbiosis. Mol Microbiol. 2014;92(2):413–34.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  43. 43.

    Connolly LR, Smith KM, Freitag M. The Fusarium graminearum histone H3 K27 methyltransferase KMT6 regulates development and expression of secondary metabolite gene clusters. PLoS Genet. 2013;9(10):e1003916.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  44. 44.

    Studt L, Rösler SM, Burkhardt I, Arndt B, Freitag M, Humpf H-U, et al. Knock-down of the methyltransferase Kmt6 relieves H3K27me3 and results in induction of cryptic and otherwise silent secondary metabolite gene clusters in Fusarium fujikuroi. Environ Microbiol. 2016;18(11):4037–54.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  45. 45.

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

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  46. 46.

    Pham KTM, Inoue Y, Vu BV, Nguyen HH, Nakayashiki T, Ikeda K-I, et al. MoSET1 (histone H3K4 methyltransferase in Magnaporthe oryzae) regulates global gene expression during infection-related morphogenesis. PLoS Genet. 2015;11(7):e1005385.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  47. 47.

    Tamaru H, Selker EU. A histone H3 methyltransferase controls DNA methylation in Neurospora crassa. Nature. 2001;414(6861):277–83.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  48. 48.

    Freitag M, Hickey PC, Khlafallah TK, Read ND, Selker EU. HP1 is essential for DNA methylation in Neurospora. Mol Cell. 2004;13(3):427–34.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  49. 49.

    Silar P, Dauget JM, Gautier V, Grognet P, Chablat M, Hermann-Le Denmat S, et al. A gene graveyard in the genome of the fungus Podospora comata. Mol Genet Genom. 2019;294(1):177–90.

    CAS  Article  Google Scholar 

  50. 50.

    de Wit E, Braunschweig U, Greil F, Bussemaker HJ, van Steensel B. Global chromatin domain organization of the Drosophila genome. PLoS Genet. 2008;4(3):e1000045.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  51. 51.

    Espagne E, Lespinet O, Malagnac F, Da Silva C, Jaillon O, Porcel BM, et al. The genome sequence of the model ascomycete fungus Podospora anserina. Genome Biol. 2008;9(5):R77.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  52. 52.

    Selker EU, Jensen BC, Richardson GA. A portable signal causing faithful DNA methylation de novo in Neurospora crassa. Science. 1987;238(4823):48–53.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  53. 53.

    Gladyshev E. Repeat-induced point mutation and other genome defense mechanisms in fungi. Microbiol Spectr. 2017;5(4):687–99.

    Google Scholar 

  54. 54.

    Gacek A, Strauss J. The chromatin code of fungal secondary metabolite gene clusters. Appl Microbiol Biotechnol. 2012;95(6):1389–404.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  55. 55.

    Studt L, Janevska S, Arndt B, Boedi S, Sulyok M, Humpf H-U, et al. Lack of the COMPASS component Ccl1 reduces H3K4 trimethylation levels and affects transcription of secondary metabolite genes in two plant-pathogenic Fusarium species. Front Microbiol. 2016;7:2144.

    PubMed  PubMed Central  Google Scholar 

  56. 56.

    Grognet P, Lalucque H, Malagnac F, Silar P. Genes that bias Mendelian segregation. PLoS Genet. 2014;10(5):e1004387.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  57. 57.

    Nguyen T-S, Lalucque H, Malagnac F, Silar P. Prions and prion-like phenomena in epigenetic inheritance. In: Tollefsbol TO, editor. Handbook of epigenetics. 2nd ed. Academic Press; 2017. p. 61–72.

    Chapter  Google Scholar 

  58. 58.

    Coppin E, Arnaise S, Contamine V, Picard M. Deletion of the mating-type sequences in Podospora anserina abolishes mating without affecting vegetative functions and sexual differentiation. Mol Gen Genet. 1993;241(3–4):409–14.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  59. 59.

    Jamet-Vierny C, Debuchy R, Prigent M, Silar P. IDC1, a pezizomycotina-specific gene that belongs to the PaMpk1 MAP kinase transduction cascade of the filamentous fungus Podospora anserina. Fungal Genet Biol. 2007;44(12):1219–30.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  60. 60.

    Lambou K, Malagnac F, Barbisan C, Tharreau D, Lebrun M-H, Silar P. The crucial role of the Pls1 tetraspanin during ascospore germination in Podospora anserina provides an example of the convergent evolution of morphogenetic processes in fungal plant pathogens and saprobes. Eukaryot Cell. 2008;7(10):1809–18.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  61. 61.

    Malagnac F, Lalucque H, Lepère G, Silar P. Two NADPH oxidase isoforms are required for sexual reproduction and ascospore germination in the filamentous fungus Podospora anserina. Fungal Genet Biol. 2004;41(11):982–97.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  62. 62.

    Honda S, Selker EU. Direct interaction between DNA methyltransferase DIM-2 and HP1 is required for DNA methylation in Neurospora crassa. Mol Cell Biol. 2008;28(19):6044–55.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  63. 63.

    Kumar A, Kono H. Heterochromatin protein 1 (HP1): interactions with itself and chromatin components. Biophys Rev. 2020;12(2):387–400.

    PubMed  PubMed Central  Article  Google Scholar 

  64. 64.

    Freitag M. Histone methylation by SET domain proteins in fungi. Annu Rev Microbiol. 2017;08(71):413–39.

    Article  CAS  Google Scholar 

  65. 65.

    Simon JA, Kingston RE. Occupying chromatin: Polycomb mechanisms for getting to genomic targets, stopping transcriptional traffic, and staying put. Mol Cell. 2013;49(5):808–24.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  66. 66.

    Voigt P, LeRoy G, Drury WJ, Zee BM, Son J, Beck DB, et al. Asymmetrically modified nucleosomes. Cell. 2012;151(1):181–93.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  67. 67.

    Frapporti A, Miró Pina C, Arnaiz O, Holoch D, Kawaguchi T, Humbert A, et al. The Polycomb protein Ezl1 mediates H3K9 and H3K27 methylation to repress transposable elements in Paramecium. Nat Commun. 2019;20:10.

    Google Scholar 

  68. 68.

    Ho JWK, Jung YL, Liu T, Alver BH, Lee S, Ikegami K, et al. Comparative analysis of metazoan chromatin organization. Nature. 2014;512(7515):449–52.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  69. 69.

    Montgomery SA, Tanizawa Y, Galik B, Wang N, Ito T, Mochizuki T, et al. Chromatin organization in early land plants reveals an ancestral association between H3K27me3, transposons, and constitutive heterochromatin. Curr Biol. 2020;30(4):573.e7-588.e7.

    Article  CAS  Google Scholar 

  70. 70.

    Alder O, Lavial F, Helness A, Brookes E, Pinho S, Chandrashekran A, et al. Ring1B and Suv39h1 delineate distinct chromatin states at bivalent genes during early mouse lineage commitment. Dev Camb Engl. 2010;137(15):2483–92.

    CAS  Google Scholar 

  71. 71.

    Widschwendter M, Fiegl H, Egle D, Mueller-Holzner E, Spizzo G, Marth C, et al. Epigenetic stem cell signature in cancer. Nat Genet. 2007;39(2):157–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  72. 72.

    Ohm JE, McGarvey KM, Yu X, Cheng L, Schuebel KE, Cope L, et al. A stem cell-like chromatin pattern may predispose tumor suppressor genes to DNA hypermethylation and heritable silencing. Nat Genet. 2007;39(2):237–42.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  73. 73.

    Schlesinger Y, Straussman R, Keshet I, Farkash S, Hecht M, Zimmerman J, et al. Polycomb-mediated methylation on Lys27 of histone H3 pre-marks genes for de novo methylation in cancer. Nat Genet. 2007;39(2):232–6.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  74. 74.

    Deleris A, Stroud H, Bernatavichute Y, Johnson E, Klein G, Schubert D, et al. Loss of the DNA methyltransferase MET1 Induces H3K9 hypermethylation at PcG target genes and redistribution of H3K27 trimethylation to transposons in Arabidopsis thaliana. PLoS Genet. 2012;8(11):e1003062.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  75. 75.

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

    Article  CAS  Google Scholar 

  76. 76.

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

    Article  Google Scholar 

  77. 77.

    Soyer JL, Clairet C, Gay EJ, Lapalu N, Rouxel T, Stukenbrock EH, et al. Genome-wide mapping of histone modifications in two species of Leptosphaeria maculans showing contrasting genomic organization and host specialization. bioRxiv. 2020. https://doi.org/10.1101/2020.05.08.084566.

    Article  Google Scholar 

  78. 78.

    Thalheim T, Herberg M, Loeffler M, Galle J. The regulatory capacity of bivalent genes—a theoretical approach. Int J Mol Sci. 2017;18(5):1069.

    PubMed Central  Article  CAS  Google Scholar 

  79. 79.

    Roper M, Ellison C, Taylor JW, Glass NL. Nuclear and genome dynamics in multinucleate ascomycete fungi. Curr Biol. 2011;21(18):R786–93.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  80. 80.

    Money NP. Mushroom stem cells. BioEssays News Rev Mol Cell Dev Biol. 2002;24(10):949–52.

    CAS  Article  Google Scholar 

  81. 81.

    Paoletti M, Saupe SJ. Fungal incompatibility: evolutionary origin in pathogen defense? BioEssays News Rev Mol Cell Dev Biol. 2009;31(11):1201–10.

    CAS  Article  Google Scholar 

  82. 82.

    Gu Q, Wang Z, Sun X, Ji T, Huang H, Yang Y, Zhang H, Tahir HAS, Wu L, Wu H, Gao X. FvSet2 regulates fungal growth, pathogenicity, and secondary metabolism in Fusarium verticillioides. Fungal Genet Biol. 2017;107:24–30.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  83. 83.

    Reyes-Dominguez Y, Boedi S, Sulyok M, Wiesenberger G, Stoppacher N, Krska R, Strauss J. Heterochromatin influences the secondary metabolite profile in the plant pathogen Fusarium graminearum. Fungal Genet Biol FG B. 2012;49:39–47.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  84. 84.

    Chujo T, Lukito Y, Eaton CJ, Dupont P-Y, Johnson LJ, Winter D, Cox MP, Scott B. Complex epigenetic regulation of alkaloid biosynthesis and host interaction by heterochromatin protein I in a fungal endophyte-plant symbiosis. Fungal Genet Biol. 2019;125:71–83.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  85. 85.

    Reyes-Dominguez Y, Bok JW, Berger H, Shwab EK, Basheer A, Gallmetzer A, Scazzocchio C, Keller N, Strauss J. Heterochromatic marks are associated with the repression of secondary metabolism clusters in Aspergillus nidulans. Mol Microbiol. 2010;76:1376–86.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  86. 86.

    Palmer JM, Perrin RM, Dagenais TRT, Keller NP. H3K9 Methylation Regulates Growth and Development in Aspergillus fumigatus. Eukaryot Cell. 2008;7:2052–60.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  87. 87.

    Soyer JL, Ghalid ME, Glaser N, Ollivier B, Linglin J, Grandaubert J, Balesdent M-H, Connolly LR, Freitag M, Rouxel T, Fudal I. Epigenetic Control of Effector Gene Expression in the Plant Pathogenic Fungus Leptosphaeria maculans. PLOS Genet. 2014;10:e1004227

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  88. 88.

    Zhang X, Qu Y, Qin Y. Expression and chromatin structures of cellulolytic enzyme gene regulated by heterochromatin protein 1. Biotechnol Biofuels. 2016;9:206.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  89. 89.

    Angel A, Song J, Dean C, Howard M. A Polycomb-based switch underlying quantitative epigenetic memory. Nature. 2011;476(7358):105–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  90. 90.

    Lewis ZA, Honda S, Khlafallah TK, Jeffress JK, Freitag M, Mohn F, et al. Relics of repeat-induced point mutation direct heterochromatin formation in Neurospora crassa. Genome Res. 2009;19(3):427–37.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  91. 91.

    Veluchamy A, Rastogi A, Lin X, Lombard B, Murik O, Thomas Y, et al. An integrative analysis of post-translational histone modifications in the marine diatom Phaeodactylum tricornutum. Genome Biol. 2015;20(16):102.

    Article  CAS  Google Scholar 

  92. 92.

    Kamei M, Ameri AJ, Ferraro AR, Bar-Peled Y, Zhao F, Ethridge CL, et al. IMITATION SWITCH is required for normal chromatin structure and gene repression in PRC2 target domains. Proc Natl Acad Sci. 2021. https://doi.org/10.1073/pnas.2010003118.

    Article  PubMed  PubMed Central  Google Scholar 

  93. 93.

    Grognet P, Bidard F, Kuchly C, Tong LCH, Coppin E, Benkhali JA, et al. Maintaining two mating types: structure of the mating type locus and its role in heterokaryosis in Podospora anserina. Genetics. 2014;197(1):421–32.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  94. 94.

    Rizet G, Engelmann C. Contribution à l’étude génétique d’un Ascomycète tétrasporé: Podospora anserina. Rev Cytol Biol Veg. 1949;11:201–304.

    Google Scholar 

  95. 95.

    Silar P. Podospora anserina: from laboratory to biotechnology. In: Horwitz BA, Mukherjee PK, Mukherjee M, Kubicek CP, editors. Genomics of soil- and plant-associated fungi. Berlin: Springer; 2013. p. 283–309. https://doi.org/10.1007/978-3-642-39339-6_12.

    Chapter  Google Scholar 

  96. 96.

    Silar P. Podospora anserina. 2020. https://hal.archives-ouvertes.fr/hal-02475488.

  97. 97.

    Ausubel F, Brent R, Kingston R, Moore D, Seidman J, Smith J, et al. Current protocols in molecular biology. New York: Wiley; 1987.

    Google Scholar 

  98. 98.

    Lecellier G, Silar P. Rapid methods for nucleic acids extraction from Petri dish-grown mycelia. Curr Genet. 1994;25(2):122–3.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  99. 99.

    Brygoo Y, Debuchy R. Transformation by integration in Podospora anserina. I. Methodology and phenomenology. Mol Gen Genet. 1985;200:128–31.

    CAS  Article  Google Scholar 

  100. 100.

    Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  101. 101.

    Grognet P, Timpano H, Carlier F, Aït-Benkhali J, Berteaux-Lecellier V, Debuchy R, et al. A RID-like putative cytosine methyltransferase homologue controls sexual development in the fungus Podospora anserina. PLoS Genet. 2019;15(8):e1008086.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  102. 102.

    Silar P, Picard M. Increased longevity of EF-1 alpha high-fidelity mutants in Podospora anserina. J Mol Biol. 1994;235(1):231–6.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  103. 103.

    Grigoriev IV, Nikitin R, Haridas S, Kuo A, Ohm R, Otillar R, et al. MycoCosm portal: gearing up for 1000 fungal genomes. Nucleic Acids Res. 2014;42(Database issue):D699–704.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  104. 104.

    Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59(3):307–21.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  105. 105.

    Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3(7):1–12.

    Article  Google Scholar 

  106. 106.

    Pfaffl MW, Horgan GW, Dempfle L. Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002;30(9):e36.

    PubMed  PubMed Central  Article  Google Scholar 

  107. 107.

    du Prel J-B, Hommel G, Röhrig B, Blettner M. Confidence interval or p-value?: part 4 of a series on evaluation of scientific publications. Dtsch Arzteblatt Int. 2009;106(19):335–9.

    Google Scholar 

  108. 108.

    Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55(4):611–22.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  109. 109.

    Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 2016;44(1):W160–5.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  110. 110.

    Krzywinski MI, Schein JE, Birol I, Connors J, Gascoyne R, Horsman D, et al. Circos: an information aesthetic for comparative genomics. Genome Res. 2009;19(9):1639–45.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  111. 111.

    Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer; 2009. (Use R!).

    Book  Google Scholar 

  112. 112.

    David FPA, Delafontaine J, Carat S, Ross FJ, Lefebvre G, Jarosz Y, et al. HTSstation: a web application and open-access libraries for high-throughput sequencing data analysis. PLoS ONE. 2014;9(1):e85879.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  113. 113.

    Boyle EI, Weng S, Gollub J, Jin H, Botstein D, Cherry JM, et al. GO::TermFinder–open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes. Bioinforma Oxf Engl. 2004;20(18):3710–5.

    CAS  Article  Google Scholar 

  114. 114.

    Qian C, Zhou M-M. SET domain protein lysine methyltransferases: structure, specificity and catalysis. Cell Mol Life Sci. 2006;63(23):2755–63.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  115. 115.

    Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  116. 116.

    Lalucque H, Malagnac F, Brun S, Kicka S, Silar P. A non-Mendelian MAPK-generated hereditary unit controlled by a second MAPK pathway in Podospora anserina. Genetics. 2012;191(2):419–33.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We acknowledge the technical assistance of Sylvie François. We acknowledge the high-throughput sequencing facility of I2BC for its sequencing and bioinformatics expertise. We are thankful to Sébastien Bloyer for providing expertise to start the project as well as for fruitful scientific discussions, to Benoît Moindrot for his help and for providing scripts for the domainogram analysis, and to Gaëlle Lelandais for her help with statistical analysis. Pierre Grognet thanks the "École de bioinformatique AVIESAN-IFB" for its excellent training. We thank Linda Sperling for editing and critical reading of the manuscript.

Funding

FC was recipient of a LIDEX BIG Paris-Saclay Ph.D. fellowship. ML is recipient of a Chinese Scholarship Council PhD fellowship. PG, FC and FM were supported by grants from UMR9198.

Author information

Affiliations

Authors

Contributions

Conceived and designed the experiments: FM. Performed the experiments: FC, ML, LM, RD, PG, FM. Analyzed the data: FC, RD, DN, PG, FM. Contributed reagents/materials/analysis tools/expertise sharing: PG, CS, DN. Wrote the paper: FC, DN, PG, FM. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to P. Grognet or F. Malagnac.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

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

Supplementary Information

Additional file 1: Figure S1.

Conserved structure and phylogenetic analysis of histone methyltransferases involved in heterochromatin assembly. Domain structure comparison of histone methyltransferase Kmt1 (A) and Kmt6 (B). Sizes in amino acid (aa) are given (right). Pre-SET (red, IPR007728), SET (orange, IPR001214) and Post-SET (green, IPR003616) conserved domains are required for H3K9 methyltransferase activity of Kmt1 homolog proteins. SRA (SET and RING finger-associated, purple, IPR003105) and SANT (blue, IPR001005) are protein–protein interaction domains. CXC (yellow, IPR026489) is a cysteine-rich conserved domain located in the H3K27 methyltransferase catalytic domain of Kmt6 homologs. Filamentous fungi: Podospora anserina (Pa), Neurospora crassa (Nc), Fusarium graminearum (Fg), Epichloë festucae (Ef), Aspergillus nidulans (An); yeasts: Schizosaccharomyces pombe (Sp) and Cryptococcus neoformans (Cn); the worm Caenorhabditis elegans (Ce), the fruit-fly Drosophila melanogaster (Dm), the mouse Mus musculus (Mm) and the model plant Arabidopsis thaliana (At). Accession numbers for proteins used in alignments are listed in Additional file 24: Table S6. C Alignment of the SET domains within SU(VAR)3–9 (upper panel) and EZH2 (lower) homologs for a set of selected species. Deleted parts of proteins are indicated by // (shaded in gray). Conserved GxG motif (shaded in green) corresponds to the substrate AdoMet/SAM binding site (motif I). Motif II displays the conserved YxG motif, motif III the YLF triplet. Folding of the histone methyltransferases in knot-like structure brings together the catalytic residues NH and Y (shaded in yellow), embedded in the conserved signature motifs IV and VI (in bold). The F residue of the FxY motif (shaded in blue) is responsible for methylation specificity of SU(VAR)3–9 [114]. If changed to a Y, H3K9me3 modification is converted to H3K9me2. Mm: Mus musculus, Dm: Drosophila melanogaster, Sp: Schizosaccharomyces pombe, Cn: Cryptococcus neoformans, Fg: Fusarium graminearum, Ef: Epichloë festucae, Nc: Neurospora crassa, Pa: Podospora anserina. D Alignment of the pre/post-SET domains within SU(VAR)3–9 homologs for a set of selected species. Deleted parts of proteins are indicated by // (shaded in gray). Conserved pre-SET motif Zn3Cys9 (shaded in yellow) and post-SET CxCx4C (shaded in blue) [114]. Mm: Mus musculus, Dm: Drosophila melanogaster, Sp: Schizosaccharomyces pombe, Fg: Fusarium graminearum, Ef: Epichloë festucae, Nc: Neurospora crassa, Pa: Podospora anserina. Accession numbers for proteins used in alignments are listed in Additional file 24: Table S6.

Additional file 2: Figure S2.

Expression kinetics PaKmt1 in wild-type strain. Upper panel: RT-PCR on RNA extraction from wild-type vegetative growing mycelium (1 day and 4 days), perithecia (2 days and 4 days after the fertilization), and input genomic DNA. PaKmt1 CDS was predicted to be made of two exons separated by a 62 bp intron (positions 49–110). However, two amplicons of distinct sizes are obtained, corresponding to both spliced (1108 bp) and unspliced (1170 bp) PaKmt1 mRNAs. MW: molecular weight. Lower panel: schematic representation of the PaKmt1 locus (DNA). mRNA1 corresponds to the spliced form of the transcripts while mRNA2 (1108 bp) corresponds to the unspliced form (1170 bp). Translation of the unspliced form would lead to a premature termination and thus to a truncated protein. Primers used for the reverse-transcription polymerase chain reaction (RT-PCR) are drawn as arrows above and below the PaKmt1 CDS.

Additional file 3: Figure S3. A

Heatmap of Spearman’s correlation coefficient comparison: clustering analysis of histone marks in wild-type background (WT). 2e and 3d are the two WT strains used for this study. They are issued from two spores from the same WT cross. Mock = IP performed with GFP antibody in absence of GFP tag in P. anserina’s genome (see “Methods”). Raw data are given in Additional file 25: Table S7. B H3K27me3, H3K4me3 and H3K9me3 proportion on P. anserina chromosomes in the WT, ΔPaKmt1, ΔPaKmt6 and ΔPaHP1 mutant strains. Plot showing the percentage of each chromosome covered with H3K4me3 (green), H3K9me3 (red) and H3K27me3 (blue). The coverage is the sum of all MACS2-predicted peak sizes. C Normalized ChIP-seq data representation for all marks on the seven P. anserina chromosomes for all conditions. ChIP-seq patterns display histone modification coverage and MACS2 detected peaks. D Domainogram representations for all marks on the seven P. anserina chromosomes for all conditions. Domainograms show significance of enrichment of H3K4me3, H3K9me3, H3K27me3 marks in windows of varying size. Color-coding of p-value is indicated (top).

Additional file 4: Figure S4.

Combined epigenetic landscapes in wild-type and heterochromatin mutant strains of P. anserina. Panorama of genome-wide peak localization for each genotype, wild-type, ΔPaKmt1, ΔPaKmt6 and ΔPaHP1 strains. Telomeres sequences were arbitrarily defined as the segment going from the end of each arm of the chromosomes to the first annotated gene (with the exception of the rDNA cluster localized on chromosome 3) and centromeres are indicated. Mat region = Non-recombining region containing the mating-type locus as defined in [93]. A segment overlapping portions of chromosomes 3 and 4 is expanded to show a zoom of the combined epigenetic landscapes.

Additional file 5: Figure S5.

H3K4me3, H3K27me3 and H3K9me3 modifications of P. anserina TEs in the WT, ΔPaKmt1, ΔPaKmt6 and ΔPaHP1 mutant strains. A. Histone marks on P. anserina TE families. Top panel: Plots of normalized ChIP signal: H3K4me3 (green), H3K9me3 (red) and H3K27me3 (dark blue) signals in the wild-type strain for five TE families, i.e., Copia, Gypsy, MITE, Tc1 Mariner, solo LTR and unclassified TEs [51] (Additional file 21: Table S3). Because MITE TEs are shorter than the other categories (< 500 bp in length), the genomic window was narrowed on the graph. Aligned sequences correspond to TE bodies ± 0.5 kbp surrounding region (see Additional file 21: Table S3 for TE’s numbers). Bottom panel: K-means built clusters representing the association versus non-association of the indicated histone modifications within a specific TE family. Histone modification levels in the heatmaps were calculated for non-overlapping 10-bp windows within the specific genomic regions and sorted by average value of each row. B Number of TEs, classified by family, marked with H3K4me3, H3K27me3 or H3K9me3, according to genetic backgrounds (wild-type strain, ΔPaKmt1, ΔPaKmt6 and ΔPaHP1 mutant strains). C Violin plots of expression of TEs classified by family. Gene expression was inferred from the TPM (Transcripts Per Kilobase Million) values calculated in [49]. Gene expression of non-repeated CDS (i.e., “genes”) was added for comparison.

Additional file 6: Figure S6.

Snapshots of a set of TEs representing all the annotated TE families in P. anserina’s genome. ChIP-seq signals were normalized as described in “Methods” and visualized using the Integrated Genomics Viewer (IGV) [115]. H3K4me3 (green), H3K9me3 (red) and H3K27me3 (blue).

Additional file 7: Figure S7.

Antibody specificity analysis. Dot blot results using H3K9me3 and H3K27me3 peptides (left) at different concentration (top). Membranes were inoculated with the corresponding antibody (right). Signal intensity comparison shows a cross-reactivity between anti-H3K9me3 antibody and H3K27me3 evaluated at 3% compared to the immunogen reaction.

Additional file 8: Figure S8.

Molecular characterization of knockout mutants by Southern blot hybridization. Schematic representations of the endogenous and disrupted loci are given (Left). Replacement by homologous recombination of the wild-type PaKmt1 allele by the disrupted ΔPaKmt1 allele results in the substitution of 2.4 kbp and 5.7 kbp EcoRV fragments by a unique 11 kbp PstI fragment as revealed by hybridization of the 5′ and 3′ digoxygenin-labeled probes (dashed rectangles PaKmt1 locus). Replacement by homologous recombination of the wild-type PaKmt6 allele by the disrupted ΔPaKmt6 allele results in the substitution of a unique 6.5 kbp KpnI fragment by two 1.8 and 2.4 kbp KpnI fragments as revealed by hybridization of the 5′ and 3′ digoxygenin-labeled probes (dashed rectangles PaKmt6 locus). A second verification has been made with the HindIII enzymes and the same probes shows the substitution of two 1.9 kbp and 4.6 kbp HindIII fragments by a 3.7 kbp HindIII fragment as revealed by hybridization of the 5′ and 3′ digoxygenin-labeled probes.

Additional file 9: Figure S9.

Localization of histone marks on specific genomic regions in the ΔPaKmt1 and ΔPaKmt6 mutant strains. Top panel: Plots of normalized ChIP-seq signal. Bottom panel: Heatmaps divided in K-means built clusters representing the association versus non-association of the indicated histone modifications with the specific genomic regions. Coding sequences or CDS were aligned by their two ends (indicated by START and STOP) ± 1 kbp of surrounding sequence (N = 10,839; Additional file 20: Table S2); repeats were defined as TE bodies, duplications and the rDNA array ± 0.2 kbp surrounding regions (N = 1680; Additional file 21: Table S3). Histone modification levels in the heatmaps were calculated for non-overlapping 10 bp windows within the specific genomic regions and sorted by average value of each row.

Additional file 10: Figure S10.

Relative expression of selected genes in the ΔPaKmt1 mutants. Caption as in Fig. 5a. The error bars represent the 95% confidence interval. No significant fold change was detected, except for Pa_1_6263 and Pa_6_7370, which both lost the H3K27me3 mark and were up-regulated and down-regulated, respectively. Pa_1_6263: expression ratio = 1.855, p-value = 0.002; Pa_4_1170: expression ratio = 1.535, p-value = 0.094; Pa_1_16300: expression ratio = 1.393, p-value = 0.092; Pa_6_7270: expression ratio = 1.397, p-value = 0.035; Pa_6_7370: expression ratio = 0.696, p-value = 0.006; Pa_5_10: expression ratio = 0.759, p-value = 0.017; Pa_1_1880: expression ratio = 1.422, p-value = 0.034; Pa_7_9210: expression ratio = 1.082, p-value = 0.252; TC1mlr represent quantification of the cDNAs from the members of the Tc1_mariner-like_rainette family: expression ratio = 1.127, p-value = 0.339; TC1mlp represent quantification of the cDNAs from the members of the Tc1_mariner-like_pelobates family: expression ratio = 0.851, p-value = 0.308. Copia_Ty1_nephelobates transcripts could not be quantified as NRT-qPCR controls were too close to RT-qPCR (see Additional file 21: Table S3 for details of analysis).

Additional file 11: Figure S11.

Growth features of the ΔPaKmt1 strain. A Experimental procedure to test ability to resume growth of the ΔPaKmt1 strain. To set up this restart test, mycelium implants issued from germination thalli were inoculated onto fresh M2 medium and incubated at 27 °C for 8 days (step 1). Two independent plugs from each location (purple cross), 8-day stationary phase (plug#3); pink cross, 4-day stationary phase (plug#2); and green cross, growing phase (plug#1) were then transferred to fresh M2 solid medium and incubated for 3 days at 27 °C (step 2). The growing phase of each thallus originating from step 2 (green margins) were transplanted again to fresh M2 solid medium and incubated at 27 °C for 3 days (step 3). Growth restart from stationary phase was impaired for ΔPaKmt1 mutants, which resulted in smaller and thinner colonies than the wild-type ones (white arrows), whereas continuous growth (plug#1) was not altered. Complemented ΔPaKmt1-PaKmt1+ strains behaved as wild-type strains. As control experiments (step 3), we then transferred mycelia from growing margins (marked in green, step 2) of thalli deriving from Plug#1, Plug#2 and Plug#3. In this case, ΔPaKmt1 mutants did not show any delay to resume growth (orange arrows), confirming that this defect was not permanent but rather linked to the disability of the ΔPaKmt1 strains to resume growth properly. B Crippled growth test for ΔPaKmt1 strain. Crippled growth (CG) process can be shown using a ‘band test’. Strains were incubated on M2 medium for 7 days at 27 °C. Two 1-mm-wide slices of agar were then inoculated onto fresh M2 media with or without yeast extract (YE) in the following method for 3 days. Picture shows actively growing apical hyphae (right) and resting stationary phase hyphae (left). The surface side of the slice on plate has been orientated at the top. CG is a degenerative process caused by C element production in the stationary phase. When inoculated on yeast extract medium, it displays slow growth, alteration of pigmentation, inability to differentiate aerial hyphae, and female sterility. The ΔPaMpk1 mutant strain is impaired for CG development while the PDC2208 strain is supposed to display CG on M2 medium without yeast extract [116]. The single ΔPaKmt1 mutants are not impaired for CG.

Additional file 12: Figure S12.

Longevity tests for ΔPaKmt1 and ΔPaKmt6 strains. Graphs show the maximal growth length on M2 medium at 27 °C in race tubes for wild-type (WT), ΔPaKmt1 and ΔPaKmt6 mutant strains for both mating types. Data correspond to the mean of three technical replicates from five biological samples.

Additional file 13: Figure S13.

Role of PaKmt6 during sexual development. A Heterozygote oriented crosses between wild-type strains (WT) and ΔPaKmt6 mutants. Fertilization of wild-type ascogonia with ΔPaKmt6 male gametes results in normal fruiting body development displaying a single neck (black arrow), while fertilization of ΔPaKmt6 ascogonia with wild-type spermatia results in fewer crippled ΔPaKmt6-like fruiting bodies (two necks are indicated by two black arrows). These features indicate that PaKmt6 is a maternal gene. B Mosaic analyses using the Δmat mutant strain. Wild-type or ΔPaKmt6 mat+ and mat− strains were mixed with or without Δmat mycelium and inoculated onto fresh M2 medium. After a week of incubation on M2 medium, ΔPaKmt6 dikaryon displayed crippled and mis-orientated perithecia, resulting in scattered and reduced ascospore production. The tricaryon ΔPaKmt6 mat+ /ΔPaKmt6 mat−mat showed nearly wild-type restoration of mycelium growth, as well as perithecia and ascospore production. These features confirm that PaKmt6 is a gene expressed in the maternal tissues of the perithecium and not in its zygotic tissues (centrum).

Additional file 14: Figure S14.

Structure and evolutionary relationships of HP1 orthologs. A Domain structure comparison of heterochromatin protein 1 homologs. Size in amino acids (aa) is given (right). PaHP1 displays the evolutionary conserved and topologically connected chromodomain (turquoise) and chromoshadowdomain (dark green), along with three disordered elements: the N-terminal extension (NTE), the Hinge region (HR) and the C-terminal extension (CTE). Chromodomain recognizes H3K9me2/3 histone modification, while chromoshadowdomain is a protein–protein interaction domain. B Alignment of the chromodomains of HP1 homologs for a set of selected species. The aromatic cage residues (three aromatic amino acids shaded in yellow, i.e., Y and WxxY) form the binding cavity for H3K9me. Residues involve in the H3 peptide recognition surface are highlighted in blue. Residues of chromo-shadow domain involved in dimerization are shown in green, those involved in folding are shown in lavender [63].

Additional file 15: Figure S15. A

Tree showing evolution of some HP1 orthologs from fungi, plant and metazoans. The P. anserina heterochromatin protein-1 homolog is marked with a rectangle. Bootstraps are given. Filamentous fungi: Podospora anserina (Pa), Neurospora crassa (Nc), Fusarium graminearum (Fg), Trichoderma reesei (Tr), Epichloë festucae (Ef) Botrytis cinerea (Bc), Magnaporthe oryzae (Mo), Zymoseptoria tritici (Zt), Leptosphaeria maculans (Lm), Aspergillus nidulans (An), Aspergillus fumigatus (Af), Penicillium oxalicum (Po), Ascobolus immersus (Ai), Puccinia graminis (Pg), Pneumocystis jirovecii (Pj), Conidiobolus coronatus (Cc), Mucor circinelloides (Mc); yeasts: Schizosaccharomyces pombe (Sp) and Cryptococcus neoformans (Cn); the worm Caenorhabditis elegans (Ce), the fruit-fly Drosophila melanogaster (Dm), the mouse Mus musculus (Mm) and the model plant Arabidopsis thaliana (At). Accession numbers for proteins used in alignments are listed in Additional file 24: Table S6. B Molecular characterization of knockout ΔPaHP1 mutants by Southern blot hybridization. Replacement by homologous recombination of the wild-type PaHP1 allele by the disrupted ΔPaHP1 allele results in the substitution of a unique 3.4 kb EcoRI/EcoRV fragment by two 1.4 and 2.5 kb EcoRI/EcoRV fragments as revealed by hybridization of the 5′UTR and 3′UTR digoxygenin-labeled probes (dashed rectangles PaHP1 locus). C Vegetative growth of ΔPaHP1 single mutants and ΔPaKmt1ΔPaHP1 double mutants compared to wild-type strains. ΔPaKmt1 and ΔPaHP1 single mutants are impaired in aerial mycelium production. This defect was more pronounced for the single ΔPaHP1 mutants. ΔPaKmt1ΔPaHP1 double mutants show ΔPaKmt1 morphologic phenotype. The strains have been grown on M2 medium 3 days.

Additional file 16: Figure S16.

Localization of histone marks on specific genomic regions in the ΔPaPH1 mutant strains. Top panel: Plots of normalized ChIP-seq signal. Bottom panel: heatmaps divided in K-means built clusters representing the association versus non-association of the indicated histone modifications with the specific genomic regions. Coding sequences or CDS were aligned by their two ends (indicated by START and STOP) ± 1 kbp of surrounding sequence (N = 10,839; Additional file 20: Table S2); repeats were defined as TE bodies, duplications and the rDNA array ± 0.2 kbp surrounding regions (N = 1680; Additional file 20: Table S3). Histone modification levels in the heatmaps were calculated for non-overlapping 10 bp windows within the specific genomic regions and sorted by average value of each row.

Additional file 17: Figure S17.

A Relative expression of selected genes in the ΔPaHP1 strain. Caption as in Fig. 5a. The error bars represent the 95% confidence interval. No significant fold change was assayed, except for Pa_1_6263 and Pa_5_10, which both lost the H3K27me3 mark and were up-regulated and down-regulated, respectively. Pa_1_6263: expression ratio = 2.463, p-value = 0.004; Pa_4_1170: expression ratio = 1.911, p-value = 0.032; Pa_1_16300: expression ratio = 0.931, p-value = 0.693; Pa_6_7270: expression ratio = 1.527, p-value = 0.004; Pa_6_7370: expression ratio = 0.766, p-value = 0.009; Pa_5_10: expression ratio = 0.490, p-value = 0.006; Pa_1_1880: expression ratio = 1.313, p-value = 0.013; Pa_7_9210: expression ratio = 0.951, p-value = 0.547; TC1mlr represents quantification of the cDNAs from the members of the Tc1_mariner-like_rainette family: expression ratio = 1.264, p-value = 0.008. Tc1_mariner-like_pelobates and copia_Ty1_nephelobates transcripts could not be quantified as NRT-qPCR controls were too close to RT-qPCR (see Additional file 21: Table S3 for details of analysis). B Vegetative growth kinetics of ΔPaHP1 mutants and ΔPaKmt1 mutants compared to wild-type strains. See “Methods” section for details. C Experimental procedure to test growth resuming capabilities of ΔPaHP1 and double ΔPaHP1ΔPaKmt1 mutants. For description of experimental settings see Additional file 8: Fig. S8A. Growth restart from stationary phase (plug#2 and plug#3) was impaired for ΔPaHP1 and double ΔPaHP1ΔPaKmt1 mutants, which resulted in smaller and thinner colonies than the wild-type ones (white arrows), whereas continuous growth (plug#1) was not altered. Complemented ΔPaHP1-PaHP1+ strains behaved as wild-type strains. As control experiments (step 3), we transferred mycelia from growing margins (marked in green, step 2) of thalli derived from Plug#1, Plug#2 and Plug#3. In this case, neither ΔPaKmt1 mutants nor ΔPaHP1ΔPaKmt1 mutants showed any delay to resume growth (orange arrows). D Crippled growth test of ΔPaHP1 single mutants and ΔPaHP1ΔPaKmt1 double mutants. For description of experimental settings see Additional file 8: Fig. S8B. The single ΔPaHP1 and the double ΔPaKmt1ΔPaHP1 mutants are not impaired for CG.

Additional file 18: Figure S18.

Heat map of Spearman’s correlation coefficient comparison. Clustering analysis of histone marks in wild-type background (WT) and in mutant backgrounds ΔPaKmt1, ΔPaKmt6 and ΔPaHP1. Mock = IP performed with GFP antibody in the absence of any GFP tag in the P. anserina genome (see “Methods”). Raw data are given in Additional file 26: Table S8.

Additional file 19: Table S1.

Main features of significantly enriched peaks. Peak numbers, sizes, and genome coverage of H3K4me3, H3K9me3 and H3K27me3 modifications, in wild-type background (WT) and in mutant backgrounds ΔPaKmt1, ΔPaKmt6 and ΔPaHP1.

Additional file 20: Table S2.

Distribution of H3K4me3, H3K9me3 and H3K27me3 on the complete set of P. anserina genes in wild-type background (WT) and in ΔPaKmt1, ΔPaKmt6 and ΔPaHP1 mutants. List of genes where H3K4me3 and H3K27me3 can be found overlapping in the wild-type background. Distribution of H3K4me3, H3K9me3 and H3K27me3 on a subset of P. anserina secondary metabolite gene clusters in wild-type background (WT) and in PaKmt1, ΔPaKmt6 and ΔPaHP1 mutants. ChIP_peaks, 0 means no peaks detected by MACS2, 1 means at least one peak detected by MACS2. ChIP_coverage: mapped reads after normalization (see “ChIP-seq sequencing data analysis and visualization” paragraph in “Methods” for details). All_BRICKS: H3K4me3, H3K9me3 or H3K27me3 domainograms were simplified as described in [50] to a set of discrete domains of enrichment called Blocks of Regulators In Chromosomal Kontext (BRICKs). BRICKS can be used to define domains with significant clustering of ChIP-seq signal and their boundaries, while domainograms provide an overview.

Additional file 21: Table S3.

Distribution of H3K4me3, H3K9me3 and H3K27me3 on the complete set of P. anserina transposable elements in wild-type background (WT) and in ΔPaKmt1, ΔPaKmt6 and ΔPaHP1 mutants. ChIP_peaks, 0 means no peaks detected by MACS2, 1 means at least one peak detected by MACS2. ChIP_coverage: mapped reads after normalization (see “ChIP-seq sequencing data analysis and visualization” paragraph in “Methods” for details).

Additional file 22: Table S4.

Cq for RT-qPCR experiments and geNorm analysis of candidate normalization genes.

Additional file 23: Table S5.

Primers used for PCR experiments. A Primers used for RT-PCR and allele construction experiments. B Primers used for RT-qPCR experiments.

Additional file 24: Table S6.

Accession numbers for proteins used in alignments to build the phylogenic trees.

Additional file 25: Table S7.

Heatmap of Spearman correlation coefficient comparison between ChIP-seq samples. Samples are labeled as follow: Genotype_strain_marks. E.g., KMT6-11d_K9 = ΔPaKmt6 strain, from spore number 11d (biological replicate), IP = H3K9me3. Spearman’s correlation coefficient confirms the reproducibility of the experiments since biological replicates cluster together with high correlation coefficients. Moreover, H3K9me3 and H3K27me3 IP share a same cluster that negatively correlates with the H3K4me3 cluster. Both display no correlation with non-IP input and negative sample (Mock and H3K27me3 in the ΔPaKmt6), which means that IP results are significant and specific.

Additional file 26: Table S8.

Raw data from the heatmap of Spearman’s correlation coefficient comparison in Additional file 18: Fig. S18.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Carlier, F., Li, M., Maroc, L. et al. Loss of EZH2-like or SU(VAR)3–9-like proteins causes simultaneous perturbations in H3K27 and H3K9 tri-methylation and associated developmental defects in the fungus Podospora anserina. Epigenetics & Chromatin 14, 22 (2021). https://doi.org/10.1186/s13072-021-00395-7

Download citation

Keywords

  • Histone methylation
  • Sexual development
  • ChIP-seq
  • Podospora anserina