Lessons on fruiting body morphogenesis from genomes and transcriptomes of Agaricomycetes

Fruiting bodies (sporocarps, sporophores or basidiomata) of mushroom-forming fungi (Agaricomycetes) are among the most complex structures produced by fungi. Unlike vegetative hyphae, fruiting bodies grow determinately and follow a genetically encoded developmental program that orchestrates their growth, tissue differentiation and sexual sporulation. In spite of more than a century of research, our understanding of the molecular details of fruiting body morphogenesis is still limited and a general synthesis on the genetics of this complex process is lacking. In this paper, we aim at a comprehensive identification of conserved genes related to fruiting body morphogenesis and distil novel functional hypotheses for functionally poorly characterised ones. As a result of this analysis, we report 921 conserved developmentally expressed gene families, only a few dozens of which have previously been reported to be involved in fruiting body development. Based on literature data, conserved expression patterns and functional annotations, we provide hypotheses on the potential role of these gene families in fruiting body development, yielding the most complete description of molecular processes in fruiting body morphogenesis to date. We discuss genes related to the initiation of fruiting, differentiation, growth, cell surface and cell wall, defence, transcriptional regulation as well as signal transduction. Based on these data we derive a general model of fruiting body development, which includes an early, proliferative phase that is mostly concerned with laying out the mushroom body plan (via cell division and differentiation), and a second phase of growth via cell expansion as well as meiotic events and sporulation. Altogether, our discussions cover 1 480 genes of Coprinopsis cinerea, and their orthologs in Agaricus bisporus, chrysosporium, Pleurotus ostreatus, and Schizophyllum commune, providing functional hypotheses for similar to 10 % of genes in the genomes of these species. Although experimental evidence for the role of these genes will need to be established in the future, our data provide a roadmap for guiding functional analyses of fruiting related genes in the Agaricomycetes. We anticipate that the gene compendium presented here, combined with developments in functional genomics approaches will contribute to uncovering the genetic bases of one of the most spectacular multicellular developmental processes in fungi

Bazídiumos gombák termőtestfejlődésében szerepet játszó és biotechnológiai jelentőséggel bíró konzervált, annotálatlan gének azonosítása

In this dissertation, we investigated the initial and final stages of fruiting body development, paying particular attention to conserved, unannotated genes and genes without known function. We investigated genes involved in fruiting body development from two different aspects. First, we examined conserved genes encoding putative cell surface receptors that are overexpressed during fruiting body initiation, and second, we used RNA sequencing to uncover the post-meiotic steps of sporulation. Deletion mutants of two predicted genes encoding proteins with 7-transmembrane domains (p7TMP) showed morphological alterations. Among these mutants, the deletion strain of Δp7TMP79 produced fruiting bodies with previously undescribed morphological characteristics. Differentiation was almost completely absent in the white, round primordia. Due to the striking phenotype of the deletion mutant, we named the deleted gene “snb1”, after “snowball”. The Δsnb1 strain produced fruiting bodies without internal tissue, which allowed us to identify genes involved in tissue differentiation. The fruiting bodies without internal tissues of the Δsnb1 strain allowed the identification of genes involved in tissue differentiation. The phylogenetic analysis revealed that the SNB1 orthogroup represents a previously unknown novel gene family. Previous research on sporulation mainly investigated the steps of meiosis that produce haploid gametes. Transcriptomic analyses not only revealed these processes, but also proved effective in identifying gene candidates whose deletion could impede sporulation, thus serving as valuable resources for the development of industrial spore-free fungal strains. We knocked out three conserved, unannotated, or functionally less known genes in Coprinopsis cinerea, resulting in spore-deficient and completely spore-free strains with no visible morphological defects

Megaphylogeny resolves global patterns of mushroom evolution

