Genome expansion and lineage-specific genetic innovations in the forest pathogenic fungi Armillaria

Armillaria species are both devastating forest pathogens and some of the largest terrestrial organisms on Earth. They forage for hosts and achieve immense colony sizes via rhizomorphs, root-like multicellular structures of clonal dispersal. Here, we sequenced and analysed the genomes of four Armillaria species and performed RNA sequencing and quantitative proteomic analysis on the invasive and reproductive developmental stages of A. ostoyae. Comparison with 22 related fungi revealed a significant genome expansion in Armillaria, affecting several pathogenicity-related genes, lignocellulose-degrading enzymes and lineage-specific genes expressed during rhizomorph development. Rhizomorphs express an evolutionarily young transcriptome that shares features with the transcriptomes of both fruiting bodies and vegetative mycelia. Several genes show concomitant upregulation in rhizomorphs and fruiting bodies and share cis-regulatory signatures in their promoters, providing genetic and regulatory insights into complex multicellularity in fungi. Our results suggest that the evolution of the unique dispersal and pathogenicity mechanisms of Armillaria might have drawn upon ancestral genetic toolkits for wood-decay, morphogenesis and complex multicellularity. Fungi of the genus Armillaria include devastating forest pathogens that cause root rot disease in many plants. Sequencing genomes and transcriptomes of several species, the authors reveal the genetic basis of dispersal, multicellular development and pathogenic mechanisms in Armillaria.

T he genus Armillaria causes root rot disease in both gymnoand angiosperms, in forests, parks, and even vineyards in more than 500 host plant species 1 across the world. Most Armillaria species are facultative necrotrophs, which, after colonizing and killing the root cambium, transition to a saprobic phase, decomposing dead woody tissues of the host. As saprotrophs, Armillaria spp. are white rot (WR) fungi, which can efficiently decompose all components of plant cell walls, including lignin, (hemi-)cellulose and pectin 2 . They produce fleshy fruiting bodies (honey mushrooms) that appear in large clumps around infected plants and produce sexual spores. The vegetative phase of Armillaria is predominantly diploid rather than dikaryotic like most basidiomycetes.
Individuals of Armillaria can reach immense sizes and include the 'humongous fungus' , one of the largest terrestrial organisms on Earth 3 , measuring up to 965 hectares and 600 tons 4 , and can display a mutation rate ≅ 3 orders of magnitude lower than most filamentous fungi 5 . Individuals reach this immense size via growing rhizomorphs, dark mycelial strings 1-4 mm wide that allow the fungus to bridge gaps between food sources or host plants 1,6 (hence the name shoestring root rot). Rhizomorphs develop through the aggregation and coordinated parallel growth of hyphae, similar to some fruiting body tissues 7,8 . As migratory and exploratory organs, rhizomorphs can grow approximately 1 m yr −1 and cross several metres underground in search for new hosts, although roles in uptake and longrange translocation of nutrients have also been proposed 1,9,10 . Root contact by rhizomorphs is the main mode of infection by the fungus, which makes the prevention of recurrent infection in Armillariacontaminated areas particularly difficult 1 . Despite their huge impact on forestry, horticulture and agriculture, the genetics of the pathogenicity of Armillaria species is poorly understood. The only -omics data published so far have highlighted a substantial repertoire of plant cell wall degrading enzymes (PCWDE) and secreted proteins, among others, in A. mellea and A. solidipes 11,12 , while analyses of the genomes of other pathogenic basidiomycetes (such as Moniliophthora 13,14 , Heterobasidion 15 and Rhizoctonia 16 ) identified genes coding for PCWDEs, secreted and effector proteins or secondary metabolism (SM) as putative pathogenicity factors. However, the lifecycle and unique dispersal strategy of Armillaria prefigure other evolutionary routes to pathogenicity, which, along with other potential genomic factors (such as transposable elements 17 ) are not yet known.
Here, we investigate genome evolution and the origin of pathogenicity in Armillaria using comparative genomics, transcriptomics

