CC-BY 4.0
Introduction
Most eukaryotes rely on sexual reproduction to produce offspring. However, the maintenance of sex remains a fundamental question in evolutionary biology. The amazing diversity of fungi makes them of particular interest for addressing this question. Fungi are found in almost every terrestrial and marine ecosystems (Stajich et al. 2009; Breyer and Baltar 2023). Some play key roles in biogeochemical cycles, while others possess traits that enable them to be pathogenic to plants or animals (Iliev et al. 2024). Understanding the biology of fungi, and in particular why many of them can produce either asexual or meiotic products depending on environmental conditions (Dyer and Kück 2017; Wang et al. 2018), is crucial for fields as diverse as ecology, agronomy and human health.
For fungi, as for all eukaryotes, sexual reproduction has well-documented advantages. In particular, meiosis, the sex-specific division, reduces the accumulation of deleterious mutations and can link beneficial alleles that are subsequently fixed in populations by natural selection (Otto 2009). In harsh environments, sexual reproduction, as opposed to asexual reproduction, can lead to improved fitness of the offspring.
Fungal sexual reproduction relies on a wide range of developmental strategies and mating systems. These include both inbreeding, in which a single individual can self-fertilize, and obligate outbreeding, in which only genetically distinct compatible partners can mate. Some fungi, referred to as pseudohomothallic, have managed to retain the best of each strategy. Pseudohomothallism extends the diversity of reproductive systems evolved in fungi by selecting breeding systems that allow self-sterility to be bypassed by maintaining two populations of compatible nuclei in the same mycelium (Peraza-Reyes and Malagnac 2016). Early forward genetic screens and more recent transcriptomic studies in Sordariomycetes (Wang et al. 2014; Lütkenhaus et al. 2019; Kim et al. 2022) have led to the identification of large sets of genes involved in specialized developmental pathways that control the fertilization and the subsequent developmental stages (Wilson et al. 2019). Some of these species were heterothallic (self-incompatible, self-sterile (Wang et al. 2014; Lütkenhaus et al. 2019)), while others were homothallic (self-compatible, self-fertile (Teichert et al. 2012; Wang et al. 2019)). However, how the development of meiotic tissues within the fruiting body is controlled and which genes are critical in this process remains unknown.
In this study, we propose to approach this question from a different angle using the pseudohomothallic species Podospora anserina. This species is used as a model system for the study of the fruiting body development (Peraza-Reyes and Malagnac 2016) but a detailed analysis of the cellular programs allowing its formation from fertilization to ascospore discharge was lacking. In this context, we performed a genome-wide gene expression profiling during the sexual reproduction of P. anserina at ten time points representing key developmental steps. To characterize the topology of differential gene expression profiles across sexual reproduction, we analyzed our transcriptomic data using two complementary approaches: i) a data-driven approach using unsupervised clustering and ii) a biological expertise-driven approach to mine the gene functions. The data-driven approach revealed five major patterns (waves) of differential expression during sexual development. The expert annotation approach has led to the identification of a number of interesting gene categories, including secondary metabolites and toxins encoding genes, genes responsible for fungal vegetative incompatibility programmed cell death, and gene related to epigenetic effectors that interfere with messenger stability (RNAi/MSUD genes). Finally, we identified several master transcriptional regulators and draw their potential target gene network. Overall, our study has produced a comprehensive database that offers numerous candidates for functional genomics experiments. This database also paves the way for the characterization of new pathways during sexual reproduction.
Materials and methods
Podospora anserina strains and methods
P. anserina is a filamentous ascomycete whose life cycle and general methods have already been described (Esser 1974; Rizet and Engelmann 1949; Zickler et al. 1995). All strains used in this study were derived from the S strain (Rizet and Engelmann 1949; Espagne et al. 2008). Growth and crosses of S mat+ and S mat- strains were carried out on a minimal synthetic medium (M2 medium) at 27°C with constant light (Silar 2020). The 136 strain is described in (Coppin and Silar 2007).
Cytology
Sexual cells were either fixed in 7.4% paraformaldehyde and processed for microscopy as previously described (Thompson-Coffe and Zickler, 1994) or directly observed in M2 medium. When required, nuclei were stained with DAPI (4’,6-diamidino-2-phenylindole, 1 mg/mL; Boehringer Ingelheim). Observations were performed with a Zeiss Axioplan microscope, and images captured with a CDD Princeton camera system.
Time-course procedure and RNA isolation
To facilitate biological materials collecting, we used cellophane sheets (cat # 1650193, Bio Rad Laboratories, Hercules) and cheesecloth pieces (Nitex 03-48/31, Sefar AG, Heiden), which were deposited on M2 medium prior to inoculation. Cellophane sheets avoid the contamination of biological material by M2 medium when scraping mycelium. The mesh of cheesecloth accommodates the mycelium, while perithecia develops above the mesh and can be scraped with minimal mycelium contamination. S mat- strain has been used as male gamete donor and S mat+ strain as female. The S mat+ strain was grown on cellophane sheet for the analysis of mycelium and on cheesecloth for analysis of perithecial development. The biological material was then collected 6, 12, 18, 24, 30, 42, 54 and 96 hours after fertilization and an additional time point 24 hours prior fertilization was taken. Three biological replicates issued from three independent time course have been collected for all time points and two for T42, as we discarded the third and the back-up replicates, which did not meet our quality requirements.
For each ten points of the time course, the number of Petri dishes has been adjusted (five were used for T-24h, T0h, T6h, T12h, T18h, T24h, four for T30h, two for T42h, one for T54h and T96h) to collect 20 to 100 mg of fertilization competent mycelium or/and perithecia. One to five Petri dishes per time course experiment were thus inoculated simultaneously with a growing S mat- strain and spermatia were collected by washing with 1.5 ml of H2O six days later. The collected spermatia (around 104 spermatia per ml) were pooled and then used to fertilize a corresponding number of Petri dishes of 4-day-old S mat+ strain. The moment just prior the fertilization process was considered as the T0 (0 hours).
The biological material collected by scraping gently either cellophane or cheesecloth with a cover glass, was dried out with filter paper, weighted, frozen in liquid nitrogen and stored at -80°C. Total RNAs of P. anserina were extracted using RNeasy Plant Mini Kit (Qiagen), including a grinding process using a Mikro-Dismembrator (Sartorius) and a DNase treatment. The quality and quantity of the total RNAs was determined by using a Nanodrop spectrophotometer (Nanodrop Technologies) and the Bioanalyzer 2100 system (Agilent Technologies) as previously described in Imbeaud (2005).
Gene expression array-based hybridization
Transcriptome microarray experimental procedure, i.e. targets preparation, hybridization and washing, was done following the two-color microarray-based gene expression analysis instructions (version 5.0, February 2007) as described by the manufacturer (Agilent Technologies). One-μg aliquots of total RNA were labelled using the Low RNA input fluorescent linear amplification (LRILAK) PLUS kit (Agilent Technologies); internal standards came from the Two-color RNA spike kit (Agilent technologies). The labelling efficiency and the product integrity were checked as described by (Graudens et al. 2006). The reference cRNA was labelled with Cy5 and sample probes with Cy3. Briefly, the reference RNA pool was prepared by mixing equal amounts of RNA extracted from mycelium at 48 h and 96 h post-inoculation, and from perithecia at 24 h, 48 h and 96 h post-fertilization (Bidard et al. 2010). This composition allows a detectable and stable expression signal on the majority of CDS (Bidard et al. 2010). Then, Cy3- and Cy5-labeled targets were mixed and incubated on an Agilent microarray slide for 17 hours at 65°C, in a rotating oven (6 rpm for a 1x44K and 10 rpm for a 4x44K array format), using an Agilent in situ hybridization kit. The slides were washed and then any traces of water were removed by centrifugation at 800 rpm for 1 min (Bidard et al. 2010, 2011).
Data acquisition from microarray experiments
Microarrays were scanned using an Agilent dual laser DNA microarray scanner (Agilent Technologies), model G2567AA, with 5-µm resolution. The resulting 16-bit image files were analyzed in the Feature Extraction processing system (FE, v9.5.3). Spot and background intensities were extracted and normalized with the Feature Extraction software using the GE2-v4_95_Feb07 default protocol (Local background subtracted and LOWESS normalization). Preliminary array quality was assessed using Agilent control features as well as spike-in controls (Agilent 2-Color Spike-in Kit for RNA experiment). Subsequent flagging was done according to the GenePix Pro software (Molecular Devices Syunnyvale, CA, USA) nomenclature, including four levels of flags (good (100), bad (-100), not found (-50), moderate (0)) and raw data normalization was performed as previously described in articles such as (Bidard et al. 2010, 2011). Briefly, it consisted in using standard bioinformatics protocols to first filter probes with good quality flags, then correct signal intensities considering global and local background noises, and finally correct for systematic bias between Cy3 and Cy5 signals.
Graphical summary of data analysis
A schematic representation of the obtained data set is shown in Supplementary Figure 1 (see number #1). Gene expression measurements were thus analyzed according to two different strategies, allowing to study either the co-expression between any pairs of annotated genes in the P. anserina genome (see number #2) or the differential expression for all annotated genes in the P. anserina genome between any pairs of time points (see number #3). For this, gene expression profiles were calculated and statistical analyses to search for differentially expressed genes were performed (see the following sections for more details).
Calculation of gene expression profiles
For each gene, “A values” were computed at each time point and for each replicate, with the following formula: \(A_{g,t,r} = \frac{1}{2}\left( \log_{2}\left( G_{g,t,r} \right) + \log_{2}\left( R_{g,ref} \right) \right)\), where \(G_{g,t,r}\) represents the normalized Cy3 intensity value measured for gene \(g\) at time point \(t\) in biological replicate \(r\), and \(R_{g,ref}\) represents the normalized Cy5 intensity value measured for gene \(g\) in the reference RNA pool (see previous sections). At each time point, triplicated values were next averaged to obtain one value per time point, noted \(A_{g,t}\). The mean and the median values, respectively noted \(A_{mean}\) and \(M_{g}\), and associated with \(A_{g,t}\) values over all time points, were finally calculated and used to adjust expression measurements for each gene around a new median value of 0, i.e. \(A_{g,t}^{*} = A_{g,t} - M_{g}\). An expression profile for a gene \(g\) is therefore made with ten \(A^{*}\) values which are: \(A_{g,Tm24}^{*};A_{g,T0}^{*};A_{g,T6}^{*};A_{g,T12}^{*};A_{g,T18}^{*};A_{g,T24}^{*};A_{g,T30}^{*};A_{g,T42}^{*};A_{g,T54}^{*};A_{g,T96}^{*}\). Notably, it represents the expression fluctuations of the gene over the time course.
Statistical analyses to search for differentially expressed genes
Gene expression measurements between pairwise time points were compared using the MAANOVA (MicroArray ANalysis Of Variance) package (version 1.14.0) (Wu et al., 2003). It allows fitting an ANOVA model as proposed by (Kerr and Churchill 2001) and provides a way to consider multiple sources of variation in the microarray experiments. P-values were computed by permutation analysis and were adjusted by the false discovery rate step-up method of Benjamini and Hochberg to account for multiple testing (Benjamini and Hochberg, 1995). From all comparisons, lists of differentially expressed genes (Table S1) were defined using several filtering strategies:
Filter #1 to select the genes with (i) an adjusted p-value < 0.05, (ii) a fold change ≥ 4 or ≤ -4 with a \(A_{mean}\) ≥ 4.5,
Filter #2 to select the genes with (i) an adjusted p-value < 0.001 and (ii) a fold change ≥ 2 or ≤ -2 (i.e. absolute value of log2(fold change) ≥1).
The two lists of genes were used respectively for unsupervised clustering (see below) and biological interpretations (see the Results section).
PCA analysis and unsupervised clustering based on gene expression profiles
Expression profiles of all genes were used for principal component analysis using the function “pca” available in the R package mixOmics (Table S2) and expression profiles of differentially expressed genes (selected based on Filter #1, see previous section) were used for cluster analyses. To this list, the expression profiles of genes coding for mating-type transcription factor SMR1, SMR2 and FMR1 were added. A hierarchical clustering (UPGMA method) and several K-means clustering were computed independently using Spotfire DecisionSite software package (Spotfire, Somerville, Mass.) with correlation coefficient as similarity measure between \(A^{*}\) values in gene expression profiles. Intersection of clusters obtained by the two-clustering method results in robust clusters (Table S1) from which the temporal waves (see Result section) were derived.
Sequence comparisons to search for orthologous genes with other species
Orthologous genes were searched between P. anserina and the other species Neurospora crassa, Sordaria macrospora, Chaetomium globosom and Trichoderma reesei. For this purpose, protein sequences were downloaded from Uniprot (https://www.uniprot.org/) under the accession numbers: UP000001805_367110 (N. crassa), UP000001881_771870 (S. macrospora), UP000001056_306901 (C. globosum) and UP000008984_431241 (T. reesei). The OrthoFinder algorithm (Emms and Kelly 2019) was used with default parameters to obtain “orthogroups” (Table S3), i.e. groups of genes that are good candidates to share a common ancestor between the compared species.
Functional enrichment of terms associated with lists of genes
Gene Ontology (Ashburner et al. 2000; The Gene Ontology Consortium et al. 2023) and Pfam annotations (Mistry et al. 2021) were used to explore the biological relevance of the lists of genes that composed the temporal waves (Tables S4 and S5). P-values were calculated using Fisher’s exact test and following the methodology described in (Boyle et al. 2004).
Expert annotation of P. anserina genes
Expert annotation (Table S6) was performed using either BLASTP reciprocal best-hit and/or Pfam conserved motif analyses. Output results were all validated after careful manual inspection including the use of the FungiBD database for orthologs verifications (Basenko et al. 2018).
Co-expression network based on expression profiles of genes coding for transcription factors
Expression profiles of differentially expressed genes (selected based on Filter #2, see previous section) for which the expert annotation was associated with “transcription factors” (Table S7) were used to search for significant correlation (> 0.95) with the expression profiles of all genes previously classified into transcriptional waves. Selected gene pairs were used to represent a co-expression network as defined in (Wolfe et al. 2005), in which all potential links of transcriptional control, between the TF and another gene (referred to as “target gene”) are shown.
Availability of data and source code, website to explore the results
The raw data have been deposited in the MIAME compliant Gene Expression Omnibus database (Edgar 2002) and is accessible through the GEO Series accession number GSE21659 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE93094. The source code has been deposited at : https://github.com/Podospora-anserina/transcriptome_kinetics_development/). To explore the results for a particular gene or for a set of genes, we have created a WEB application which is available at: https://pierregrognet.shinyapps.io/cinetique_app/.
Results
High-resolution kinetics of transcriptional fluctuations during sexual reproduction
To ensure an accurate, high-density representation of transcriptional variation, we collected RNA samples at ten key time points during the sexual reproductive cycle of P. anserina (Figure 1). P. anserina is a pseudohomothallic fungus, namely natural isolates are self-fertile because they contain nuclei of each mating-type. Self-fertilization does not start at a specific time point and therefore natural isolates cannot be used for a time course analysis of sexual development. However, the genetic basis of the P. anserina mating system is heterothallism with two mating-type idiomorphs called mat+ and mat- (Turgeon and Debuchy 2007). Typical ascospores contain mat+ and mat- nuclei, but a few asci contain monokaryotic ascospores which produce self-sterile individuals. Those latter individuals were used in the experiments presented here. Fertilization was performed by spreading mat- male gametes on a mat+ strain. This experimental design allows the time course to start at a well determined time point for synchronization purpose. Under this experimental condition almost all fertilization events occurred within the first six hours (Figure S2A and S2B). However, despite this experimental design, RNA samples were extracted from biologically non-homogeneous tissues because P. anserina fructifications are composed of different cell types (Figures 1A and 1B).
We measured the expression of 10,507 nuclear genes at ten time points which are designated as Tm24 for mycelium 24 h before fertilization and T0, T6, T12, T18, T24, T30, T42, T54 and T96 for time 0 to 96 h postfertilization, respectively (Figure 1C, Table S1). We performed a Principal Component Analysis (PCA) on the complete dataset. More than 75% of the total variance was captured by only three principal components (PCs, Figure 2A and 2B, Table S2). With the exception of Tm24, we observed that all samples were distributed along the PC1 axis (38%) according to the sexual development time line (Figure 2A). On PC2 axis, we found T30 and T42 samples clustered apart from the others (24%). Interestingly, this period corresponds to the formation of the dikaryotic cells, which are the progenitors of fungal zygotic lineage (Figure 1). The peripheral position of the Tm24 sample (Figure 2B) which appeared on the PC3 axis (16%), indicated that the gene expression profiles derived from vegetative mycelium (Tm24) were clearly different from those of all samples collected after fertilization. Overall, this unsupervised exploratory data analysis validated both the relevance of our experimental design and the robustness of our differential analysis (see below). The expression profiles for the 50 genes whose expression measurements contributed the most to the calculation of the principal components 1, 2 and 3 are shown Figure 2C (and detailed Table S2). The data indicated that the greatest fold changes, both positive and negative, occurred after the formation of dikaryotic cells (>T30). This was also the case for PC2, but only after the formation of the ascospores (>T54). We observed that the positive fold change values on PC1 exhibited a higher dispersion after T42.
Figure 1 - P. anserina developmental time course from fertilization to ascospore maturation. A. Schematic representation of the fertilization and the main steps of zygotic development. B. Pictures of the main steps of sexual development. Sexual reproduction began with the production of male (spermatium) and female (ascogonium) gametes. A pheromone/receptor signaling system allows the female gametes to attract specialized hyphae (trichogynes) that recognize and fuse with male gametes of compatible mating type. This event initiates the multi-step process of sexual development that begins with the formation of fructifications and ends with the discharge of ascospores. The nucleus of the male gamete enters the trichogyne, migrates along this specialized hypha, reaches the ascogonium and undergoes several mitotic divisions. Fertilization (T0) triggers the development of protoperithecia (T6 to T24), which recruit multiple layers of protective maternal hyphae to grow into perithecia (T24 to T45). The developing fructifications shelter the fertilized ascogonial cells (T24 to T30), which contain multiple haploid parental nuclei of either mat+ or mat- genotype. Differentiation of paraphyses begins (T24) just before three-celled hook-shaped structures called croziers emerge from the ascogonial cells (T30 to T42). Croziers are dikaryotic because they contain two haploid nuclei derived from both parents. Karyogamy (T30 to T44) occurs in the upper cells of croziers. The resulting diploid zygote immediately undergoes meiosis (T30 to T44), producing four haploid nuclei, which then undergo mitosis prior to ascospore delineation (T54 to T96). Ascospores are formed around two non-sister nuclei within the developing asci. In rare cases, two ascospores form around a single haploid nucleus, resulting in five-spored asci. C. Time line of the sexual cycle of P. anserina with the steps from which RNAs have been extracted on the left.
Figure 2 - PCA analysis reveals a strong temporal effect on gene expression variations relevant to the key developmental steps in P. anserina. Principal Component Analysis (PCA) on the complete dataset (i.e. all gene expression profiles, see Methods). A. Axis 1 and 2 of the PCA explain 38% + 24 % = 62% of the variance. B. Axis 2 and 3 of the PCA explain 24% + 16 % = 40 % of the variance. Altogether, these first three axes of the PCA explain 78 % of the variance. C. Gene expression profiles for the 50 genes, which contribute the most to the definition of the first (PCA1), second (PCA2) and third (PCA3) principal components, as shown in A and B. Detailed lists of these genes are given in Table S2. They were separated according to the type of regulation they experience during sexual development (positive or negative). Tm24, T0, T6, T12, T18, T24, T30, T42, T54 and T96 : time points of the kinetics were the biological samples were collected, described in Figure 1.
Topology of global expression variation reveals five major transcriptional waves during sexual development
Using k-means clustering, we defined 13 representative expression patterns (Figure S1 and Figure S3), encompassing 1,186 genes. Further data integration outlined five successive transcriptional waves (Figure 3, Figure S1 and Table S1). Notably, they all shared a common pivot point (Figure S4), corresponding to the mid-course of the kinetics (T30-T42). This observation indicated a pause in the regulation of transcriptional activity at the time of dikaryon formation. The 1,186 genes were not evenly distributed across the five waves (Figure 3 and Table S1), resulting in two periods of intense and equal transcriptional fluctuation (waves II and IV) interspersed with three periods of moderate transcriptional fluctuation (waves I, III, and V).
Figure 3 - Temporal transcriptional waves over the course of sexual development. A. Expression profiles of the five major successive transcriptional waves (numbered in Roman figures) that take place during P. anserina sexual reproduction (Bottom timeline, see panel B for zoom). Number of differentially expressed (DE) genes present in each wave is indicated below the corresponding graph. Waves II (early up/late down pattern, 31.51%) and IV (early down/late up pattern, 36.37%) accounted for 67.87% of them (33.98% and 37.42%, respectively), while wave I (early down pattern, 7.41%), wave III (early-mid up pattern, 9.62%) and wave V (late up pattern, 15.02%) contained only 32.13% of them, respectively. B. Major developmental events according to the sexual reproduction timeline described in Figure 1.
The topology of the transcriptional waves therefore suggests that the genes of waves II and III are down-regulated after a peak of expression during protoperithecium development, while the genes of waves IV and V are up-regulated after perithecium development. These patterns are in accordance with the results obtained in N. crassa, during fruiting body development (Wang et al. 2014).
Evolutive conservation and functional annotation of differentially expressed genes during sexual reproduction
Across the entire time-course experiment, we found 3,466 genes that were differentially expressed (DE, see Methods). They represent 33% of the total gene set for which gene expression measurements were monitored (Table S1, Figure 4A), including 1,133 of the genes clustered in the five transcriptional waves (more than 95% of the clustered genes, Table S1 and Figure 3 and Figure S1). Overall, this analysis revealed that in P. anserina, one-third of the complete gene set is transcriptionally regulated during sexual reproduction. Of this, one-third follows one of the five transcriptional waves disclosed by the clustering analysis. To test for evolutive conservation of these DE genes, we built group of orthologs (see Methods), using four species of filamentous fungi, i.e. the homothallic species Chaetomium globosum and Sordaria macrospora and the heterothallic species N. crassa and Trichoderma reesei. A total of 2,957 DE genes (more than 85%) were identified as belonging to an ortholog group shared with at least one of the four species (Table S3). A core set of 1,496 DE genes is conserved in all four fungi (Figure 4, Table S3).
Figure 4 - Distribution of differentially expressed genes along the kinetics of sexual reproduction and their evolutive conservation. A. Heatmap of pairwise comparisons of differentially expressed (DE) genes across the complete time-course experiment. Numbers correspond to the numbers of genes with an adjusted p-value lower than 0.001 and an absolute fold change value between the two compared time points, higher than 2 (see Methods). B. Venn diagram showing the conserved orthologs between the indicated fungal species. Cg: C. globosum (homothallic), Sm: S. macrospora (homothallic), Nc: N. crassa (heterothallic) and Tr: T. reesei (heterothallic).
A pairwise comparison of the DE gene set of P. anserina with each gene set of the four fungi, showed that the number of conserved DE genes is consistent with fungal evolutionary history. Closely related species (C. globosum) share more DE genes than the more distantly related species (N. crassa, S. macrospora and T. reesei). The remaining 506 genes (14.78%) were classified as orphan genes, with 247 (21.84%) belonging to the five transcriptional waves. Notably, the set of DE genes, including that of the five waves, were found to be enriched with orphan genes (p-values: 2.55x10-9 and 4.60x10-20 respectively).
We then search for functional annotation of the proteins encoded by the DE genes (see Methods). We found that 1,279 (36.95%) of the 3,466 putative proteins encoded by the DE genes lacked Pfam annotation (Table S1, (Mistry et al. 2021)), which is significantly more (p-value: 8.55x10-9) than in the complete list of P. anserina’s proteins. The lack of annotated conserved domains was further enriched (p-value: 1.48x10-13) in wave-embedded DE genes, as 515 out of these 1,133 putative proteins (45.44%) lacked Pfam annotation. We also performed Gene Ontology (GO, (Ashburner et al. 2000; The Gene Ontology Consortium et al. 2023)) analysis on the DE genes (Table S4, Figure 5A).
Figure 5 - GO term and expert annotation categories enriched in the transcriptional waves. A. Descriptions of GO terms significantly enriched (p-value ≤ 0.05) in at least one wave of the transcriptional waves. B. Enrichments of genes of interest in at least one wave of the transcriptional waves. TF: transcription factors, SM: secondary metabolites, HET: fungal specific vegetative incompatibility HET domain. Points are scaled and colored according to the number of genes in each category for each wave. No score is shown for zero genes.
In general, the enriched GO terms and Pfam domains (Tables S4 and S5) described intense catalytic (GO:0003824) and metabolic (GO:0008152) activities, associated with binding (GO:0005488, PF00187). Accordingly, we observed enrichments in genes involved in redox catalytic activities (GO:0003824, PF00067, PF03055, PF01565, PF00106, PF01494, PF00199, PF08240) and transporter activities (GO:0005215, Pfam PF07690.19). We also uncovered enrichments in DE genes involved in cellular differentiation programs more specifically related to sexual reproduction. They either deal with cellular anatomical entity (GO:0110165), localization (GO:0051179) or encode zinc-binding domain transcription factors (PF13695), necrosis inducing factor (PF14856) and TUDOR domain containing protein (PF11160). However, if the repertoire of GO terms (Table S4) and Pfam domains (Tables S1 and S5) generated by this study should provide some clues to select candidate genes for in vivo functional characterization, our analysis essentially showed that more than one third of the DE genes have no associated function yet.
Searching for genes of interest, the expertise-driven approach
We then performed a biological expertise-driven approach to mine the functions of the full set of DE genes (Table S6, Figure 5B).
Developmental processes at work
Fertilization of P. anserina requires exposure to blue light, but the molecular basis of this signal transduction is still unknown. In N. crassa, the blue light signal is transduced by the photoreceptors WC-1 and WC-2 (Froehlich et al. 2002; He et al. 2002). In this study, we found five blue light responsive genes that were differentially expressed (Table S6A, Figure S5A-B). Consistent with a role in fertilization three of them (including Pa_ 5_6130/WC-1) were expressed before meiosis, while the other two were expressed after meiosis, suggesting that blue light may also play a role during ascospore formation.
Orthologs of genes involved in plasmogamy in yeast models (Pa_3_11320/Prm1 (Curto et al. 2014), Pa_3_4770/CDC3 (Kurahashi et al. 2002), Pa_5_9510/Fig1 (Erdman et al. 1998), Pa_1_1930/KEX2 (Julius et al. 1984), Pa_4_6420/KAR2 (Rose et al. 1989), Pa_7_10810/SEC63 (Young et al. 2001)) were found to be up-regulated until T42 and then downregulated (Table S6A, Figure S5C). Development of the zygotic lineage likely involves compartmentalization through active cell wall remodeling, as we found several DE genes encoding different types of cell wall proteins and growth factor receptors (Table S6A). One of them was the expansin-like EEL2, involved in cell wall expansion, whose gene Pa_2_310 was downregulated from karyogamy to ascospore maturation (T30-T96). This expression pattern was similar to that of the C. globosum homolog (CHGG_00523) but divergent from that of the N. crassa homolog (NCU04603) (Hutchison and Glass 2010; Wang et al. 2019). Cell-cell interactions may also be tightly regulated, as exemplified by the strict downregulation (T0-T96) of Pa_1_17480, which encodes the adhesin protein MAD1 (Wang et al. 2007).
Since sexual development is rich in cell divisions, we searched for DE cell cycle related genes and found this category slightly enriched (p-value = 2.9x10-2, Table S5B). In addition to the TOR pathway inhibitor Pa_3_3700/FgFkbp12, it included homologs of genes required for proper fructification formation (Pa_1_17170/KIN3, Pa_1_23300/cel-2, Pa_2_6460/Pah1, Pa_2_890/asm-1), axial budding (Pa_2_290/axl2) and genes encoding G-coupled receptors (Pa_4_5350, Pa_1_590) and ras-like protein (Pa_1_19030/krev-1). The expression profile of the Pa_2_290/axl2 gene, encoding an integral plasma membrane protein, was of particular interest: it was down-regulated when the two haploid parental nuclei actively divide into multinucleated syncytial ascogonial cells (Figure S5D).
As expected, we also found DE genes encoding cytoskeletal elements (Figure S5E). Interestingly, Pa_5_5720, encoding a pyruvate decarboxylase (PDC) was found strongly up-regulated during ascospore maturation (wave V). This observation further enforces the hypothesis that pyruvate decarboxylase is involved as a structural protein (PDC-filaments) associated with the cytoskeleton during the sexual development of filamentous fungi (Thompson-Coffe et al. 1999). Conversely, Pa_4_7260, encoding the nucleoprotein Cik1, required for both karyogamy and mitotic spindle organization (Cottingham et al. 1999), was found down-regulated at the same time. Some of the cytoskeleton DE genes showed co-expression with the plasmogamy DE genes (Figure S5F), making them good candidates to be involved in fusion morphogenesis processes.
Three genes essential for ascus development were also found DE (Figure S5G, Table S6A): P. anserina’s orthologs of Pa_4_9730/asd-1 (wave IV (Nelson et al. 1997)), Pa_4_9750/asd-3 (Galagan et al. 2003) and Pa_2_13270/con-6 (wave V (White and Yanofsky 1993)). P. anserina does not form clonal propagules. However, the orthologs of genes required for conidiation in N. crassa and Aspergillus nidulans such as Pa_2_13270/con-6, Pa_4_9730/asd-1, Pa_4_9750/asd-3, Pa_1_12720/Pyg1, Pa_1_30/TmpA and Pa_0_1540 were found to be up-regulated after meiosis, suggesting that genes primarily involved in conidiation, i.e. asexual reproduction, may have been recruited to participate in ascospore formation. On the contrary, early conidial development-2 Pa_7_780/ecd-2 was down regulated during the first half of the kinetics (Sun et al. 2012).
Producing energy
The reproductive process is energy consuming and therefore the mobilization of carbon resources is one of the keys to its success. The genome of P. anserina contains 276 genes encoding CAZymes (Espagne et al. 2008; and Table S6D). Among them, 141 were found DE during sexual development, showing an enrichment of this metabolic activity (p-value = 4.36x10-5). The expression profile of 42 CAZyme genes corresponded to one of the five transcriptional waves, especially waves II and IV. The remaining 99 showed contrasting expression patterns (Figure S5 H-T), but up-regulation was enriched over down-regulation, as 62 of them were up-regulated during the kinetics (Chi², p-val < 0.001). We also found that several CAZymes genes had some of the highest positive fold changes measured in this study (Table S1). In addition, three glycolysis enzymes genes (6-phosphofructo-2-kinase Pa_3_6610, fructose-bisphosphate aldolase Pa_0_710, and phosphoglycerate kinase Pa_6_10610) were up-regulated with a similar pattern of overexpression when asci begin to differentiate (Figure S5U, Table S6E).
In P. anserina, long-chain fatty acids are shortened in peroxisomes then degraded in mitochondria. The Pa_6_3500/PaPex7 gene, which encodes the peroxisomal-targeting-signal-2 receptor, is found to be up-regulated during ascospore maturation (T54-T96, Table S1 and Figure S5V). This is in line with the defects in nuclear positioning of the post-meiotic asci observed in pex7 knock-out mutants, which can result in either abnormal ascospore formation or to ascospores containing an abnormal number of nuclei (Bonnet et al. 2006). The Pa_1_23790/PaFox2 gene (Table S1 and Figure S5V), which encodes the peroxisomal multifunctional protein 2, a key enzyme in the fatty acid ß-oxidation pathway, was found to be down-regulated during the early stages of sexual development (T12-T54) and up-regulated during ascospore maturation (T54-T96). This expression profile may explain the reduced acetyl-CoA availability within the ascospores that results in the pigmentation defect observed in PaFox2 knock-out mutants (Boisnard et al. 2009).
Secondary metabolites: major actors of sexual development?
Prokaryotic and eukaryotic micro-organisms produce a wide variety of secondary metabolites (SMs). These molecules are not intended to perform housekeeping functions, but rather are used as adaptive tools to cope with environmental fluctuations or to compete with other living organisms. The P. anserina genome contains 35 biosynthetic gene clusters (BGCs) (Lamacchia et al. 2016), (Table S6F, Figure S5W), containing 474 genes. SM encoding genes were enriched among the complete set of DE genes (p-value = 8.31x10-10) as 258 of them were found to be DE (33 out of 35 BGCs contained at least two DE genes). These genes were also significantly up-regulated (Chi2, p-val < 0.001). BGCs N°1, 11, 18, 19, 24, 27 and 34 showed an expression profile topology (Figure S5W) similar to the general pattern described in N. crassa (Wang et al. 2022b), which was a coordinated downregulation from fertilization to meiosis followed by a coordinated upregulation during ascospore formation and maturation.
We also examined whether co-expression was shared by multiple SM genes within the 34 BGCs. This was the case for the sterigmatocystin BGC, where 17 of the 23 genes were differentially expressed during sexual reproduction (Table S6F, cluster 32, Figure S5X, Slot and Rokas (2011)), including the transcription factor PaAflr (Shen et al. 2019). In addition, we identified four positive regulators of the sterigmatocystin BGC (Table S6F, Figure S5Y, Pa_2_7220, Pa_2_7240, Pa_2_7250, Pa_2_7260 (Bhetariya et al. 2011)), that showed the same expression pattern across sexual reproduction. Similarly, 13 of the 20 genes in the terrequinone BGC were differentially expressed (Figure S5Z).
In addition, we identified several isolated genes from different BGCs that showed similar peak expression to each other but no co-expression with other genes from their respective BCG (Figure S5AA-AB). These expression peaks of genes encoding mostly proteins of unknown function corresponded to two key moments in early development (T6 and T12), just before and during dikaryon formation. Taken together, our analyses further support the key role of SM in sexual reproduction and highlight the role of previously uncharacterized clusters and isolated SM genes.
From regulation of biotic interactions to that of sexual development
Genes encoding STAND proteins (Signal-transducing ATPases with Numerous Domains, (Leipe et al. 2004)) have been described as involved in fungal biotic interactions with living organisms, either from the same species or from different species but mostly during vegetative growth. The P. anserina genome contains 34 genes encoding STAND proteins with a NACHT domain (PF05729), eight of which were found to be down-regulated either throughout the sexual cycle kinetics or prior to perithecium development (Table S6G, Figure S5AC). This group includes Pa_3_10930/nwd1 (Chevanne et al. 2010), Pa_7_10370/sesB-like, Pa_2_7940/het-d (Espagne et al. 2002) and Pa_1_11380 (Bidard et al. 2013). Fungal interactions are also under the control of a programmed cell death called vegetative incompatibility that prevents hyphal fusion of two genetically distinct wild-type strains. Most of the effectors of this incompatibility reaction (IR) contain HET domains (PF06985). Strikingly, 59 of the 128 het genes annotated in the P. anserina’s genome were found to be DE, including Pa_4_1190/hnwd1 and Pa_7_3460/Pin-C (Chevanne et al. 2010; Kaneko et al. 2006; Bidard et al. 2013). The majority of the DE het genes were down-regulated either throughout the complete sexual development or from dikaryon formation (Figure S5AD), but seven of them were up-regulated during some steps of the kinetics, mostly after karyogamy (Figure S5AE). In addition, three het genes showed a peculiar behavior (Figure S5AF): Pa_5_2540 expression oscillated during sexual reproduction, Pa_3_5550 expression peaked during protoperithecium development and Pa_4_1800 expression was down-regulated until ascospore maturation and over-expressed during this process.
In the DE gene set, we also found gene families encoding toxins that can induced host cell death or effector-like proteins. All of the six genes encoding the NEP1-like proteins pore-forming toxin NPP1 (necrosis-inducing Phytophthora protein, (Ottmann et al. 2009)) were down-regulated at least during the second half of the sexual development, i.e. before karyogamy (Figure S5AG). Six out of the eight genes encoding Hce2 pathogen effectors (Stergiopoulos et al. 2013) and one gene annotated as Ptu bug toxin encoding gene (Pa_6_9450) were also DE. In contrast to the NPP1 genes, all of them showed up-regulation during the course of sexual development, five of them belonging to the wave IV pattern of expression (Figure S5AH).
Epigenetic effectors at work in sexual reproduction
In recent years, epigenetic effectors have been identified as key players in fructification morphogenesis (Carlier et al. 2021; Nowrousian 2022). Therefore, we investigated the expression patterns of genes encoding proteins that can either methylate cytosines, write or erase histone post-translational modifications or interfere with messenger RNA synthesis and stability.
PaRid and PaDnmt5, two of the three DNA methyltransferase encoding genes present in the P. anserina genome, were found to be DE during this time-course experiment. PaRid expression, which is essential for dikaryon formation (Grognet et al. 2019) peaked at T12, while PaDnmt5 expression was increasingly up-regulated from T12 to T42 (Figure S5AI). Pa_5_2900, which encodes a RING-type E3 ubiquitin transferase URF1-like, was found to be strongly up-regulated from T42. In mammals, this protein binds both hemi-methylated DNA and tri-methylated lysine 9 of histone H3 (H3K9me3), thereby bridging DNA methylation and chromatin modification (Foster et al. 2018).
Pa_7_10020/KAT1 and Pa_1_7460/Hos-2 encodes two chromatin modifiers with antagonistic enzymatic activities, respectively a histone acetyltransferase encoding gene (Tong et al. 2012) and a histone deacetylase (Wang et al. 2002a) respectively. However, these two genes showed co-expression over the entire kinetics, i.e. they were up-regulated from fertilization to karyogamy and then down-regulated (Figure S5AJ). This pattern was in contrast to that of Pa_7_5890 encoding the Sin3-associated polypeptide SAP18, which enhances transcriptional repression through histone deacetylation (Cheng and Bishop 2002). Two genes encoding sirtuins NAD(+)-dependent histone deacetylases were also found to be DE: Pa_4_8570/PaSir2 showed a wave III expression pattern that resembled that of Pa_1_19440/PaRid (Figure S5AI), while the gene encoding Pa_4_7390/SIRT5 was down-regulated throughout sexual reproduction. Of the 32 genes encoding SET domain proteins, seven were found to be DE (Figure S5AK), although none has yet been shown to have histone lysine methyltransferase activity. Four genes encoding putative histone demethylase containing JmJC domains were also found DE (Figure S5AL).
Promoter accessibility also depends on nucleosome positioning. The ATP-dependent chromatin remodeling protein SWI/SNF Pa_5_4310/SMARCA3 (Oh et al. 2013) was one of the most highly up-regulated wave III genes (Figure S5AM). Its transcription pattern is similar to that of Pa_1_19440/PaRid (see Figure S5AI) and Pa_4_8570/PaSir2 (Figure S5AJ). Transcription requires a large number of proteins assembled into conserved complexes. Two genes (Pa_1_17860 and Pa_6_1460) encoding components of the SAGA complex were found to be up-regulated from T6 to T54 with very similar expression patterns (Figure S5AM). However, we also found similarly co-regulated genes although encoding factors with antagonistic effects on transcription, such as the SWD2 subunit of the COMPASS complex (Pa_7_11070) and the extra sexcombs subunit of the Polycomb Repressive Complex 2 (Pa_3_4080).
There is no experimental evidence that post-transcriptional RNAi mechanisms Quelling (Romano and Macino 1992) and MSUD (Shiu et al. 2001) are functional in P. anserina. However, genes encoding essential components of Quelling (Figure S5AN), the RNA-dependant RNA polymerase (RdRP) Pa_7_4790/Qde1 (Cogoni and Macino 1999) and the dicer Pa_6_6150/Dcl-1 (Yang et al. 2015), were downregulated from T18 and T30, respectively. On the contrary, expression of Pa_7_9210 encoding the Qde-2 interacting exonuclease QIP was up-regulated from T6 to T30. N. crassa ∆qde-1 quelling mutant strains produced normal-sized perithecia but no ascospore (Wang et al. 2014). Four genes involved in MSUD were also found to be DE (Figure S5AO). Pa_2_9830/Sad-2 and Pa_5_11380/Sad-5 were up-regulated whereas Pa_1_19900/Sad-1, encoding the RdRP and Pa_1_16040/Sad-3, encoding the helicase, were down-regulated. Interestingly, Pa_6_11690/rrp-3, which encodes an RdRP that could not be clearly assigned to Quelling or MSUD pathways in N. crassa (Borkovich et al. 2004), was co-expressed with the Quelling genes Pa_1_19900/Sad-1 and Pa_1_16040/Sad-3 in P. anserina. Knockout of the N. crassa Sad-1 and Sad-3 resulted in arrest of ascus development in meiotic prophase and production of normal-sized perithecia and young asci without ascospores, respectively (Wang et al. 2014). Sad-2 is essential to N. crassa fructification development (Wang et al. 2014). Our results suggest that RNAi genes may be silenced during meiosis and sporulation, which may explain the apparent lack of MSUD in P. anserina but also that Pa_2_9830/Sad-2 and Pa_5_11380/Sad-5 might have been recruited to play a role during ascospore dormancy.
Repeated sequences present in both haploid parental genomes of P. anserina’s dikaryotic cells are subjected to cytosine-to-thymine mutations by repeat-induced point mutation (RIP, (Bouhouche et al. 2004)). Modified 5-methylcytosines are known to be hotspots for cytosine-to-thymine transitions, because, unlike unmethylated cytosines, spontaneous deamination of 5-methylcytosines cannot be repaired by uracil excision (Pearl 2000). To counteract this genetic threat, members of the AlkB alpha-ketoglutarate-dependent dioxygenase (PF13532) family oxidize 5-methylcytosines to restore undamaged bases (Bian et al. 2019). Two genes encoding AlkB-like enzymes are present in the P. anserina genome. The expression of Pa_1_12280 was down-regulated throughout the kinetics, while Pa_1_17260 was up-regulated during two segments of the kinetics: i) from Tm24 to T12, corresponding to the peak of expression of PaRid and ii) from T24 to T54, when PaRid is no longer expressed (Figure S5AP). Notably, both genes showed a superimposable flat expression profile at the time of RIP.
In filamentous fungi, a conserved post-transcriptional A-to-I mRNA editing process recodes transcription units (Teichert et al. 2017; Blank-Landeshammer et al. 2019), mainly during fructification and ascospore morphogenesis. We identified seven DE genes encoding adenosine deaminase acting on tRNA (ADAT), a class of enzymes which are responsible for A-to-I mRNA editing. Pa_1_21170, a wave III gene, showed a peak of expression similar to that of PaRid, while Pa_3_6600 and Pa_6_7750 were up-regulated from perithecium development to ascospore discharge (Figure S5AQ).
Co-expression Network: In Search of Potential Transcription Factor Targets During Sexual Reproduction
We then focused on the 96 TFs encoding genes that were found to be DE (Table S6I). The expression profile of each TF gene was used to search for significant correlations with the expression profiles of all genes previously classified as transcriptional waves I, II, III, IV or V (see Methods). Using this strategy, we build a network of 2002 potential transcriptional links between 76 TFs and 683 target genes (Figure 6 and Table S7A-C).
Figure 6 - Co-expression network based on expression profiles of differentially expressed (DE) genes coding for transcription factors. Genes encoding transcription factors (TF) are represented by beige circles whose size represents the number of co-regulated genes. These coregulated genes, which can be considered as potential TF target genes, were classified in transcriptional waves as such: wave I: grey color (61 genes), wave II: red color (1402 genes), wave III: orange color (58 genes), wave IV: green color (426 genes) and wave V: blue color (55 genes).
Some of the prominent nodes of this network were expected TFs (FMR1, SMR1, SMR2, Pa_1_10140/Pro1, see below), confirming its relevance. However, most of the 76 TFs (72%) had no characterized biological function, including those with the broader set of target genes. The distribution of the TFs encoding genes in the five transcriptional waves differed from that of the total set of DE genes (Figure 6), with the majority belonging to wave II. As with TFs, most of the target genes were of unknown function. Importantly, some TFs were also target genes, e.g. Pa_1_10140/Pro1 is a target of Pa_1_12070/HOM3, Pa_1_11050/HMG4, Pa_1_17610, and Pa_2_13840 (Table S7A).
Among the master regulators of this network, we identified HMG-domain and homeodomain TFs (Espagne et al. 2008; Aït Benkhali et al. 2013). The HMG-domain TFs include the mating-type proteins FMR1 and SMR2, which are both required for the development of the perithecium after fertilization (Zickler et al. 1995; Arnaise et al. 2001). Physical interaction between FMR1 and SMR2 was demonstrated and proposed to be required for their post-fertilization functions (Arnaise et al. 1995; Turgeon and Debuchy 2007; Bidard et al. 2011). Our study shows that FMR1 and SMR2 have a similar expression profile throughout sexual reproduction (Figure S5AR). Their co-expression networks revealed also an almost identical number of potential target genes (FMR1, n = 79; SMR2, n = 78) (Table S7). They share 62 of their 79 potential target genes (Table S7A). These observations support the idea that they act together in identical regulatory pathways by forming a heterodimer (Arnaise et al. 1997). The third and last mat- mating-type gene, SMR1, does not belong to the HMG-box family. It encodes a HPG protein (Turgeon and Debuchy 2007), which has been proposed to provide a scaffold for physical interaction of mating-type proteins (Zheng et al. 2013). This role is consistent with its presence in wave IV, along with FMR1 and SMR2 (Figure S5AR). Critical genes for the development of ascogenous hyphae are likely to be found among target genes of FMR1, SMR2 and SMR1 (Figure 6, Table S7A). The only non-mating-type HMG protein present in the network is Pa_1_11050/PaHMG4, related to Nhp6p and Hmo1p from S. cerevisiae (Aït Benkhali et al. 2013). Its deletion in P. anserina resulted in an increased number of male organs. The lack of effect on fructification development is consistent with the down-regulation of Pa_1_11050/PaHMG4 at T42-T96 (Figure S5AS). Proper shaping of the developing fructifications requires the homeobox protein Pa_2_6460/Pah1 (Coppin et al. 2012). The expression of the corresponding gene was accordingly up-regulated from T0 to T42 (Figure S5AS). In contrast, Pa_3_4650/Pah6 was down-regulated between T0 and T54 (Figure S5AS), although no fructification defect was associated with its deletion.
The Zn(II)2Cys6 class of regulators is one of the largest families found in fungi (MacPherson et al. 2006) and consequently they represented almost half of the TFs present in the network (Table S7A). Three of them are of particular interest and have been extensively studied because of their properties. Pa_1_10140/PaPro1 is a regulator of the expression of the mating-type and IDC genes, both required for female fertility (Gautier et al. 2018). In S. macrospora, Pro1 has been identified as a master regulator of fructification formation (Steffens et al. 2016; Teichert et al. 2020). Its temporal and spatial expression is consistent with the presence of Pa_1_10140/PaPro1 in wave III (Figure S5AS). Pa_7_7000 was characterized in wave IV as a coregulator of FMR1, SMR2 and SMR1 target genes (Figure 6). Pa_7_7000 is revealed as one of the most prominent nodes in the TF network (Figure 6 and Figure S5AT), but has not yet been characterized in P. anserina. Pa_7_7000 has orthologs in the seven species of the Podospora complex (Boucher et al. 2017) and in some fungi closely related to this species complex (Cercophora samala, Triangularia verruculosa and Apiosordaria backusii). However, more distantly related species (Podospora setosa, Podospora australis, Podospora appendiculata, Thermothelomyces thermophilus, Neurospora tetrasperma and Neurospora discreta) contain a protein that shares some identity with Pa_7_7000, but they all lack the Zn(II)2Cys6 domain and are therefore unlikely to be functional orthologs. Wave V contains Pa_6_11860/Fluffy (Bailey and Ebbole 1998) (Figure S5AT), which is required in N. crassa for the formation of conidia. These vegetative spores are absent in P. anserina, suggesting that this gene has been recruited for ascospore formation in this fungus, as in others members of the P. anserina / P. pauciseta / P. comata species complex (Boucher et al. 2017).
Comparison of mating-type target genes at different developmental stages
A transcriptomic analysis of P. anserina was performed previously to determine the FMR1 and FPR1 target genes required for fertilization (Bidard et al. 2011). The analysis focused on the comparison of mat- and mat+ mycelia competent for fertilization, namely the stage preceding the development of the perithecium. It was performed using the same methods as described here. The analysis identified 69 genes upregulated in a mat- mycelium and 88 genes upregulated in a mat+ mycelium. Among these genes, we search for the DE genes identified in this study. A total of 32 DE genes (46 %) were identified among the 69 genes which were found to be upregulated in the mat- mycelium (Table S8) and 45 DE genes (51%) are found among the mat+ upregulated genes (Table S9). A statistical analysis indicates that the DE genes are over-represented in both categories, when compared with the 3,466 DE genes in the 10,507 P. anserina’s gene models. This overrepresentation was higher for the mat+ upregulated genes (Chi2, p-value = 0.0003) than for the mat- upregulated genes (Chi2, p-value = 0.0183). This result is likely due to the use of a mat+ maternal strain for the fruiting body development presented here and indicates that half of the genes upregulated in fertilization competent mycelium are also up-regulated in maternal tissues during perithecial development.
Discussion
Mating-type gene expression
In this study, we analyzed on a genome-wide scale the transcriptional changes associated with sexual development in P. anserina over a 10-time points course, from fertilization to ascospore maturation. We found that more than 32% of the P. anserina’s genes are differentially expressed during at least one interval covered by this kinetics. Interestingly, more than 36% of the proteins encoded by these genes do not have a Pfam domain, suggesting that a significant proportion of them remain functionally uncharacterized. Overall, the data generated by our study are rich in novel players involved in fungal sexual reproduction and are likely to provide an exploratory framework to analyze the pathways that are required for the differentiation of the 15 different tissues identified in a perithecium (Bistis et al. 2003).
Among these 15 tissues, the ascogenous hyphae plays a central role in the sexual cycle, since it is where karyogamy, meiosis and ascospore formation take place. Its formation starts at T24 (Figure 1) and requires the expression of the four mating-type genes: FMR1, SMR1, SMR2 and FPR1 (Zickler et al. 1995; Arnaise et al. 2001). The ratio of ascogenous tissue to the other tissues of the perithecium is low, potentially masking the zygotic expression profile. This detection limit of our experiments can be assessed by searching the transcripts of FMR1, SMR1 and SMR2. These three mating-type genes are specific to the mat- nuclei, which come from the male organs and are present only in the ascogonial cell and in the ascogenous hyphae. These three mating-type genes are listed in the DE genes (Table S1), indicating unambiguously that transcripts specific to the ascogenous hyphae are detected in our experiments, in spite of its scarcity. Moreover, upregulations of FMR1 and SMR2 transcripts are detected as soon as T24, correlating with the initial development of ascogenous hyphae (Table S1). In addition, 11 genes upregulated in mat- fertilization competent mycelium are among the DE genes (Table S8), confirming the detection of transcript from mat- male nuclei in the female organ. However, no significant upregulation of the mat- mating-type transcripts is detected in the early steps of perithecial development, i.e. from T0 to the formation of the ascogonial cell at T18 (Figure 1). Correspondingly, no wave I DE genes are identified in genes that accumulate transcripts in the mat- mycelium (Table S8). The earliest DE genes upregulated in mat- competent mycelium belong to wave II (n = 4). Wave IV, which corresponds to upregulation of mat- mating-type genes and their correlated genes (Figure 6), includes 412 genes (Figure 3), among which six belong to genes upregulated in mat- mycelium (Table S8). Overall, these data suggests that critical target genes required for the development from fertilization to the formation of ascogonial cells are not included in the DE genes, because their expression level is below the detection threshold of microarrays. In contrast, the expression of genes required for the subsequent developmental stages (differentiation of ascogenous hyphae, meiosis and ascospore formation) are indeed detected in our experiments.
Surprisingly, FPR1/Pa_1_20590, the mat+ specific mating type gene, does not belong to the DE genes, although this gene is also required for the formation of ascogenous hyphae (Zickler et al. 1995; Arnaise et al. 2001). It is likely that any changes in expression of this gene in zygotic tissues is hidden by its transcription in the maternal tissues, which are predominant in the fruiting body. We have however found 45 DE genes that are specifically upregulated in mat+ mycelium competent for fertilization (Table S9). Among these DE genes, 32 are direct or indirect target genes of FPR1 (Bidard et al. 2011). Their detection proved that FPR1 target genes do not escape our method, although FPR1 itself is not present in DE genes. Moreover, detection starts in the early stages of perithecium development as suggested by FPR1 target genes belonging to wave II (n = 9). Later stages are represented by FPR1 target genes belonging to wave IV (n = 3) and V (n = 2). However, our experiments do not allow us to determine unambiguously whether the differential expression of a gene occurs in the maternal or in the zygotic tissues.
The analysis presented above indicates that our data are pertinent to reveal new pathways in the development of meiotic tissues of P. anserina. Fertilization, which initiates the development of meiotic tissues, is the best-known step of this process. Its molecular description has beneficiated from the yeast model and has revealed that pheromone/pheromone receptor interactions and downstream cascade signal were both conserved in Pezizomycotina and Saccharomycotina (Wilson et al. 2021). Subsequent steps have been cytologically well described in numerous fungi and appear specific to Pezizomycotina. In contrast to fertilization, the molecular mechanisms underlying these various cell programmes are unknown. These blind spots may be exemplified by two steps, which are critical in Pezizomycotina to produce numerous genetically different asci from a single fertilization event : i) several models were proposed for the internuclear recognition following mitotic divisions after fertilization, but no one brings a decisive understanding of this phenomenon (Turgeon and Debuchy 2007), ii) internuclear recognition is followed by the formation of dikaryotic cells and croziers, but the molecular mechanisms controlling their formation are unknown. The similarity of dikaryotic cells and croziers to the basidiomycete clamp cells, which require homeodomain transcription factors for their formation, might suggest similar underlying regulatory circuits, but any functions of homeoboxes in dikaryotic cell and crozier formation were dismissed (Coppin et al. 2012). A third instance of blind spot is the differentiation of croziers into meiocytes. Although the role of peroxisomes in this event is well documented, the signal triggering the change in cellular programme is unknown (Peraza-Reyes and Berteaux-Lecellier 2013). Our analyses revealed co-expression networks, which provide candidates for regulatory circuits of perithecium development and may be involved in the processes presented above. Some particular pathways are discussed below.
Secondary metabolite gene are actively regulated during sexual development.
As in N. crassa (Wang et al. 2022b), C. globosum (Wang et al. 2019) and A. nidulans (Liu et al. 2021), our results further support the active role of SMs in fungal sexual reproduction. The evolutionary conservation of their expression profile, i.e. downregulation from fertilization to meiosis followed by a coordinated upregulation during ascospore formation and maturation suggests that SMs protect fungal progeny, i.e. ascospores, from animal predation. However, we identified several genes of unknown function that were overexpressed during early phases. It is therefore possible that these metabolites protect newly-formed fructification but also that they are recruited to mediate developmental processes. This is reminiscent of the Pa_5_2090/pks-6 gene, which encodes the insecticidal polyketide neurosporin A produced at a similar developmental stage by S. macrospora (Nowrousian 2009) and N. crassa (Zhao et al. 2017) to protect fructifications from arthropod attack. Whether this gene product is also involved in the process of morphogenesis is still an open question (Wang et al. 2022b).
Our data reveals a relationship between sterigmatocystin synthesis and mating-type genes by correlating expression of PaAflR/Pa_2_7360, HMG5/Pa_1_13940 and HMG8/Pa_6_4110 (Table 6 and Figure 6). PaAflR/Pa_2_7360 is a regulator of the sterigmatocystin synthesis cluster (Shen et al. 2019), while HMG5/Pa_1_13940 and HMG8/Pa_6_4110 are HMG-domain TFs acting as hubs of regulation for the mating-type genes and mating-type target genes (Aït Benkhali et al. 2013; Gautier et al. 2018). Our results suggest that PaAflR/Pa_2_7360 controls positively the transcription of the two other genes. Over expression of PaAflR/Pa_2_7360 results in sterility, as does the over-expression of HMG5/Pa_1_13940 (Shen et al. 2019; Aït Benkhali et al. 2013). This result is in agreement with the positive transcriptional control of HMG5/Pa_1_13940 by PaAflR/Pa_2_7360. To date, PaAflR/Pa_2_7360 has not been deleted in P. anserina, and although its orthologue in A. nidulans has been knocked out, the effect on the sexual cycle remains poorly understood (Wang et al. 2022a).
Het genes operates important but yet unknown roles in sexual reproduction
This study unveiled an interesting regulation of the het genes, which are involved in allorecognition and thus control the viability of heterokaryons formed after hyphal fusion during vegetative growth (Daskalov 2023). However, heterokaryotic syncytia are not restricted to vegetative thalli, as fertilization also leads to the formation of such multinucleate cells, and accordingly, several het loci have been shown to control hybrid fertility in P. anserina (Ament-Velásquez et al. 2020). When we compared the het genes expression profiles during sexual reproduction with those assayed during IR (Bidard et al. 2013) and bacterial defense response (Lamacchia et al. 2016) (Table S6 G), we found that they fall into two categories and four classes.
The first category consists of het genes that were always down-regulated during sexual reproduction but were either up-regulated during IR and/or bacterial confrontation (class 1) or down-regulated during IR and up-regulated during bacterial confrontation (class 4). This category therefore contains het genes that play a role only in vegetative growth and/or defense against pathogens, but must otherwise be silenced.
The second category consisted of het genes that were either up-regulated during IR but also at some time points during sexual reproduction (class 2) or were not DE during IR and bacterial confrontation but during sexual reproduction (class 3). Class 2 genes may be recruited to control cell fusion during morphogenesis of fructifications and/or of early zygotic lineage for wave II genes, or of the developing asci and ascospores for waves IV and V genes. Class 3 genes may have been recruited to carry out specific programmed cell death events, not related to hyphal fusion but necessary for sexual morphogenesis. In this context, Pa_3_5550 is of particular interest since it shows a peak of expression during dikaryon formation. Hce2 proteins also belong to this category, suggesting a specific role during both the pre-zygotic (Pa_2_4765) and post-zygotic (Pa_1_17065) stages.
The het-6 and het-13 genes were found to be DE during sexual reproduction in the heterothallic fungus N. crassa, whereas no such observation was made in the homothallic fungus C. globosum (Wang et al. 2019). In P. anserina, both orthologs Pa_6_6160/het-6 and Pa_0_150/het-13 were also DE but showed opposite expression patterns compared to N. crassa. This suggests that the roles of het genes, although central in the sexual developmental network of fungi endowed with two different mating types, are polymorphic.
Rewiring genes in the RNAi pathway?
In P. anserina (Grognet et al. 2019), as in T. reesei (Li et al. 2018) and A. immersus (Malagnac et al. 1997), the rid (RIP defective) ortholog is required for proper sexual reproduction, which is not the case in N. crassa. Pa_1_19440/PaRid showed a short peak of expression at T12 (Figure S5AI), which is consistent with its requirement for the individualization of the dikaryotic cells leading to meiocytes (Grognet et al. 2019). In contrast, N. crassa rid up-regulation started at dikaryon formation and increased continuously throughout the second half of sexual reproduction (Wang et al. 2014). Interestingly, the N. crassa rid expression pattern was similar to that of the P. anserina DNA methyltransferase encoding gene PaDnmt5, which is absent in N. crassa. These results were in agreement with our previous data (Grognet et al. 2019) and strengthened the body of evidence pointing to a role for DNA methylation during sexual development.
Although there is no experimental proof that either somatic or meiotic RNAi pathways (Quelling or MSUD, respectively) works efficiently in P. anserina, some of the corresponding genes are found to be DE during sexual reproduction. Globally the expression patterns of P. anserina genes involved in the RNAi pathways are opposite to those of N. crassa, except for Pa_2_9830/Sad-2 and Pa_5_11380/Sad-5. During MSUD, Sad-2 is thought to interact with the RdRP Sad-1 as a scaffold. Yet, Sad-1 and Sad-2 showed opposite patterns of expression during N. crassa sexual development, which likely discards this hypothesis (Wang et al. 2014). However, Sad-5, which role is yet unknown, has an expression profile that is similar to that of Sad-2. The evolutive conservation of Pa_2_9830/Sad-2 and Pa_5_11380/Sad-5 profiles of expression in both N. crassa and P. anserina, suggests that they have been recruited together to play a role in sexual reproduction. In the homothallic ascomycete S. macrospora, Dcl1, Dcl2, Sms2 and Qde2 are involved in meiocyte formation (Girard et al. 2021), which is not consistent with the expression pattern of the Pa_6_6150/Dcl-1 but rather with that of Pa_5_8910/Sms-2. Altogether, these results suggest a complex role for the RNAi machinery during fungal sexual reproduction.
Conclusion
Our work reveals 3,466 DE genes during the differentiation of fruiting bodies of P. anserina. Only 86 of these DE genes have been previously functionally analyzed (Table S10), leaving a large number of genes for future investigations. Although our analysis does not allow us to determine the DE genes specifically expressed in each tissue of the fruiting body, our method is sensitive enough to detect DE genes in ascogenous hyphae, at a stage where this tissue is represented by few cells in the mass of the fruiting body. The sensitivity of our method is emphasized by the detection of a peak of transcription of Pa_1_19440/PaRid at a time of the early perithecial development corresponding to its expected function. Among the DE genes, we have identified at least two pathways that may control mating-type genes. The first one involves Pa_7_7000, possibly revealing a new regulatory pathway upstream of mating-type genes. The second one is the relationship between PaAflR/Pa_2_7360, HMG5/Pa_1_13940 and HMG8/Pa_6_4110. The latter two TFs control mating-type genes, thus revealing an additional putative pathway regulating indirectly these genes. Our analysis reveals also the intriguing expression patterns of the genes from the MSUD and RNAi pathways, although MSUD has not been found in P. anserina and RNAi not investigated. Further research should be focused in P. anserina on the role of the genes involved in these pathways in other species.
Acknowledgements
We are grateful to S. Saupe for the exchange of data and ideas. Preprint version 1 of this article has been peer-reviewed and recommended by Peer Community In Genomics (https://doi.org/10.24072/pci.genomics.100427; Bloyer 2025).
Funding
This work and the salary of F. Bidard were funded by the French National Research Agency (L’Agence Nationale de la Recherche, ANR) grant number ANR-05-BLAN-0385, project SexDevMycol, coordinator R. Debuchy.
Conflict of interest disclosure
The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.
Data, scripts, code, and supplementary information availability
Data are available online: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE93094. Supplementary data abd figures are available on bioRxiv: https://doi.org/10.1101/2024.12.09.627484 (Bidard et al., 2024). Scripts can be accessed on github (https://github.com/Podospora-anserina/transcriptome_kinetics_development) or via Zenodo: https://doi.org/10.5281/zenodo.15718154 (Lelandais and Grognet 2025).