Mushroom-forming fungi (Agaricomycetes) have the greatest morphological diversity and complexity of any group of fungi. They have radiated into most niches and fulfil diverse roles in the ecosystem, including wood decomposers, pathogens or mycorrhizal mutualists. Despite the importance of mushroom-forming fungi, large-scale patterns of their evolutionary history are poorly known, in part due to the lack of a comprehensive and dated molecular phylogeny. Here, using multigene and genome-based data, we assemble a 5,284-species phylogenetic tree and infer ages and broad patterns of speciation/extinction and morphological innovation in mushroom-forming fungi. Agaricomycetes started a rapid class-wide radiation in the Jurassic, coinciding with the spread of (sub)tropical coniferous forests and a warming climate. A possible mass extinction, several clade-specific adaptive radiations and morphological diversification of fruiting bodies followed during the Cretaceous and the Paleogene, convergently giving rise to the classic toadstool morphology, with a cap, stalk and gills (pileate-stipitate morphology). This morphology is associated with increased rates of lineage diversification, suggesting it represents a key innovation in the evolution of mushroom-forming fungi. The increase in mushroom diversity started during the Mesozoic-Cenozoic radiation event, an era of humid climate when terrestrial communities dominated by gymnosperms and reptiles were also expanding

Maximum Likelihood analysis of the 5284 taxa dataset

Multiple sequence alignment was carried out for the LSU, ef-1a and RPB2 loci separately using the Probabilistic Alignment Kit (PRANK release 140603). An iterative alignment refinement strategy as described in Tóth et al. was employed: ML gene trees computed from preliminary alignments (using RAxML, see below) were used as guide trees for the next round of multiple alignment for PRANK. After three rounds of iterative refinement, the alignments were further corrected manually using a text editor. Manual curation was restricted to correcting homologous regions erroneously juxtaposed by PRANK. Alignments of individual sequences were concatenated into a superalignment. Maximum likelihood trees for the 5,284-taxon dataset were inferred using the parallel version of RAxML 8.1.2 under the GTR model with gamma distributed rate heterogeneity (4 categories) with three partitions corresponding to the LSU, ef1-α and rpb2 loci. The phylogenomic tree was used as a backbone monophyly constraint. We performed 245 ML inferences and tested whether these trees adequately represented the plausible set of topologies given the alignment. This was done to ensure that phylogenetic uncertainty is properly taken into account in subsequent comparative analyses. If our tree set contains all plausible topologies, then the rolling average of pairwise Robinson-Foulds (RF) distances should show a saturation as a function of increasing the number of trees. To this end, we computed RF distance for each pair of trees for incrementally larger numbers of trees using R package “phangorn” v.2.0.2. We then plotted the rolling average, maximum and minimum values as a function of the number of trees in R. The following files are included: 5284taxa_ML_alignment.fas: input alignment. 5284taxa_ML_alignment_partitions.txt: a file containing the coordinates of the three partitions in the input multiple alignment. 5284taxa_ML_excluded_regions.txt: some of the regions of the alignment were excluded before proceeding to the ML analysis. This file describes the position of the excluded regions. 5284taxa_ML_alignment_exclude.txt: Alignment after excluding regions. 5284taxa_ML_alignment_exclude_partitions.txt: Parition file of the final alignment. 5284taxa_ML_constraint_backbone_tree.tre: A phylogenomic tree which was used to constrain the topology of the backbone of the 5284 taxa phylogeny. 5284taxa_ML_alignment.fas: 245 phylogenetic tree inferred by maximum likelihood analysis

Input files and the results of the molecular clock analysis (PhyloBayes & FastDate) of the 5,284 taxa data set