Results
We report the genomes of A. ostoyae, A. cepistipes, A. gallica and A. solidipes sequenced using a combination of PacBio and Illumina technologies ( Table 1). Genomes of Armillaria species were assembled to 103-319 scaffolds comprising 58-85 Mb and were predicted to contain 20,811-25,704 genes. In comparison, other sequenced species of the Physalacriaceae, Flammulina velutipes and Cylindrobasidium torrendii have 12,218 (35.6 Mb) and 13,940 (31.5 Mb) genes, respectively, while the sister genus of Armillaria, Guyanagaster necrorhiza, has 14,276 (53.6 Mb) (Fig. 1). Armillaria species share significant synteny, comprising macro-to microsynteny ( Supplementary Fig. 1), whereas mesosynteny, which is characteristic of certain fungal groups 18 , was not observed. The transposable element (TE) content of Armillaria genomes shows a modest expansion relative to other Agaricales and an even distribution along the scaffolds, suggesting that their genome expansion is not driven by transposon proliferation, as observed in other plant pathogens 17  Phylogenomic analysis based on 835 conserved single copy genes (188,895 amino acid sites) confirmed the position of Armillaria in the Physalacriaceae, with Guyanagaster and Cylindrobasidium as their closest relatives (Fig. 1a). We estimated the age of pathogenic Armillaria spp. at 21 million years (Myr) and their divergence from Guyanagaster at 42 Myr ( Supplementary Fig. 4), coincident with decreasing temperatures and the spread of deciduous forests in the Eocene. Reconstruction of genome-wide gene duplication and loss histories in 27 Agaricales species revealed an early origin for most genes, followed by lineage-specific gene losses in most family and genus level groups, except Armillaria, which showed a net genome expansion: 15,787 protein-coding genes were inferred for the most recent common ancestor (MRCA) of Armillaria (2,012 duplications, 945 losses), as opposed to 14,720 and 14,687 for the MRCA of Armillaria and Guyanagaster and that of Armillaria, Guyanagaster and Cylindrobasidium, respectively (Fig. 1a,d, Supplementary Fig. 5). Further expansion to 19,272 genes was inferred for the MRCA of A. solidipes, A. ostoyae, A. epistipes and A. gallica (3,192 duplications, 607 losses), although the highly fragmented A. mellea assembly might cause some duplications to map to this instead of the preceding node.
Duplicated genes were enriched in functions related to chitin and cellulose binding, polysaccharide metabolism (peroxidase, lyase, hydrolase and oxidoreductase activity), peptidase activity, transmembrane transport, extracellular region and gene expression regulation (p < 0.05, according to a hypergeometric test, HT, Supplementary Table 2). Two hundred and fourteen domains were significantly overrepresented in Armillaria genomes (p < 0.05, HT) relative to other agarics, including several peptidase, glycoside hydrolase (GH) and pectinase domains, LysM (CBM50-s) domains, expansins, multicopper oxidases as well as several pfams related to secondary metabolism (Supplementary Table 3). We found 570 protein clusters specific to Armillaria and Guyanagaster or subclades therein (Supplementary  Table 4). Although most of these (70%) had no functional annotations, they included CE4 chitooligosaccharide deacetylases, CBM50-s, iron permeases (FTR1), and 19 transcription factor (TF) families, among others. Taken together, these results suggest that gene family expansion was the predominant mode of genome evolution in Armillaria and that the genome expansion is largely concerned with diverse extracellular functions, including several lineage-specific innovations, some of which had previously been associated with pathogenicity.
Putative pathogenicity-related genes in Armillaria. Plant pathogenic fungi possess diverse gene repertoires for invading host plants and modulating their immune systems [18][19][20] . We catalogued 20 families of putative pathogenesis-related genes to assess whether Armillaria shares expansions of these families with other plant pathogens (Supplementary Table 5). Armillaria species are enriched in expansins (p = 4 × 10 −5 , Fisher exact test, FET) and possess many cerato-platanin genes, which contribute to unlocking cell-wall polysaccharide complexes and cause cell death in host plants, respectively, and might act as first-line cell-lysis weaponry during invasion (Fig. 1a). Carboxylesterases and salicylate hydroxylases show moderate enrichment in both Armillaria and Moniliophthora species, along with other weakly pathogenic taxa (such as Marasmius fiardii). Salicylate hydroxylases have been implicated in developing a tolerance to salicylic acid, which is released to block the jasmonic acid defense pathway of plants on infection. In contrast, cutinases (CE5) are missing from Armillaria species but are present in Moniliophthora spp. 14 , which primarily infect through leaves. CBM50 domains are overrepresented in Armillaria compared to other Agaricales (p = 2 × 10 −8 , FET), while GH75 chitosanases are exclusively found in Armillaria and Moniliophthora species (Fig. 1a, Supplementary Table 5). In other plant pathogens, these are involved in modifying fungal cell wall composition and/or capturing chitin residues to mask chitin-triggered immune signals and evade detection by the host plants [21][22][23][24][25] , suggesting similar roles in Armillaria.
Homologs of SM genes reported from Heterobasidion 15 are overrepresented in Armillaria species (p = 0.03 × 10 −13 , FET, Fig. 1a, Supplementary Tables 3 and 5). These include terpene cyclases, non-ribosomal peptide synthase-like prenyl transferases, and halogenases, as well as trichodiene and polyprenyl synthases, which have been linked to fungal pathogenesis, virulence and competition with other microbes 26 . Small secreted proteins (SSPs, < 300 amino acids with a secretory signal) can act as effectors in mutualistic and pathogenic interactions in various fungal groups 19,27 . On average, we found more SSPs (n = 669) in Armillaria species (Fig. 1a) than in saprotrophs (n = 552) or ECM fungi (n = 563), although this is consistent with the larger genomes of Armillaria species (p> 0.05, FET). In contrast, Major Facilitator Superfamily 1 (MFS1) and cytochrome p450 families are expanded in Armillaria species (p = 4 × 10 −10 and p = 3 × 10 −30 ). Both superfamilies have been associated with fungal    30 and might indicate links to dicot pathogenicity 20 . The PCWDE repertoire of Armillaria species underpins their ability to act as powerful WR decayers and provides a pool to draw upon as necrotrophic pathogens. It might enable them to gain access to wood and avoid competition with other microbes by damaging live trees, a strategy that is unavailable to most WR fungi.
Expression profiles of rhizomorphs. Armillaria species spread either clonally by rhizomorphs or with sexual spores produced on fruiting bodies. Rhizomorphs are unique multicellular structures of Armillaria species, but their functions and morphogenetic origins are debated 10 . Differential expression analysis of actively growing rhizomorph tips (0-5 cm from the apex) identified 1,303 and 1,610 genes over-and underexpressed relative to vegetative mycelium grown on the same medium (FC VM > 4, p ≤ 0.05, Supplementary Table 6), respectively, marking one of the largest expression changes in our experiments (Fig. 2). Similarly, the highest number of unique proteins (n = 729) was detected in rhizomorphs compared to vegetative mycelia (Fig. 2c [8][9][10][11][12], whereas that of many putative morphogenesis-related genes was shared with fruiting bodies. Rhizomorphs express an evolutionarily young transcriptome, with most upregulated genes specific to A. ostoyae or the Armillaria clade (including Guyanagaster Fig. 2b). Most of these young genes lack functional annotations but were conserved in Armillaria species, suggesting important functions in the genus. Genes belonging to older phylostrata had comparatively more functional annotation terms. We found 414 genes that had expression maxima in rhizomorphs and were significantly upregulated relative to vegetative mycelium (fold change: FC VM ≥ 4, Supplementary Table 7). The most highly upregulated genes included expansins (log 2 FC VM = 10.44), bZip and C 2 H 2 transcription factors, three caspase genes, hydrophobins, cytochrome p450s, GHs as well as several unannotated genes. We found an overexpression of 5-oxoprolinases downstream factors of neutralizing intracellular H 2 O 2 and signs of intensified biogenesis and cargo of extracellular proteins (Supplementary Figs 13,14). A moderate set of PCWDE genes was expressed in rhizomorphs, suggestive of assimilative or invasive properties, although the highest PCWDE suite was observed in vegetative mycelium .
We observed significant overexpression of three cerato-platanins, three expansins, two carboxylesterases, as well as SM-related trichodiene (five genes), polyketide (six genes) synthases and a polyprenyl synthase, among others (Supplementary Table 8). Notably, all expansins showed upregulation in fruiting bodies too, which could indicate a role in fungal cell wall remodelling instead of cellulose degradation, to which fungal expansin-related genes have mostly been linked 31 . A gene (ARMOST_01259) coding for a family 6 bacterial extracellular solute-binding protein (SBP) showed an expression peak in rhizomorphs, but low or negligible expression in other stages. Bacterial SBPs associate with ABC transporters to function in Fe 3+ transport and are required for the pathogen's survival in the host. This gene is a member of an Armillaria-specific cluster with homologs in all Armillaria species and bacterial proteins as closest BLAST hits (although maximum likelihood (ML) gene trees don't support horizontal gene transfer). Its Fe 3+ transporting ability and nearly exclusive expression in rhizomorphs suggests a role in substrate and/or host exploitation.

