DOI: 10.5061/dryad.gc2k9r9
title: Data from: Megaphylogeny resolves global patterns of mushroom evolution

description: 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.
description: Results of the BAMM
 analysisWe 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.Varga_2019_NEE_BAMM.zipInput files and the results of the molecular clock analysis (PhyloBayes & FastDate) of the 5,284 taxa data setThere 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.Varga_2019_NEE_PhyloBayes_FastDate.zipMolecular 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.Varga_2019_NEE_Genome_MolecularClock.zipMaximum Likelihood analysis of the 5284 taxa datasetMultiple 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.Varga_2019_NEE_5284taxa_MaximumLikelihood_analysis.zipInput files and the results of the phylogenomic analysisFollowing 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 &gt;=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.Varga_2019_NEE_Genome_Phylogeny.zipSpecies names to track names across analysisA three columned data matrix. Species_names_final column contains names occurring in the article. Alignment_names column contains species names found in in the alignment of 5284 taxa dataset. BayesTrait_names column contains species names used in BayesTraits analysis and BAMM analysis.Species_name_match.csv