There are three sub-directories: PhyloBayes_input, PhyloBayes_output and FastDate_analysis. The whole analysis started with PhyloBayes. The PhyloBayes input files are the following. fb_align_*.phy: 10 alignments, containing 543 taxa after randomly deleted the ~90% of the tips from the 5284taxa dataset. fb_calib_*.cal: 10 files, containing species pairs which define the constrains on the MRCA. fb_tree_*.tre: 10 phylogenies, containing 543 taxa after randomly deleted the ~90% of the tips from the 5284taxa dataset. PhyloBayes analyses were run using the 10% subsampled dataset, a birth-death prior on divergence times, an uncorrelated gamma multiplier relaxed clock model and a CAT-poisson substitution model with a gamma distribution on the rate across sites. A uniformly distributed prior was applied to fossil calibration times. All analyses were run until convergence, typically 15,000 cycles. Convergence of chains was assessed by visually inspecting the likelihood values of the trees and the tree height parameter. We sampled every tree from the posterior and after discarding the first 7,000 samples as burn-in we summarized the posterior estimates using the readdiv function of PhyloBayes. The results can be found in the PhyloBayes_output directory. The directory FastDate_anaysis contains the input files. calib_final_tree_*_.cal: 10 files, containing species pairs which define the constrains on the MRCA. FastDate was run on the complete trees (5,284 species) with the node ages constrained to the values of the 95% highest posterior densities of the ages inferred by PhyloBayes. tree_original_*.tree2: 10 phylogenies, containing 5284 taxa. These trees came from the 5284taxa ML analysis. FastDate analyses were run with time discretized into 1,000 intervals and the ratio of sampled extant individuals set to 0.14. The output files are the followings. fastdate_kronogram_*.tree: 10 chronograms inferred by FastDate analysis. transform_to_ultrametric_script.R: An R script which transforms trees to ultrametric. Because rounding issues, a negligible length was added to some of the tips to achieve ultrametric trees. fastdate_kronogram_*.tree2: 10 chronograms used in further analysis

Megaphylogeny resolves global patterns of mushroom evolution

Mushroom-forming fungi (Agaricomycetes) have the greatest morphological diversity and complexity of any group of fungi. They have radiated into most niches and fulfil diverse roles in the ecosystem, including wood decomposers, pathogens or mycorrhizal mutualists. Despite the importance of mushroom-forming fungi, large-scale patterns of their evolutionary history are poorly known, in part due to the lack of a comprehensive and dated molecular phylogeny. Here, using multigene and genome-based data, we assemble a 5,284-species phylogenetic tree and infer ages and broad patterns of speciation/extinction and morphological innovation in mushroom-forming fungi. Agaricomycetes started a rapid class-wide radiation in the Jurassic, coinciding with the spread of (sub)tropical coniferous forests and a warming climate. A possible mass extinction, several clade-specific adaptive radiations and morphological diversification of fruiting bodies followed during the Cretaceous and the Paleogene, convergently giving rise to the classic toadstool morphology, with a cap, stalk and gills (pileate-stipitate morphology). This morphology is associated with increased rates of lineage diversification, suggesting it represents a key innovation in the evolution of mushroom-forming fungi. The increase in mushroom diversity started during the Mesozoic-Cenozoic radiation event, an era of humid climate when terrestrial communities dominated by gymnosperms and reptiles were also expanding

Input files and the results of the phylogenomic analysis