Morphogenesis.
The morphogenetic machinery underlying rhizomorph development is among the least known aspects of the biology of Armillaria. As multicellular structures, rhizomorphs express a variety of genes encoding cell-wall proteins, including hydrophobins, pore-forming toxins, two CBM67 and four ricin-B-lectins, an annexin and a cell-wall integrity sensor, among others, indicating several fruiting body-like functions, such as hyphal adhesion, communication or defense (Supplementary Table 9). We found cell-wall biosynthesis genes to be generally but moderately upregulated in rhizomorphs (and stipes) ( Supplementary Fig. 15). Further, several GMC oxidoreductases, two mating-type pheromone receptors, CBM50-s, a CBM5_12 and chitooligosaccharide deacetylase genes were significantly overexpressed (the latter two also in stipes). Four hydrophobins reached their highest expression in rhizomorphs whereas two showed high expression in rhizomorphs and stipes, but not in caps or vegetative mycelium (Supplementary Table 9, Supplementary Fig. 16). A homolog of the Cryptococcus red and far-red sensing Tco3 photoreceptor was expressed at high levels (log 2 FC = 4.36). Although its function in Cryptococcus is unknown 32 , its overexpression in rhizomorphs and in brown-film forming mycelia of Lentinula 33 could suggest morphogenesisrelated functions.
We detected 19 significantly upregulated TFs, ten of which had peak expression in rhizomorphs across our experimental conditions. Although based on global TF expression, rhizomorphs are most similar to vegetative mycelium ( Supplementary Fig. 17), and several TFs showed shared overexpression in rhizomorphs and fruiting bodies. For example, the expression of ARMOST_01275 peaked in rhizomorphs and stipes, whereas a zf-Mynd TF was highly upregulated in rhizomorphs and all fruiting body stages, and is thus a candidate for governing complex multicellular development in A. ostoyae.
Fruiting bodies. The production of fruiting bodies is probably the largest morphogenetic transition in the fungal lifecycle. It involves the reorganization of hyphal growth patterns and the execution of a complex developmental program. Indeed, we detected the largest number of differentially expressed genes (DEGs) and the second largest change in proteome in stage I primordia (P1, Fig. 2, Supplementary Table 6). Upregulated genes were enriched for gene expression regulation, lipid metabolism, amino acid transport and ribonuclease activity (Fig. 2a, Supplementary Table 10) and included ten TFs and two Dicer-like genes. Following stage I primordia, we tracked expression levels in two cell lineages. While the stipe lineage showed minor changes throughout development (< 150 DEGs, < 300 unique proteins), cap differentiation included up to 1,037 DEGs and 646 unique proteins related to signal transduction, carbohydrate metabolism or the regulation of biological processes (Supplementary Table 6). Genes upregulated in gills (n = 502, 56 unique proteins) were enriched for functions related to protein phosphorylation, carbohydrate metabolism, protein kinase activity and ribonucleotide binding.
Hydrophobin genes showed stage and tissue-specific expression patterns, but were generally not expressed in caps ( Supplementary  Fig. 16), suggesting alternative sources of hydrophobin-related functionalities there. An analysis of cysteine-rich, hydrophobic proteins (excluding hydrophobins) revealed 33 and 13 developmentally regulated genes with expression maxima in fruiting bodies and caps, respectively (Supplementary Table 11). In addition, three cell wall galactomannoproteins show high expression in

NATuRE ECOlOgy & EvOluTiON
fruiting bodies (up to logFC VM = 13.5, Supplementary Table 12), two of which were expressed only in caps, whereas ARMOST_19505 was highly expressed throughout development. These genes are homologous to the Aspergillus hydrophobic cell-surface protein HsbA, which has been hypothesized to recruit hydrolytic enzymes during lignocellulose degradation 34 and appressorium attachment to host surfaces 35 . As lignocellulose-related processes are inactive in fruiting bodies, their expression suggests specific roles during cap development, possibly in cell-wall remodelling. Several chitin metabolism-related genes were upregulated in fruiting bodies, including GH88 chitinases, a CBM50-containing gene, a GH75 and six chitooligosaccharide deacetylases (Supplementary  Table 12), which might be related to generating developmentspecific cell-wall architectures. Interestingly, certain HTPs, GMC oxidoreductases (AA3) and pectinolytic genes, generally linked to lignocellulose degradation, were also upregulated in fruiting bodies ( Supplementary Fig. 18). Several defense-related genes (poreforming toxins, lectins and so on) 36,37 also showed significant upregulation through development (Supplementary Figure 19,  Supplementary Table 13), indicating a phylogenetically and functionally diverse defense arsenal expressed in both fruiting bodies and rhizomorphs.