Following the all-versus-all blast using mpiBLAST 1.6.0 with default parameters, we identified gene families using the Markov clustering algorithm MCL 14-137 with an inflation parameter of 2.0. We performed multiple sequence alignment using PRANK 140603. Ambiguously aligned regions were removed from the alignments using Gblocks 091b with the settings -p=yes -b2=26 -b3=10 -b4=5 -b5=h -t=p -e=.gbl. we screened for gene families that contains a single representative gene of each species, or ones that contained inparalogs but no deep paralogs. Deep paralogs were identified following Nagy et al. Gene trees were inferred using the PTHREADS version of RAxML 8.1.2 using the PROTGAMMAWAG model and the standard algorithm. A single inparalog, closest to the root based on root-to-tip patristic distances was retained for each species. Gene families in which >=75% of the species were represented were concatenated into a supermatrix. We used RAxML 8.1.2 to perform ML analysis and bootstrapping on the concatenated dataset, under a WAG model with gamma-distributed rate heterogeneity partitioned by input gene. We ran 100 bootstrap replicates using the rapid hill climbing algorithm. The robustness of the dataset was tested by eliminating incrementally higher numbers of fast-evolving sites using six levels of stringency in Gblocks 091b. Using these parameters, we eliminated 8.5% (-b1=78 -b2=78 -b3=10 -b4=10), 24.3% (-b1=78 -b2=78 -b3=8 -b4=15), 36.7% (-b1=88 -b2=88 -b3=8 -b4=15), 46.4% (-b1=95 -b2=95 -b3=8 -b4=15), 50.1% (-b1=100 -b2=100 -b3=8 -b4=15) and 55.2% (-b2=104 -b3=5 -b4=20) of the least reliably aligned regions of the alignments, resulting in trimmed concatenated datasets with 129.886, 107.496, 89.732, 76.153, 70.862 and 63.309 amino acid sites, respectively. We performed ML phylogenetic inference for each of the reduced datasets in RAxML as described above. The following files were included. RAxML_Genomic_Backbone_Alignment.phy: input alignment for the ML analysis. RAxML_Genomic_Backbone_Phylogeny.tre: ML tree with bootstrap values. MultipleAlignment_GenomeTree_sensitivity_gbl_*_bip_SupplFig2_A-F.phy: six alignment for the sensitivity test. RAxML_GenomeTree_sensitivity_gbl_*_bip_SupplFig2_A-F.tre: six phylogeny, the results of the sensitivity test. These trees were depicted in the Supplement Figure 2 of the article

Results of the BAMM analysis

We used BAMM 2.5.0. (Bayesian Analysis of Macroevolutionary Mixtures), to examine rate heterogeneity across lineages and detect shifts in diversification rates. We analyzed 10 chronograms and ran MCMC analyses for 100 million generations using four independent chains per analysis with 50 million generations as burn-in. Prior parameters were optimized using the setBAMMpriors function in BAMMtools 2.1.6., except for the prior on the expected number of shifts, which was set to 270 based on preliminary runs

Molecular clock analysis of genomic data.