Shared morphogenetic machineries between rhizomorphs and fruiting bodies.
Rhizomorphs share a complex multicellular organization with fruiting bodies. Consistent with this, we observed several genes with similar expression patterns in rhizomorphs and various fruiting body tissues (Fig. 3). These included two matingtype pheromone receptors and the white collar 1-2 genes, which mediate the initiation of fruiting body development and could point to shared developmental origins of rhizomorphs and fruiting bodies. There were 442 genes that showed greater than fourfold elevated expression over vegetative mycelium and relatively constant expression across rhizomorph and fruiting body samples (Fig. 3a, Supplementary Table 14), suggesting they are linked to complex multicellularity in A. ostoyae. A systematic analysis identified 2,225 genes that showed higher expression in rhizomorphs and stipes than in corresponding cap tissues and vegetative mycelium, of which 63 were at least fourfold more abundant in all four pairs of successive developmental stages (Fig. 3b). These included most of the yeast cell-wall biosynthesis pathway homologs (Supplementary  Table 15, Supplementary Fig. 20), pore-forming toxins, hydrophobins, TFs and SM genes, among others. For caps and rhizomorphs 1,728 and 28 such genes were found, respectively, including several GMC oxidoreductases, a GH88, a ricin-B-lectin and other SM genes

NATuRE ECOlOgy & EvOluTiON
( Fig. 3c, Supplementary Table 15). This indicates that rhizomorph development draws extensively on fruiting body genes, and makes us speculate that fruiting body development could have been the cradle for the evolution of rhizomorphs in Armillaria.
To assess whether co-regulated genes possess a common r-regulatory signature, we searched for putative TF binding sites (TFBS) by de novo motif discovery 38 in the promoters of 63 stipe/rhizomorph and 28 cap/rhizomorph co-expressed genes as well as of those that showed constant high expression in all complex multicellular structures. The predicted motifs and their UPGMA trees is shown on Fig. 3. All predicted motifs matched experimentally determined yeast TFBS sequences in the JASPAR database 39 , which suggests strong biological relevance to our predictions, despite the large phylogenetic distance between S. cerevisiae and A. ostoyae. Taken together, the detected co-expressed genes and the presence of shared motifs in their promoters could suggest common regulatory mechanisms for rhizomorphs and fruiting bodies and probably represent elements of the regulatory networks that govern multicellular development as well as cap and stipe differentiation in Armillaria.

Discussion
Forest pathogens in the genus Armillaria evolved from saprotrophic ancestors in the Agaricales. They have unusually large genomes for WR saprotrophs, which evolved mostly by gene family diversification, in contrast with genome expansions in other plant pathogenic fungi, which are primarily driven by TE proliferation 17 . We found many lineage-specific genes (including putative pathogenicity factors and pectinolytic families) and coincident overrepresentation of several pathogenicity-related genes with other Basidiomycota (particularly in Moniliophthora species), suggesting convergent origins of pathogenicity in the Agaricales. Armillaria species encode a full complement of PCWDE genes, which comes as no surprise given their ability to cause WR, but provides a resource to draw on in the evolution of pathogenicity. This could give Armillaria spp. early access to dead wood and a strategy to bypass competition with other microbes.
Rhizomorphs are some of the most unique structures of Armillaria spp. that enable them to become the largest terrestrial organism on Earth 3,5 . They express a wide array of genes involved in secondary metabolism, defense, PCW degradation and to a lesser extent pathogenesis, which indicate active nutrient uptake and adaptations to a soil-borne lifestyle where competition with other microbes and defense against predators are crucial. This underpins their role as exploratory and assimilating organs adapted to bridge the gaps between food sources or potential host plants. In terms of morphogenesis, rhizomorphs resemble fruiting bodies, with many stipe-and to a smaller extent cap-upregulated genes showing concomitant upregulation and shared cis-regulatory signatures. Both rhizomorphs and fruiting bodies show key traits of complex multicellularity, such as three-dimensional organization, cell adhesion or a highly integrated developmental program. We hypothesize that the observed similarity in gene expression indicates common developmental origins and suggest that the evolution of rhizomorphs may have extensively drawn upon the genetic toolkit of fruiting body development in the Agaricomycotina. We identified genes putatively involved in cap and stipe morphogenesis, as well as some co-expressed in all complex multicellular stages and candidate TF-binding motifs within their promoters. These represent putative building blocks of the gene regulatory circuits that govern mushroom development and enabled us to zero in on the genetic bases of complex multicellularity in fungi. This study has provided comparative genomics insights into the evolution, pathogenicity and multicellular development of a group of devastating forest pathogens. It should facilitate further understanding of the biology of Armillaria, which, combined with new genomic resources and in planta interrogation of its pathogenic behaviour, could soon bring the development of efficient strategies for containing the spread and damage of Armillaria root-rot disease in various forest stands within reach.

Strains and fungal material used for genome sequencing. The diploid A. ostoyae
C18 is a field isolate from Switzerland 40 and the haploid C18/9 is derived from C18 as a single spore isolate. The A. cepistipes B5 haploid 41 was originally isolated from a fruiting body on Fagus sylvatica in Italy. The A. ostoyae C18/9 and A. cepistipes B5 haploid isolates were cultured on Roth and Shaw plates (for 1 l RS: 40 g malt extract, 20 g dextrose, 5 g bacto peptone, 19 g agar) covered with cellophane sheets, and incubated at 25 °C for three weeks. Before genomic DNA extraction, fungal mycelia were detached and harvested from the cellophane sheets, and frozen in liquid nitrogen. The A. gallica 21-2 and A. solidipes 28-4 strains were maintained on 2% malt extract, 2% agar medium. To produce mycelium for DNA extraction, strains were grown in liquid CYM medium (d-glucose 20 g, yeast extract 2 g, peptone 2 g, MgSO 4 .7H 2 O 0.5 g, KH 2 PO 4 0.46 g, K 2 HPO 4 1 g, 1 l water) in Petri dishes with 8-10 small pieces of inoculum floated on the air-water interface. Genomic DNA was extracted using the PowerMax MOBIO DNA isolation kit according to the manufacturer's instructions. Note that a previous study 42 proposed that the A. solidipes and A. ostoyae are one species and that the earlier name A. solidipes from a North American collection should always be used in place of A. ostoyae. For consistency with recent historical usage, however, we choose to consider A. ostoyae in Europe as a separate, but closely related species to A. solidipes in North America We developed an in vitro fruiting system for A. ostoyae C18 based on the protocol published for A. mellea 43 . A modified RST medium, termed RSTO, was prepared in 720 ml jars and inoculated with A. ostoyae C18, incubated for 28 days at 24 °C in the dark, then placed in a growth chamber at 15 °C with a 10/14 hour light/dark cycle (light intensity: 11 µ E m −2 s −1 ). Vegetative mycelium, rhizomorphs, stage 1 and 2 primordia, young and mature fruiting bodies were harvested as shown on Fig. 1. Fruiting body stages were defined to conform to the general notation 44 of mushroom developmental stages as closely as possible. Stage 1 primordia were defined as fruiting body initials up to 1-2 mm tall without a clear differentiation of a pileal section. Stage 2 primordia had a developed buttonshaped cap, up to 7-10 mm tall. Young fruiting bodies were 30-50 mm tall with a developed hymenial cavity but prior to cap expansion. From stage 2 primordia, caps and stipes were separated and RNA was extracted separately. Mature fruiting bodies were separated into stipe, gills and cap. RNA was extracted by using the RNEasy Midi kit (QIAGEN) following the manufacturer's instructions. RNA quality was checked by gel electrophoresis and the Agilent 2200 TapeStation before library preparation. Biological triplicates were analysed for all sample types.
For RNA sequencing (RNA-seq), A. cepistipes was grown on MEA (for 1 l MEA: 20 g malt extract, 0.5 g yeast extract, 15 g agar), RST (15 g Picea abies sawdust, 30 g rice, approximately 100 ml water mixed together in 720 ml jars, sterilized and overlaid with a layer of homogenized tomato approximately 2 cm thick followed by another round of sterilization), RSTO (RST medium with additional 100 g minced orange), Orange media (3 roughly chopped oranges in 720 ml jars) and RNA was extracted using the RNEasy Midi kit (QIAGEN).
Genome sequencing and assembly. The haploid genomes of A. ostoyae and A. cepistipes were sequenced employing the PacBio (Pacific Biosystem) RS II platform (Functional Genomics Center; http://www.fgcz.ch). For preparing the sequencing libraries, 10 μ g of gDNA aliquots were mechanically sheared to an average size of 10 kb using the Covaris gTube (KBiosciences p/n 520079) in an Eppendorf microcentrifuge, and the fragment size distributions were assessed applying the Bioanalyzer 2100 12K DNA-Chip assay (Agilent p/n 5067-1508). Five micrograms of the sheared gDNA aliquots were DNA damage repaired, end-repaired and the final SMRTbell templates were created by blunt end ligation and exonuclease treatments. The libraries were quality inspected on an Agilent Bioanalyzer 12K DNA Chip and quantified on a Qubit.1 Fluorimeter (Life Technologies). The SMRTbell was set up by using the DNA Template Prep Kit 2.0 (3 kb to < 10 kb) (Pacific Biosciences p/n 001-540-835). A ready to sequence SMRTbell-Polymerase Complex was arranged by applying the P4 DNA/ Polymerase binding kit 2.0 (Pacific Biosciences p/n 100-236-500) according to the manufacturer instructions. Libraries were sequenced on 15 SMRT cell v3.0 (Pacific Biosciences p/n100-171-800), taking 1 movie of 120 minutes each per SMRT cell. The MagBead loading (PacBio p/n 100-133-600) technique served to improve the enrichment for the longer fragments. Final sequencing reports were generated for every cell, via the SMRT portal, to assess the adapter dimer contamination, the sample loading efficiency, the obtained average read length and the number of filtered sub-reads.
The HGAP3 45 workflow of the SMRT Analysis suite v2.3 was used to create an initial assembly. After removal of redundant contigs, a scaffolding using PBJelly2 46 and FinisherSC 47 followed by polishing via applying the RS Resequencing protocol (SMRT Analysis suite) was performed in four iterations. The final scaffold set was checked for miss-assemblies using the RS BridgeMapper protocol (SMRT Analysis

NATuRE ECOlOgy & EvOluTiON
suite) and corrected if necessary. The mitochondrial scaffolds were first identified using a BLASTn search and then circularized by merging and truncating the overlapping ends. To correct PacBio reads, Illumina sequencing was carried out by shotgun sequencing of a 350-450 bp library with paired-end 100 bp reads to 180fold coverage using HiSeq 2000 (Functional Genomic Center). Finally, Pilon 48 was applied to further polish the scaffolds with the genomic Illumina reads.
Draft gene models for A. ostoyae and A. cepistipes were generated by three de novo prediction programs: (1) Fgenesh 49 with different matrices (trained on Aspergillus nidulans, Neurospora crassa and a mixed matrix based on different species); (2) GeneMark-ES 50 ; and (3) Augustus 51 with RNA-seq-based transcripts as training sets. Annotation was aided by exonerate 52 hits of protein sequences from Armillaria species to uncover gene annotation gaps and to validate de novo predictions. Transcripts were assembled on the RNA-seq data sets using Trinity 53 . The different gene structures and evidences (exonerate mapping, RNA-seq reads and transcripts) were visualized in GBrowse 54 allowing manual validation of coding sequences with a focus on chitin, cellulose, pectin, lignin, SM key genes and other genes of interest. The best fitting model per locus was selected manually and gene structures were adjusted by splitting or fusion of gene models and redefining exon-intron boundaries if necessary; tRNAs were predicted using tRNAscan-SE 55 . The predicted protein sets were searched for highly conserved single-(low) copy genes to assess the completeness of the genomic sequences and gene predictions. Orthologous genes to all 246 single-copy genes were searched for all proteomes by blastp comparisons (eVal: 10-3) against the single-copy families from all 21 species available from the FunyBASE 56 . Additionally, the proteomes were searched for the 248 core genes commonly present in higher eukaryotes (CEGs) by Blastp comparisons (eVal: 10-3) 57 . All genomes were analysed using the PEDANT system 58 .
The genomes and transcriptomes of A. gallica and A. solidipes were sequenced using the Illumina platform at the Joint Genome Institute (JGI). Genomic DNA was sequenced as pairs of Illumina standard and Nextera long mate-pair (LMP) libraries. For the standard libraries, 500 ng of DNA was sheared to 270 bp using the Covaris E220 (Covaris) and size selected using SPRI beads (Beckman Coulter). The fragments were treated with end-repair, A-tailing, and ligation of Illumina compatible adapters (IDT, Inc) using the KAPA-Illumina library creation kit (KAPA biosystems). For LMP, 1 µ g of DNA was used to generate the library using the Nextera LMP kit (Illumina). DNA was fragmented and ligated with biotinylated linkers using the Tagmentation enzyme. The fragments were circularized via ligation followed by random shearing using the Covaris LE220 (Covaris). The mate pair fragments were purified using Strepavidin beads (Invitrogen) and treated with end repair, A-tailing, and ligation of Illumina adaptors. The final product was enriched with ten cycles of PCR.
For transcriptomics, stranded cDNA libraries were generated using the Illumina TruSeq Stranded RNA LT kit. Messenger RNA (mRNA) was purified from 1 µ g of total RNA using magnetic beads containing poly-T oligos, fragmented and reverse transcribed using random hexamers and SSII (Invitrogen), followed by second strand synthesis. The fragmented cDNA was treated with end-pair, A-tailing, adapter ligation, and ten cycles of PCR.
All libraries were quantified using KAPA Biosystem's next-generation sequencing library qPCR kit and run on a Roche LightCycler 480 real-time PCR instrument. Except for A. gallica LMP, the quantified libraries were prepared for sequencing on the Illumina HiSeq sequencing platform utilizing a TruSeq pairedend cluster kit, v3 (A. gallica) or v4 (A. ostoyae), and Illumina's cBot instrument to generate clustered flowcells for sequencing. Sequencing of the flowcells was performed on the Illumina HiSeq2000 or HiSeq2500 sequencer using a TruSeq SBS sequencing kit, v3 or v4, respectively, following a 2 × 150 indexed run recipe. Sequencing of the A. gallica LMP library was performed on the Illumina MiSeq sequencer using a MiSeq Reagent kit, v2, following a 2 × 150 indexed run recipe.
Genomic reads from each pair of libraries were QC filtered for artifact/process contamination and assembled together with AllPathsLG 59 . Illumina reads of stranded RNA-seq data were used as input for de novo assembly of RNA contigs, assembled into consensus sequences using Rnnotator (v. 3.4) 60 . Both genomes were annotated using the JGI Annotation Pipeline and made available via the JGI fungal portal MycoCosm 61 .
CEGMA and BUSCO were used to assess the completeness of the assemblies. We used BUSCO Version 3.0.1 with the lineage specific profile library basidiomycota_odb9 (species:selected 41:25 with 1.335 BUSCO groups) downloaded from http://busco.ezlab.org. We performed the whole analysis in Gene set (proteins) assessment mode. The results were combined in a folder and plotted with the script generate-plot.py.
Phylogenomic analysis. We used whole genomes of five Armillaria spp. and 22 Agaricales encompassing white and brown rot wood decayers, litter decomposers and ECM fungi; species from the Russulales and Boletales were included as outgroups (Supplementary Table 5). An all-versus-all protein BLAST was performed using mpiBLAST-1.6.0 62 with default parameters to find homologs based on similarity. Homologs were clustered into protein families using the Markov Cluster Algorithm 63, 64 with an inflation parameter of 2.0. For each cluster, a multiple sequence alignment was inferred using PRANK v.150803 65 , run with default settings. Subsequently, spuriously aligned regions were excluded with trimAl v1.4.r15 66 , with a gap threshold of 0.2. Next, we inferred ML gene trees for each cluster of a minimum of four proteins using FastTree v2.19 67 . FastTree was run with CAT/GAMMA20 model for rate heterogeneity and an improved LG model 68 for substitution probabilities. Gene trees were mid-point re-rooted using Phyutility v 2.2.6 69 .
Next, we screened the gene trees with a custom Perl script 70 for the presence of deep paralogs and inparalogs. Gene trees with deep paralogs were eliminated, while in trees with inparalogs, the paralog closest to the root was retained in the corresponding alignments. Finally, 835 clusters passed these criteria and had > 50 amino acid sites and higher than 60% taxon occupancy were concatenated into a supermatrix with 188,895 sites. We used this supermatrix to infer species tree with RAxML PTHREADS version 7.2.8 71 using a partitioned WAG+ G model, where each data partition represented a single input gene family. To assess branch support of the tree we performed a bootstrap analysis in RAxML with 100 replicates under the same model.
We used penalized likelihood (PL) implemented in the program r8s version 1.70 72 for molecular clock dating based on the optimal ML tree, two fossils and one secondary calibration. Archaeomarasmius legettii   73 was used to define minimum age constraint (92 Myr) for the origin of marasmioid clade (as used in a previous study 74 ). Palaeoagaricites antiquus (100-110 Myr) is the oldest fossil can be placed within Agaricales 75 , hence we constrained the MRCAs of this order to a minimum age of 108 Myr. To define the minimum age constraint of the origin of Boletales to 84 Myr, we referred to published analyses 27 . Maximum age constraints were defined with a wide safety window to allow for the calibrated clades to be inferred at least twice as old as their minimum ages. The MRCA of the Boletales, that of the Agaricales and that of the marasmioid clade were constrained between 84-250 Myr, 108-250 Myr and 92-180 Myr, respectively, as minimum and maximum ages. A cross-validation (CV) analysis was used to identify an optimal smoothing parameter between 10 −20 and 10 9 . All analyses were run in four replicates (random number seeds). We applied the additive penalty function and run optimization 25 times initialized from independent starting points. Furthermore after reaching an initial solution of an optimization step, the solution was perturbed and the truncated Newton (TN) optimization was rerun 20 times. We found that the optimal smoothing parameter (λ) is between 10 −7 and 10 −3 , thus we estimated divergence time with five λ ranges between 10 −7 and 10 −3 . By checking the gradients and the estimated ages we found a high similarity between the results of the analyses obtained with different values of λ ( Supplementary Fig.  2). By inspecting the CV analyses and the gradient values, λ = 10 −6 was chosen as the most adequate and was used to estimate divergence times.

COMPARE analysis.
We analysed the evolutionary history of gene families of Armillaria species and closely related Agaricales using the COMPARE pipeline 76 . To reconstruct gene duplications and losses, a genome-wide collection of 13,821 gene trees (see previous section) was first reconciled with the species tree using Treefix v1.1.10 77 . Treefix was run with RAxML site-wise likelihood model, Maximum Parsimony Reconciliation model (MPR) and an -alpha threshold of 0.001 to find any alternate gene tree topologies that minimize duplication/loss costs but are statistically equivalent to the ML gene tree. Orthologs were identified and recoded into a presence/absence matrix, as described in a previous study 76 . We then inferred duplications and losses for each orthogroup along the species tree using Dollo parsimony. Gene trees with less than four proteins were excluded. We also reconstructed the genome size for a given node by summing over gains and losses to the genome size of the MRCA. GO enrichment analysis based on the Fisher exact test with Benjamin-Hochberg correction was performed using Pfam domains mapped to protein clusters and creating GO annotations with pfam2go version 2016/10/01 78 followed by enrichment analysis at p < 0.05.

Copy number analyses.
We analysed the copy number distribution of selected gene families in Armillaria, Moniliophthora, and Heterobasidion annosum and, non-pathogenic members of Agaricales. For this, predicted proteomes of all the taxa under study were annotated with Pfam domains and IPR signatures. We scanned the protein sequences for pfams using pfamscan.pl v1.5 against the Pfam 30.0 database and InterPro signatures using InterProscan v5. Next, we built search terms for 42 PCWDE families and 25 putative pathogenicity related genes based on evidence from the literature ('Search_Terms' in Supplementary Table 5). In few cases, Pfam signatures were either absent or did not yield consistent results. CBM67 copy numbers were obtained based on BLAST hits by counting homologs obtained through blasting proteins annotated by JGI as CBM67. Genomes were searched iteratively for significant hits to avoid the impact of phylogenetic distance on the detection of homologs. We annotated cellobiose dehydrogenases (CDH, AA3_1), alcohol oxidases (Aox, AA3_3) and GMC oxidoreductases (GMT, AA3_2) using a tree-based approach. To this end, proteins of A. gallica (Armga1) were combined with homologs in other genomes, followed by multiple sequence alignment construction using MAFFT v7.273 79 with -auto option and the estimation of a gene tree with FastTree v2.19 (as above). Classifications were based on the occurrence as monophyletic groups in the phylogenetic tree. SSPs were defined as proteins with less than 300 amino acids and sequence-based evidence for secretion as inferred by signalP v4.1 80 . Both signalP-noTM and signalP-TM models were used for the detection of signal peptides.

NATuRE ECOlOgy & EvOluTiON
Analysis of TEs. De novo element discovery and annotation. We used the REPET package version 2.5 81,82 to identify, classify and annotate TEs within the genomes studied. Since the A. mellea genome contained large amounts of 'chaff ' contigs, we thresholded this assembly to only include contigs larger than 1,000 bp before further analysis.
The REPET de novo pipeline was run using genome self-alignment as well as a search for structural features using LTRHarvest 83 . Consensus sequences were clustered using Piler v1.0 84 , Recon 1.08 85 and Grouper 86 and classified using PASTEC 82 . The resulting consensus sequences were filtered for low-frequency elements, short sequence repeats and sequences that were classified as putative host genes and then clustered into families using Markov clustering 64 . For annotation, the consensus sequences collected from each species were combined into a panspecies TE library. TE identification was then carried out using the REPET anno pipeline, implementing BLASTER, RepeatMasker 4.06 and CENSOR 4.2.29 87 . Configuration files for the pipelines containing detailed parameter settings as well as a sample command line script can be found at https://github.com/JackyHess/ Armillaria_TE_annotations.
TE organization along the genome in core Armillaria species. To gain an understanding of how TEs are organized among the core Armillaria genomes we investigated TE and gene content density in 50 kb genome windows. Each genome was segmented into 50 kb partitions using bedtools (http://bedtools.readthedocs.io) and for each partition, TE coverage and genic coverage were estimated. Before estimating genic coverage, we filtered annotated CDS for which 20% of the sequence overlapped with TE annotations to remove putative TE-derived genes. This reduced the number of protein-coding genes considered by up to 7,000 (Supplementary Table 1). All analysis scripts can be found at https://github.com/ JackyHess/Armillaria_TE_annotations.
In-depth transcriptome analysis of A. ostoyae. Whole transcriptome sequencing was performed using TrueSeq RNA Library Preparation Kit v2 (Illumina) according to the manufacturer's instructions. Briefly, RNA quality and quantity measurements were performed using RNA ScreenTape and Reagents on TapeStation (all from Agilent) and Qubit (ThermoFisher); only high quality (RIN > 8.0) total RNA samples were processed. Next, RNA was DNaseI (ThermoFisher) treated and the mRNA was purified and fragmented. First strand cDNA synthesis was performed using SuperScript II (ThermoFisher) followed by second strand cDNA synthesis, end repair, 3'-end adenylation, adapter ligation and PCR amplification. All purification steps were performed using AmPureXP Beads (Backman Coulter). Final libraries were quality checked using D1000 ScreenTape and Reagents on TapeStation (all from Agilent). Concentration of each library was determined using the QPCR Quantification Kit for Illumina (Agilent). Sequencing was performed on Illumina NextSeq instrument using the NextSeq Series 2 × 150 bp high-output kit (Illumina) generating more tha 20 million clusters for each sample.
Bioinformatic Analysis Draft genome sequence together with genome annotation file was used as a reference for A. ostoyae RNA-seq analysis. Paired-end Illumina NextSeq reads were quality trimmed using CLC Genomics Workbench tool version 9.5.2 (CLC Bio/Qiagen) removing ambiguous nucleotides as well as any low-quality read end using an error probability cutoff value of 0.05 (corresponding to a Phred score of 13). Trimmed reads were mapped using the RNA-Seq Analysis 2.1 package in CLC allowing intergenic read mapping and requiring at least 80% sequence identity over at least 80% of the read lengths; strand specificity was omitted. Reads with less than 30 equally scoring mapping positions were randomly mapped; reads with more than 30 potential mapping positions were considered as uninformative repeat reads and were excluded from the analysis (Supplementary Table 6).
'Total gene read' RNA-seq count data was imported from CLC into R version 3.0.2. Genes were filtered on the basis of their expression levels, keeping only those features that were detected by at least five mapped reads in at least 25% of the samples included in the study. Subsequently, 'calcNormFactors' from package edgeR version 3.4.2 88 was used to perform data normalization based on the trimmed mean of M-values (TMM) method 89 . Log transformation was carried out by the 'voom' function of the limma package version 3.18.13 90 . Linear modelling, empirical Bayes moderation and the calculation of differentially expressed genes were carried out using limma. Genes showing at least four-fold gene expression change with an FDR-corrected p-value below 0.05 were considered as significantly differentially expressed. Multidimensional scaling ('plotMDS' function in edgeR) was also applied to visually summarize gene expression profiles to reveal similarities between samples. In addition, unsupervised cluster analysis with Euclidean distance calculation and complete-linkage clustering was carried out on the normalized data using the 'heatmap.2' function from R package 'gplots' .
LC-MS, data processing and interpretation. Peptide mixtures were analysed using a Thermo Fisher Q-Exactive mass spectrometer coupled to a Dionex RSLCnano for LC-MS/MS analysis. LC gradients operated from 3-40%B over 40 min, with data collection using a Top15 method for MS/MS scans 93 . LC-MS spectra chromatograms were analysed manually using rawMeat software. Raw files corresponding to the aforementioned spectra were then processed against an A. ostoyae predicted protein database using MaxQuant (version 1.3.0.5) and further filtered and visualized in Perseus (Version1.4.1.3) 95 . MaxQuant parameters were as described in another study 94 . Peptide intensity values were normalized to log 2 values in Perseus. Only samples represented in all replicates of a sample group were taken (n = 1-4). Proteins common to all stages were z-score normalised, averaged and hierarchical clustering was performed to generate heat maps and profile plots. Venn diagrams were prepared for all proteins (unique and abundant proteins inclusive) in the entire data set (http://bioinformatics.psb.ugent.be).
Proteomic analysis revealed the presence of 2,549 proteins, of which 39.7% (1,012 of 2,549) proteins were common to all A. ostoyae developmental stages and 25.2% (643 of 2,549) were unique to individual developmental stages. Proteins were unique to rhizomorphs (n = 286), vegetative mycelium (n = 163), young fruiting bodies (n = 97), stage 1 and 2 primordia (n = 73) and mature fruiting bodies (n = 24). Quantitative changes in protein abundance were also evident between different development stages, albeit to different extents. For instance, comparing rhizomorph to vegetative mycelium stages, individual proteins underwent increased (n = 184) and decreased (n = 125) abundance (total detected n = 2,190), whereas comparative analysis of young fruiting bodies against mature fruiting bodies revealed only eight proteins with increased abundance (total detected n = 1,203). Comparative analysis of all other developmental stages yielded intermediate changes in the abundance of specific proteins.
Motif discovery. We used 1 kb sequences of co-expressed genes (which are most likely to contain their promoter), located upstream of the start codon, including 5' UTR, to predict putative TFBS. We performed de novo motif discovery in our co-expressed gene sets using Weeder 2.0 38,96 . We created frequency files from promoter regions of all genes of A. ostoyae using the Weeder2.0 frequencymaker. Frequencies were counted on both strands for motif lengths 6, 8 and 10, allowing 1, 2 and 3 mismatches, respectively. Using the inferred frequency counts as reference, both strands of the sequences were scanned using Weeder, allowing the detection of a maximum of 100 motifs. Next, we grouped the motifs based on Pearson correlation co-efficient, aligned using the ungapped Smith-Waterman algorithm and clustered using UPGMA, all executed via the STAMP webtool 97 . Finally, we constructed familial binding profiles for each motif group and searched for matching DNA motifs in the JASPAR Core (Fungi) 2016 database using the TOMTOM webserver 98 . Corresponding author(s): Laszlo G. Nagy

Initial submission Revised version Final submission
Life Sciences Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form is intended for publication with all accepted life science papers and provides structure for consistency and transparency in reporting. Every life science submission will use this form; some list items might not apply to an individual manuscript, but all fields must be completed for clarity.

Data exclusions
Describe any data exclusions. RNA-Seq and proteomics: no samples were excluded

Replication
Describe whether the experimental findings were reliably reproduced.
Of our experiments, replication was needed for RNA-Seq and proteomic analyses: for these 3 biological replicates were done.

Randomization
Describe how samples/organisms/participants were allocated into experimental groups.
Samples were grouped based on developmental progression and tissue type. Armillaria has a well-known developmental progression from vegetative mycelia to mature fruiting bodies.

Blinding
Describe whether the investigators were blinded to group allocation during data collection and/or analysis.

N/A
Note: all studies involving animals and/or human research participants must disclose whether blinding and randomization were used.

Statistical parameters
For all figures and tables that use statistical methods, confirm that the following items are present in relevant figure legends (or in the Methods section if additional space is needed).

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement (animals, litters, cultures, etc.) A description of how samples were collected, noting whether measurements were taken from distinct samples or whether the same sample was measured repeatedly A statement indicating how many times each experiment was replicated The statistical test(s) used and whether they are one-or two-sided (note: only common tests should be described solely by name; more complex techniques should be described in the Methods section) A description of any assumptions or corrections, such as an adjustment for multiple comparisons The test results (e.g. P values) given as exact values whenever possible and with confidence intervals noted A clear description of statistics including central tendency (e.g. median, mean) and variation (e.g. standard deviation, interquartile range)

Clearly defined error bars
See the web collection on statistics for biologists for further resources and guidance.