Phylogenomic dataset We used a smaller, more conserved subset of the 568-gene and 104-species phylogenomic dataset (which was computationally not tractable in these analyses). First, we selected the first 70 most conserved genes of the 568-gene dataset by calculating the mean genetic distances for each gene using the dist.alignment function of the seqinR R package v.3.4-5. To enable a more accurate placement of fossil calibration points we added additional three species (Cyathus striatus, Pycnoporus cinnabarinus, Suillus brevipes) to this dataset (Supplementary Table 1) and excluded two taxa that harbored ambiguous positions. We searched homologous sequences in the additional genomes using blastp v.2.7.1 with one randomly selected gene from each of the 70 gene families as query. We selected the best hit (smallest E-value) as a 1-to-1 ortholog if the second best hit had a significantly worse E-value (by 20 orders of magnitude). Protein clusters were aligned by PRANK v.100802 using default settings. Next, conserved blocks of the alignments were selected using Gblocks V.0.91b with default settings except for the minimum length of a block which was set to 5 and gap positions in half of the sequences were allowed. A phylogenomic tree was constructed by RAxML v.8.2.11 under WAG+G substitution model partitioned by gene. Calibrations To dissect sources of differences in molecular age estimates, we ran analyses under 3 fossil calibration schemes (Supplementary Data 2) and the 105-species phylogenomic tree. First we used the same fossil calibration scheme as for the 5,284-species phylogenetic dataset (“Default calibration scheme”). Next, we replicated the analyses of Kohler et al. on our tree, using the fossil calibration points from Kohler et al. (“Kohler et al. calibration scheme 1”). For this, we placed the suilloid ectomycorrhiza fossil in the split of Suillinae/Paxillinae/Sclerodermatinae and Archaeomarasmius leggettii in the mrca of Gymnopus luxurians and Schizophyllum commune with uniformly distributed 40–60 mya and 70–110 mya time priors, respectively. Finally, we used the calibrations used by Kohler et al. (2) but placed the two fossils in the mrca-s of the Suillaceae and marasmioid clade, respectively (Kohler et al. calibration scheme 2”). In all analyses we constrained the age of the root to be between 300 mya and 600 mya. Penalized likelihood analysis in r8s We ran a series of molecular clock analyses in r8s v.1.81. A cross-validation analysis was performed to determine the optimal smoothing parameter (λ) by testing values across 7 orders of magnitude starting from 10-3. The additive penalty function was applied and the optimization was run 25 times starting from independent starting points. In one optimization step, after reaching an initial solution, the solution was perturbed and the truncated Newton (TN) optimization was rerun 20 times. We compared the results of previous studies to that of analyses across seven ancestral nodes in Agaricomycotina (Supplementary Data 2). Bayesian molecular clock dating We used the mcmctree method implemented in PAML version 4.8a. The independent-rates clock model, a WAG substitution model and approximate likelihood calculation were used. The birth rate, the death rate and the sampling fraction of the birth-death process were set to 1, 1 and 0.14 respectively. The shape and the concentration parameter of the gamma-Dirichlet prior for the drift rate coefficient (σ2) was set to 1 and three different scale parameters were tested (10, 100, 1,000) to see their effect on the time estimates. The substitution rates of each gene were estimated by codeml under a global clock model, to set the parameters of the gamma-Dirichlet prior for the overall rate. By calculating the mean substitution rate of the loci and examining the density plot of the rates we set up a prior which reasonably fitted the data: the shape parameter, the scale parameter and the concentration parameter were set to 5, 90.7441 and 1, respectively, resulting in an average substitution rate per site per time unit of 0.055. We set the time unit to 100 myr and applied uniform priors on 8 fossil calibrations with lower and upper hard bounds. MCMC (Markov chain Monte Carlo) analysis was run for 80,000 iterations, discarding the first 20,000 iterations as a burn-in and sampling every 30th tree from the posterior. After three independent analyses were run the convergence of log-likelihood values was visually inspected and the estimated ages were compared between replicates. The following files are included.Genome_MolClock_ML_input_alignment.fasta: Maximum Likelihood input alignment. Genome_MolClock_ML_input_alignment_partitions.txt: Maximum Likelihood input alignment partitions. Genome_MolClock_ML_output_tree.tree: Genome MolClock ML output tree. Genome_MolClock_mcmctree_input_alignment.phy: mcmctree input alignment. Genome_MolClock_mcmctree_input_tree.tree: mcmctree input tree. Genome_MolClock_mcmctree_Hessian_matrix.inBV: mcmctree input Hessian matrix. Genome_MolClock_mcmctree_output.out: Genome MolClock mcmctree output. r8s_default_calibration_scheme.tre: R8S output1. r8s_kohler_etal_calibration_scheme1.tre: Genome MolClock R8S output2. r8s_kohler_etal_calibration_scheme2.tre: R8S output3

Data from: Megaphylogeny resolves global patterns of mushroom evolution

Mushroom-forming fungi (Agaricomycetes) have the greatest morphological diversity and complexity of any group of fungi. They have radiated into most niches and fulfill diverse roles in the ecosystem, including wood decomposers, pathogens or mycorrhizal mutualists. Despite the importance of mushroom-forming fungi, large-scale patterns of their evolutionary history are poorly known, in part due to the lack of a comprehensive and dated molecular phylogeny. Here, using multigene and genome-based data, we assemble a 5,284-species phylogenetic tree and infer molecular clock ages and broad patterns of speciation/extinction and morphological innovation in mushroom-forming fungi. Agaricomycetes started a rapid class-wide radiation in the Jurassic, coinciding with the spread of (sub)tropical coniferous forests and a warming climate. A possible mass extinction, several clade-specific adaptive radiations, and morphological diversification of fruiting bodies followed during the Cretaceous and the Paleogene, convergently giving rise to the classic toadstool morphology, with a cap, stalk, and gills. This morphology is associated with increased rates of lineage diversification, suggesting it represents a key innovation in the evolution of mushroom-forming fungi. The accumulation of mushroom diversity started during the Mesozoic-Cenozoic radiation event, an era of humid climate when terrestrial communities dominated by gymnosperms and reptiles were also expanding