Genome reduction and relaxed selection is associated with the transition to symbiosis in the basidiomycete genus Podaxis.

Insights into the genomic consequences of symbiosis for basidiomycete fungi associated with social insects remain sparse. Capitalizing on viability of spores from centuries-old herbarium specimens of free-living, facultative, and specialist termite-associated Podaxis fungi, we obtained genomes of 10 specimens, including two type species described by Linnaeus >240 years ago. We document that the transition to termite association was accompanied by significant reductions in genome size and gene content, accelerated evolution in protein-coding genes, and reduced functional capacities for oxidative stress responses and lignin degradation. Functional testing confirmed that termite specialists perform worse under oxidative stress, while all lineages retained some capacity to cleave lignin. Mitochondrial genomes of termite associates were significantly larger; possibly driven by smaller population sizes or reduced competition, supported by apparent loss of certain biosynthetic gene clusters. Our findings point to relaxed selection that mirrors genome traits observed among obligate endosymbiotic bacteria of many insects.

Introduction

Engaging in a symbiotic relationship provides a multitude of ecological and evolutionary advantages (Kiers and West, 2015; Oliver et al., 2010), allowing range expansions and occupation of new niches or providing an escape from competition and abiotic stress (Cornwallis et al., 2017; Oliver et al., 2010). These relationships can benefit both (mutualism) or only one (parasitism/commensalism) of the partners (Wilkinson, 2001). The longer and more specific symbioses are, the more dependent the partners are expected to become (Bennett and Moran, 2015; Ebert, 1998; Legros and Koella, 2010), often leading to genomic change driven by drift or selection (Bennett and Moran, 2015). This is best known in beneficial endosymbionts of insects where small effective population sizes and frequent bottlenecks during vertical transmission can result in extreme genome reduction, the loss of non-essential genes, and rapid protein evolution (McCutcheon and Moran, 2012). However, genomic consequences in fungal symbiosis are less well understood.

Obligate symbiotic mycorrhizal and plant pathogenic fungi can exhibit some of the largest and smallest genome sizes in fungi (Raffaele and Kamoun, 2012; Kohler et al., 2015; Lofgren et al., 2021; Reinhardt et al., 2021). Increases in genome size are generally associated with the proliferation of repetitive and non-coding regions of DNA (Raffaele and Kamoun, 2012), while reduced genome sizes can be the result of reduced transposon content or intron losses (Raffaele and Kamoun, 2012). Although gene gains and losses can correlate with genome expansions, this is not always the case (Raffaele and Kamoun, 2012; Kohler et al., 2015; Lofgren et al., 2021; Reinhardt et al., 2021). However, rapid genetic turnover and evolution in gene content and proteome composition, driven by drift or selection, is frequently reported in mycorrhizal and plant pathogenic fungi with both increased and reduced genome sizes (Kohler et al., 2015; Lofgren et al., 2021; Reinhardt et al., 2021).

Obligate social insect-fungal symbioses have independently evolved multiple times in ants, termites, beetles, and bees (Nygaard et al., 2016; Paludo et al., 2018; Poulsen et al., 2014; Vanderpool et al., 2018). In each of these, the fungi show morphological differences to non-symbiotic relatives and at least one of the partners possesses genes which the other lacks (Nygaard et al., 2016; Paludo et al., 2018; Poulsen et al., 2014). However, the paucity of well-characterized, closely related, free-living relatives makes it challenging to deduce genomic change resulting from the evolution of obligate symbioses in these lineages. In contrast, the basidiomycete genus Podaxis (Desv., Family: Agaricaceae; Figure 1A) covers the entire range of host dependence. In Sub-Saharan Africa and Australia, fruiting bodies of Podaxis grow from the mounds of grass-harvesting Trinervitermes (Holmgren; Africa) and Nasutitermes (Dudley; Australia) termites (both sub-family: Nasutitermitinae) (Conlon et al., 2016; Priest and Lenz, 1999), while free-living species are found in semi-arid grasslands and arid deserts in Africa, Asia, Australia, and the Americas (Conlon et al., 2019). The role of Podaxis with termites is unknown, but some specimens appear specialized and only produce fruiting bodies on termite mounds (Conlon et al., 2016). The genus Podaxis was long considered as a single species, but there is now ample evidence that these different life histories segregate into clades—likely corresponding to several species (Figure S1)—that either exclusively associate with termites (apparently obligate), are found on both termite mounds and free-living (apparently facultative), or are exclusively free living, often in deserts (Conlon et al., 2016).

Figure 1

Sampling of the 10 Podaxis strains with genome assembly and annotation statistics

(A) The Podaxis pistillaris type specimen LINN1287.7 (India) described by Linnaeus in 1771. Scale bar: 1 cm. (Podaxis photo: © The Linnaean Society of London; Linnaeus image from Wikimedia Commons).

(B) Podaxis specimen India (LINN1287.7) growing in culture at least 250 years after its collection (© B.H.C.).

(C) Genome assembly lengths varied significantly between all life histories with termite-specialist clades having smaller genomes than free-living. pANOVA: F = 14.13, p = 0.010; pairwise padj < 0.05.

(D) The number of annotated genes varied significantly between life histories, exhibiting the same pattern as the genome assembly lengths. pANOVA: F = 24.80, p = 0.006; pairwise padj < 0.05.

(E) Mitochondrial genomes were significantly longer in specialists than in free-living or facultative clades. pANOVA: F = 17.45, p = 0.005; free-living vs. facultative padj = 0.653, specialist vs. free-living/facultative padj < 0.012.

Beyond indications that spores from termite-specialist clades are smaller and with thinner walls than those of facultative or free-living clades (Conlon et al., 2016), nothing is known about potential genomic changes in these fungi as a consequence of living with termites or how these changes differ from free-living desert clades. Using a combination of genome analyses and bioassays guided by these analyses, we investigate genomic changes associated with the transition from a free-living life history to becoming a specialist termite symbiont.

Results and discussion

Germination of spores from herbarium specimens, including two described by Linnaeus

We obtained spores from dried fruiting bodies of Podaxis (Figure 1A) from the South African National Collection of Fungi (PREM), the Natural History Museum of Denmark (C-F), the Queensland Herbarium, Australia, (AQ) and the Linnaean Society of London (LINN) (Table S1). The specimens from the Linnaean Society represented the first two type specimens described in the genus (Podaxis pistillaris and Podaxis carcinomalis, respectively) by Linnaeus in 1771 (LINN1287.7, described as Lycoperdon pistillare, hereafter the Indian strain; Figure 1A) (Linnaeus, 1771) and 1781 (LINN1287.8, as Lycoperdon carcinomale, hereafter LINN Cape) (Linnaeus, 1781), making these specimens at least 250 and 240 years old, respectively. However, as LINN Cape was collected by Thunberg in the South African Cape region (Linnaeus, 1781), we place its likely collection date between 1772 and 1775: the years Thunberg was in South Africa. Due to the taxonomic confusion surrounding Podaxis (Conlon et al., 2016; Conlon et al., 2019), we refrain from using species names for all other studied strains.

Despite two of our specimens having been collected in the 18th century and the historical use of mercury as a pesticide in the herbarium collections (Kataeva et al., 2009), we successfully germinated spores from Podaxis strains across all known clades (Figures 1B and S1; Table S1), including the Linnean type specimens (Figure 1B). This reveals an extraordinary capacity for Podaxis spores to remain viable, despite prolonged periods without moisture, and suggests that spores could remain dormant in the environment for centuries before germinating once conditions allow. The dark spore pigmentation in Podaxis (Figures 1A and S2; Table S8) (Conlon et al., 2016) indicates heavy melanization, likely providing protection from environmental stressors (Eisenman et al., 2020), as melanization is known to positively correlate with spore longevity in basidiomycetes (Nguyen, 2018).

Genome reduction in the transition to a symbiotic lifestyle

From our successful cultures of Podaxis, we chose 10 specimens for whole-genome sequencing, based on their Internal Transcribed Spacer (ITS) sequences, to cover all six previously described (Conlon et al., 2016) and three new phylogenetic clades which spanned the different life histories: three free-living, three facultative, and three termite-specialist clades (Figure S1). Haploid assembly lengths were somewhat small (31–39 Mbp; Figure 1C; Table S2) for Agaricales (typically 35–65 Mbp) (Gupta et al., 2018). However, genome completeness was estimated to be >90% for all strains and >97% for assemblies <35 Mbp (Table S2). N50s varied from 61–140 Kbp (Table S2), and genomes sequenced from presumed homokaryons (heterozygosity minus anticipated sequencing error ≤0.001; Table S2) had significantly larger N50s than those of presumed heterokaryons (heterozygosity minus anticipated sequencing error >0.001; Table S2), even after reduction to a haploid assembly (linear model: F1 = 15.07, p = 0.005; Table S2). GC content ranged from 48.2–49.9% (Table S2), with no significant variation between life histories (phylogenetic analysis of variance [ANOVA], hereafter “pANOVA”: F = 1.293, p = 0.442).

Genomes of termite-specialist clades were significantly shorter than those of other clades (Figure 1C; pANOVA: F = 14.13, p = 0.010; pairwise padj < 0.05) and were well below the 35–65 Mbp assembly lengths typical for members of the Agaricales (Gupta et al., 2018), in contrast to all free-living and three out of four facultative strains. While there was much more variation in genome length between specimens from facultative clades, one of the specimens found on, or close by, a termite mound (Kuranda) had a longer assembly than the free-living Linnean specimen from India (Figure 1C; Table S2). The variation in genome size for strains from facultative clades (Figures 1C and S1; Conlon et al., 2016), where fruiting bodies are found both on termite mounds and free-living species on the ground, suggests that a stable, long-term association, seen in the specialist clades, is generally associated with genome reduction in Podaxis. These overall genome differences were also reflected in the number of annotated genes (Figure 1D; Table S2), which varied significantly with life history (pANOVA: F = 24.80, p = 0.006). While the difference in the number of genes was substantially larger for termite specialist vs. other life histories (pairwise padj < 0.012) than between facultative and free-living clades (pairwise padj = 0.036), the number of introns per gene did not differ significantly between life histories (Table S2; pANOVA: F = 2.75, p = 0.220; pairwise padj < 0.05). The length of coding sequences per gene (Table S2; pANOVA: F = 3.62, p = 0.168; pairwise padj: free living-facultative = 0.117, free living-specialist = 0.195, facultative-specialist = 0.646) and mean intron length (Table S2; pANOVA: F = 1.25, p = 0.472; pairwise padj: free living-facultative = 0.994, free living-specialist = 0.639, facultative-specialist = 0.639) did not differ significantly between life histories either; suggesting overall gene loss, rather than a reduction in coding sequence length or a decrease in intron content, is one of the factors driving genome reduction. Genome reduction, and subsequent gene loss, is common during transitions to obligate symbiosis (McCutcheon and Moran, 2012) and associated with increased specialization by a pathogenic (de Man et al., 2016; Schuelke et al., 2017) or mutualistic (Kooij and Pellicer, 2020; Nagendran et al., 2009) symbiont. Our findings are consistent with this and thus suggest that these clades are most likely obligate specialist symbionts of termite hosts.

Mitochondrial genome expansion in the transition to a symbiotic life history

Contrary to the whole-genome assemblies, mitochondrial genomes were significantly longer in specialists than in facultative and free-living specimens (Figure 1E; Table S3; pANOVA: F = 17.45, p = 0.005; pairwise padj = 0.012), which did not differ significantly from each other (pairwise padj = 0.653). This did not appear to be due to variation in assembly quality, as we successfully annotated 11–14 genes and 22 tRNAs for all but one strain (Moombra: 4 genes and 16 tRNAs; Table S3). Annotations were thus close to the predicted 15 core genes and the range of 21–35 tRNAs present in other Agaricomycetes (Xu and Wang, 2015). This, combined with the fact that five of the 10 genomes were circularized (Table S3), supports high assembly quality. While previous work has documented that other genera of fungi can exhibit interspecific variation of up to 16 Kbp in mitochondrial genome length (Sandor et al., 2018), this is still substantially below what we observe in Podaxis (Figures 1E and S3). The significant increase in mitochondrial genome size of termite specialists is conceivably driven either by reduced selection on the streamlining of genome size, relative to the scenario that free-living lineages experience (Conlon et al., 2016; Sandor et al., 2018; Selosse et al., 2001), or as a result of reduced effective population size for specialist lineages (Lynch et al., 2006).

A free-living lineage is sister to all other Podaxis

Phylogenomic analysis of predicted protein sequences for 2,339 single-copy orthologs placed the free-living Ethiopian strain as the sister group of all other strains and termite specialist clades as a terminal monophyletic group (Figures 2A and S4), suggesting that the genome reduction observed in these strains likely occurred in a single ancestor. The clade containing the type specimen for Podaxis rugospora (PREM 43879, Figure S1), represented by the Malagasy specimen (Figure 2A), was the sister lineage to the termite-specialist life history (Figure 2A). This contrasts previous ITS phylogenies that placed the specialist termite clades as sister to all other Podaxis (Figure S1) (Conlon et al., 2016). To verify the phylogenomic placement, we performed phylogenetic analyses using non-coding intron sequences from identified single-copy orthologs. Although the Podaxis rugospora clade (represented by the Madagascar specimen) became the sister lineage to the clades represented by Mexico, India, and Kenya, the intron phylogeny supported the placement of the Ethiopian clade and the monophyly of the termite-specialists (Figures 2A and S5).

Figure 2

Phylogenomic and comparative analysis of identified orthologs in Podaxis genomes

Single-copy ortholog (A) protein (left) and intron (right) sequence phylogenies support an origin of Podaxis in the deserts of Sub-Saharan Africa with ocean crossings to Australia and the Americas. Termite specialists form a monophyletic group with an apparent African origin. All nodes in the protein sequence tree had a bipartition support value of 1, and nodes in the intron sequence tree had a posterior probability of 1, unless specified within the tree. Termite-specialist clades exhibit a higher rate of protein but not intron sequence evolution compared to non-specialists. Leucoagaricus leucothites and Agaricus bisporus were used as outgroups with Coprinopsis cinerea (not shown here but see Figures S4 and S5). Full phylogenies available as Figures S4 and S5.

(B) Most unique orthogroups (present in 2 or more genomes) were found shared between the free-living and facultative clades while the termite-specialist clades had the second highest number of unique orthogroups.

The phylogenomic analyses support a Sub-Saharan African origin of Podaxis as all members of the clade containing the Ethiopian strain are from this region (Figure S1) (Conlon et al., 2016). This suggests that Podaxis was introduced to the Americas at least once, while the separation of the three Australian specimens into two separate, and non-monophyletic clades implies at least two independent introductions to Australia. Although it is not uncommon for taxa to have originated on the supercontinent Gondwana and, as a result, to be present in Africa, Australia, and the Americas (Beck et al., 2008), the separation of the African and Australian landmasses ∼160 Mya (Peace et al., 2020) and of Africa and South America ∼120 Mya (Peace et al., 2020) would require the last common ancestor (LCA) of Podaxis to have existed for more than 160 My. The divergence of the Agaricaceae (containing Podaxis) from the Psathyrellaceae (including Coprinopsis cinerea which roots our phylogenies) is estimated to be ∼125 Mya (He et al., 2019), making it extremely unlikely that Podaxis originated before the breakup of Gondwana, and provides strong evidence that free-living Podaxis crossed the Atlantic Ocean at least once.

The phylogenetic origin of the termite host (Nasutitermitinae) 20–30 Mya (Bourguignon et al., 2016) implies that the LCA of the termite-specialist clades existed long after the separation of Africa and Australia. The LCA of Trinervitermes and Nasutitermes is believed to have been in Africa with the earliest arrival to Australia 10–20 Mya (Bourguignon et al., 2016). While both phylogenies provide strong support for an African origin of termite specialization, how termite specialists came to Australia remains unclear. After their African origin, termite hosts are expected to have arrived in Australia via South America and Asia (Bourguignon et al., 2016). While it is possible that Podaxis reached Australia traveling with its termite hosts, the lack of termite-specialist Podaxis in Asia or South America suggests that Podaxis arrived in Australia independent of its hosts, potentially through spores carried by currents across the Indian Ocean (Lumpkin and Johnson, 2013). The evolution of termite association corresponded with an increase in the per base substitution rate of the predicted proteome compared to other Podaxis clades (Figure 2A) and to non-coding regions (Figure 2A). This implies rapid proteomic evolution, through adaptive change or reduced effective population sizes (Lynch et al., 2006) during the transition to a termite-specialist lifestyle.

In contrast to the derived termite-specialist clades, the free-living and facultative lifestyles are not monophyletic (Figure 2A). Our phylogenomic reconstruction shows that facultative clades have multiple independent origins and likely co-opt adaptations for survival in deserts to overcome defenses within termite mounds; something which is supported by their genomic similarity to free-living sister clades (Figures 1C–1E; Figure 2A). While both free-living and facultative clades may survive in deserts and semi-arid grasslands, they apparently rely on rainfall, or other ephemeral water sources, for mushroom formation (Conlon et al., 2019). The variation in water availability between life histories could thus cause massive differences in generation times between environments, such as the Namib Desert (Van Zinderen Bakker, 1975) with decades-long intervals between rain, semi-arid grasslands such as the South African Highveld with yearly rains (Mason and Tyson, 2000), and termite mounds with permanent access to groundwater (Lee and Wood, 1971). Increased evolutionary rates in environments with shorter generation times could result in population structures that help drive speciation (Hamilton et al., 1990; Huyse et al., 2005). Consequently, it is plausible that biogeography, reproductive rates, and overlapping species ranges between Podaxis and termite hosts collectively influence the evolution and designation of clades as facultative or free living.

Termite-specialist proteomes differ from other life histories in functional capacity

Orthology analysis revealed 1,110 orthogroups which were unique to Podaxis and present in at least two of the annotated genomes (Figure 2B; Table S4). The largest number of unique Podaxis-specific orthogroups (273) belonged to those shared between free-living and facultative clades, which is perhaps unsurprising given the higher number of predicted genes in these genomes. However, despite their significantly reduced protein-coding capacity, the specialist termite-associated clades exhibit the second highest number of unique Podaxis-specific orthogroups (67) and very few shared with the other life histories, indicating some divergent proteome evolution with the transition to become specialists (Conlon et al., 2016). Identification of carbohydrate-active enzymes (CAZymes) from these genes using Peptide Pattern Recognition (PPR) revealed significant variation in the total number of CAZymes (Table S5; pANOVA: F = 40.10, p = 0.002); a trend which was predominantly driven by significant differences between specialist and other life histories (pANOVA pairwise padj = 0.003; supported by a phylogenetically-weighted PCA, Figure S6).

Gene Ontology enrichment analysis of Podaxis-specific orthogroups, unique to each life history or combination of life histories (Figure 2B; Table S4; Table S6), revealed significant enrichment (pairwise padj = 0.05) of terms linked to cellular responses to stress (particularly oxidative stress) and resistance to antibiotics, temperature, and toxins in Podaxis-specific orthogroups unique to free-living and facultative clades (Table S6). Similarly, oxidoreductases that are involved cellular oxidative stress responses and represented by EC group 1 were significantly reduced (pANOVA: F = 29.53, p = 0.002) in termite specialists in our CAZyme analysis (Figure 3A). This was mainly due to a significant reduction in laccases (EC 1.10.3.2) (Figure 3B; pANOVA: F = 66.39, p < 0.001) but also other oxidoreductases (Figure 3C; pANOVA: F = 10.69, p = 0.016).

Figure 3

Termite-specialist Podaxis have fewer oxidoreductase enzymes and exhibit reduced tolerance to oxidative stress

(A) The total number of oxidoreductases was significantly reduced in termite-specialist clades compared to other life histories. pANOVA: F = 29.53, p = 0.002.

(B) Laccase exhibited the greatest variation between termite-specialist clades and the other life histories. pANOVA: F = 66.39, p < 0.001.

(C) The number of oxidoreductases still varied significantly, even when laccases were excluded from the analysis. pANOVA: F = 10.69, p = 0.016.

(D) Lignin was the only plant structural polymer with a significantly reduced number of degradative enzymes in termite-specialist clades. pANOVA: F = 66.39, pairwise padj < 0.001.

(E) All Podaxis clades could reduce manganese (shown here by a blue halo), compared to controls (C). While there was variation between clades, this did not correlate with life history. Photographs taken after 14 days (© M.S.L.S.) correspond to (i) Machattie, (ii) Ethiopia, (iii) Mexico, (iv) Madagascar, (v) Kenya, (vi) Kuranda, (vii) PREM Cape, (viii) LINN Cape, (ix) Moombra.

(F) Specialist Podaxis clades exhibit significantly reduced tolerance to elevated H2O2 concentrations (linear model: F2 = 17.19, p < 0.001; Tukey's test: padj < 0.05; Table S7).

Laccases are involved in a wide range of biotic processes including interactions between fungi and other microorganisms (Baldrian, 2004) and the neutralization of toxic phenolic compounds (De Fine Licht et al., 2013), including lignin (Youn et al., 1995). Categorization by target substrates of CAZymes (Table S5) confirmed that although genes coding for enzymes targeting lignin, cellulose, and hemicellulose were all reduced in absolute numbers in termite specialists, the only substrate for which the number of CAZymes differed significantly was lignin (Figure 3D; pANOVA: F = 66.39, pairwise padj < 0.001; Table S5; Figure S7). The loss of oxidoreductase-coding genes may be driven by relaxed selection due to reduced stress, caused by more consistent abiotic conditions (Conlon et al., 2016) or reduced competition due to lower microbial diversity in termite mounds compared to free-living environments (Lee and Wood, 1971). However, laccases and oxidoreductases may also have been reduced in numbers due to the consistent presence of dried grass within termite mounds, from which Podaxis grows (Conlon et al., 2016). Grass has some of the lowest lignin content of any land plant (Novaes et al., 2010), and a significant reduction in lignin-degrading oxidoreductase enzymes in specialist clades may thus represent relaxed selection on these genes due to dietary specialization on grass stores of the termite hosts.

Clades all exhibit lignin-degradative potential but differ in ROS tolerance

We investigated the role of oxidoreductase enzymes, including laccase, in lignin degradation and oxidative stress responses using in vitro growth assays. The presence of lignin-degrading enzymes, inferred by the oxidation of manganese, was apparent for strains from all life histories, suggesting that all Podaxis clades retain enzymatic capacity for lignin degradation (Figure 3E). By contrast, the growth rate of termite specialists, relative to growth on 0 mM H2O2, was significantly reduced compared to free-living and facultative specimens (linear model: F2 = 17.19, p < 0.001; Tukey's test: padj < 0.001) that did not differ significantly from each other (Tukey's test: padj = 0.318) (Figure 3F; Table S7). This suggests that the reduction in oxidoreductase enzymes is associated with reduced tolerance to oxidative stress and other cellular processes, rather than a loss of the potential for lignin degradation. With reactive oxygen species, such as H2O2, secreted by microorganisms during respiration, reduced growth in high H2O2 environments would represent a competitive disadvantage. This, combined with the reduced microbial diversity in Nasutitermitinae mounds (Lee and Wood, 1971), indicates that Podaxis likely benefits from its association with termites by escaping from competitors which, over time, appears to have led to reduced tolerance for elevated H2O2 levels in the termite specialists.

No evidence for variable secondary metabolite potential across clades

Secondary metabolites, encoded by Biosynthetic Gene Clusters (BGCs), play a vital role in the competitive arsenals of microorganisms, providing protection from biotic and abiotic stressors (Keller, 2015). With the evolution of the termite-specialist life history apparently reducing the competitive ability of Podaxis, we explored whether this also influenced their biosynthetic potential. We identified 15–21 BGCs per genome, irrespective of life history (pANOVA: F = 2.191, p = 0.289). Only six BGC networks were restricted to one or two of the three life histories (Figures 4, S8, and S9). Four of these clusters were missing from the termite-specialist clades, and one BGC was only identified in two of the three termite specialists. BGC lengths are frequently over 10 Kbp (Maplestone et al., 1992), challenging accurate predictions in silico, especially when low N50s cause them to bridge scaffolds. However, we did not identify a significant link between N50 and the number of predicted BGCs (linear model: F1 = 0.004, p = 0.954). That four BGC networks are missing in the termite-specialist lifestyle could thus indicate that they are beneficial when free living but become redundant inside termite mounds. This is consistent with the historic use of Podaxis as a medicinal fungus (Conlon et al., 2019), supporting that it produces a range of compounds likely involved in interactions with other microorganisms.

Figure 4

Biosynthetic gene cluster (BGC) network analysis of sequenced Podaxis

There was a high degree of similarity in predicted BGCs across life histories. Of the twenty networks containing BGCs from three or more strains, four networks (1 terpene, 1 NRPS-like, 1 PKS, and 1 NRPS; highlighted in bold) were restricted to free-living and facultative lifestyles, while one network (a terpene; in bold) was restricted to facultative and specialist lifestyles. One NRPS (also in bold) was only predicted from two of the termite specialist strains.

Podaxis appears to benefit from symbiosis with termites

The association of Podaxis with termites appears to carry some fitness benefits for the fungus. The interior of a termite mound is more stable than a desert, providing constant temperatures, humidity, nutrition, and protection from UV radiation (Conlon et al., 2016). This, combined with the chemical defenses of Nasutitermitinae (Šobotník et al., 2010) that reduce microbial diversity (Lee and Wood, 1971), conceivably frees Podaxis from a degree of abiotic stress and biotic competition. Although it appears that Podaxis could benefit from reduced biotic and abiotic stress through its association with termites, it is not clear if its presence could be beneficial for the termites. However, the production of antimicrobial secondary metabolites could represent a potential benefit (Conlon et al., 2019).

The presence of Podaxis in termite mounds has only previously been confirmed by the production of fruiting bodies, and it is thus unclear if Podaxis is permanently or transiently present (Conlon et al., 2016, 2019). Our analyses indicate that some Podaxis strains likely have an obligate association with termites, although it is unclear if the termites have an obligate association with Podaxis. The production of fruiting bodies implies that Podaxis is not generally vertically transmitted between termite generations but potentially has an obligate yet horizontally transmitted association similar to Termitomyces, the fungal symbiont farmed by termites in the Macrotermitinae. In this symbiosis, termite mating flights and colony foundation ensue at the start of the rainy season, and fungal fruiting bodies occur toward the end of the rainy season, when the first termite foragers leave the nest to obtain fungal spores to start their fungal gardens (Koné et al., 2011). This timing is similar to apparent Podaxis-Trinervitermes associations in South Africa, where mating flights occur at the start of the summer rains while the peak time for Podaxis fruiting body formation is toward the end of the summer (Figure S10).

The reduced size and coding capacity of genomes for the termite-specialist Podaxis clade, combined with morphological differences in spores (Conlon et al., 2016), could suggest an obligate association with termites that has persisted over evolutionary time (Bennett and Moran, 2015; Fisher et al., 2017; Fu et al., 2020; Hess et al., 2018; Nicks and Rahn-Lee, 2017; Veselská and Kolařík, 2015). A reduced genome size with the transition from facultative to obligate symbiosis has also been observed in the genomes of the attine ant obligate fungal symbiont (Kooij and Pellicer, 2020) and its specialist parasite (de Man et al., 2016). Like the termite hosts of Podaxis, ants are superorganismal (Boomsma and Gawne, 2018; Wheeler, 1911) with sterile and morphologically distinct worker castes. These castes can only pass on their genes indirectly by raising reproductive siblings, in a role comparable to the germline and the soma, meaning that selection shifts to act more at a colony than an individual level (Boomsma and Gawne, 2018; Wheeler, 1911). This places superorganismal symbiotic fungi, like Podaxis, in the unusual position of having an ectosymbiotic relationship to individual insects but an endosymbiotic relationship with the colony, and unit of selection, as a whole (De Fine Licht et al., 2014). It has been suggested that variation in the symbiotic gut microbiota among castes of superorganismal insect colonies may be analogous to the variation in microbiota observed between tissues or body regions of multicellular organisms (Sinotte et al., 2020). Our results indicate a similar analogy between apparently obligate fungal symbionts living within a superorganismal insect colony and the bacterial endosymbionts of non-superorganismal insects.

Limitations of the study

While the ability to germinate spores from centuries-old herbarium specimens can be a valuable tool for studying biodiversity, the quality and thoroughness of the accompanying descriptions can be highly variable, meaning specific environmental data are not available for all specimens. While it is clear that Podaxis can produce fruiting bodies on termite mounds and some clades reproduce exclusively on them, our understanding of the life histories of Podaxis and its termite hosts restricts us to speculating about the precise nature of the interaction at this time. While we accounted for phylogenetic similarity between our samples, the monophyly of the three termite-specialist strains implies that they may have been subject to phylogenetic constraints while the genome size and content of the free-living and facultative strains vary. Greater sampling and sequencing would therefore be beneficial to test the trends we uncover here. More research into the ecology of Podaxis and its hosts could help us understand the prevalence of Podaxis on termite mounds and whether partners receive further fitness benefits from the association.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Michael Poulsen (mpoulsen@bio.ku.dk).

Materials availability

Spores were sampled from The South African National Collection of Fungi (PREM; Agricultural Research Council at Roodeplaat, Pretoria, South Africa), The University of Copenhagen (C-F; Copenhagen, Denmark), The Linnean Society of London (LINN; London, United Kingdom), the Queensland Herbarium (AQ; Brisbane Botanic Gardens, Mt. Coot-tha, Toowong, Queensland, Australia) or sent after collection. Fungal cultures have been deposited in the culture collection at the Forestry and Agricultural Biotechnology Institute (FABI), University of Pretoria, South Africa.

Data and code availability

Sequence data have been deposited with NCBI under the BioProject number: PRJNA681736

New ITS sequences for the ITS SI phylogeny are deposited with GenBank under the accession numbers: MW430027-MW430047.

Experimental model and subject details

Podaxis sp.

Podaxis is a globally-distributed genus of Basidiomycete fungi, found in hot deserts, grassland and in association with termites from the Nasutitermitinae subfamily (Conlon et al., 2016). The taxonomy of Podaxis is unclear after its reduction to a single-species genus: Podaxis pistillaris (Morse, 1933). Despite the subsequent naming of new species, albeit without consideration of the older descriptions, the vast majority of herbarium specimens are still labeled as P. pistillaris (Conlon et al., 2019). We therefore refrain from using species names for our specimens, aside from those which are designated as species-types.

Spores (Table S1) were germinated by first heating them to 40°C for 1 hr before spreading them onto Yeast-Malt Extract Agar (YMEA: 4 g/L yeast extract, 10 g/L malt extract, 4 g/L dextrose, 20 g/L agar) and incubated at 40°C. For strains from Queensland, 50 μg Ampicillin, 100 μg Streptomycin (both dissolved in H2O) and 35 μg Chloramphenicol (dissolved in EtOH) were also added per 1 L of YMEA. Plates were then checked regularly and any mycelial growth sub-cultured. Cultures were confirmed as Podaxis through sequencing of the Internal Transcribed Spacer (ITS) region and comparison to sequences generated through direct DNA extractions from the spores. Pure Podaxis cultures were then stored at 25 or 40°C.

Method details

ITS barcoding

Aside from the Australian (AQ) specimens (Machattie, Kuranda & Moombra), PREM 57485 (PREM Cape), MADm22 (Madagascar), LINN1287.7 (India) and LINN1287.8 (LINN Cape), all sequences had been published previously (Conlon et al., 2016). DNA was extracted from spores and mycelium using a chelex protocol (200 μL Chelex 100 (Sigma Aldrich), vortexed then incubated at 99.9°C for 30 min). The ITS region was amplified in a Polymerase Chain Reaction (PCR). PCR reactions (total 25 μL) contained 8.4 μL sterile distilled water, 0.1 μL bovine serum albumin (BSA), 2 μL of template, 12.5 μL of VWR Red Taq DNA Polymerase Master Mix (VWR International, USA) and 1 μL each of the Basidiomycete-specific ITS primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS4B (5′-CAGGAGACTTGTACACGGTCCAG-3′) (Gardes and Bruns, 1993). PCR reactions were run with the following conditions: 94°C for 4 min followed by 35 cycles of 94°C for 30 s, 58°C for 30 s, and 72°C for 30 s with a final extension step at 72°C for 4 min. PCR products were checked using agarose gel electrophoresis and purified by either MSB Spin PCRapace (STRATEC Molecular, Germany) or Exosap-IT (Affymetrix Inc., USA) according to the manufacturer's instructions. Purified PCR products were sent to Eurofins MWG Operon (Ebersberg, Germany) for sanger sequencing with both forward and reverse primer.

Whole genome sequencing

DNA extractions for whole-genome sequencing were performed using a scaled-up CTAB extraction (Poulsen et al., 2014). Whole-genome sequencing was performed using a combination of 100bp/150bp paired-end shotgun (BGIseq/DNBseq) and long-read (PacBio sequel) sequencing by BGI (Table S2).

In vitro assays

In vitro growth assays were performed for nine of the sequenced strains (all but India). Resistance to oxidative stress was tested using Potato-Dextrose Agar (PDA) media (VWR, PA, USA) enriched with 5 mM (0.17 g/L), 20 mM (0.68 g/L) or 40 mM (1.36 g/L) H2O2. Plates were photographed after 14 days with mycelial area calculated using ImageJ version 1.53g (Schneider et al., 2012). In addition to this, manganese oxidation due to H2O2 and the presence of a manganese peroxidase (MnP) enzyme was inferred using PDA enriched with 10 mL/L of 10 mM Manganese chloride tetrahydrate (MnCl2 . 4H2O; CAS: 13,446-34-9) solution and 40 mg/L leucoberberlin blue (CAS: 52,748-86-4). Leucoberberlin blue reacts with Mn(III) and Mn(IV) but not Mn(II). The reduction of Mn(II) by MnP could therefore be inferred by the presence of a blue halo around the inoculated mycelium (Janusz et al., 2017).

Quantification and statistical analysis

ITS phylogenetic analysis

Forward and reverse sequences were aligned in Geneious prime version 2019.1.1 (Biomatters Ltd., New Zealand) using Geneious' own algorithm after primers and low-quality ends were trimmed. Spore and culture sequences were checked to ensure they were identical. Resulting ITS consensus sequences were aligned using MUSCLE version 3.8.425 (Edgar, 2004) with 12 iterations. The alignment was inspected, ends were trimmed and sequences were realigned. A maximum-likelihood phylogenetic analysis was performed using the IQ-Tree webserver version 1.6.12 (Hoang et al., 2017; Kalyaanamoorthy et al., 2017; Trifinopoulos et al., 2016). The best fitting model was identified by Modelfinder and ultrafast bootstrap support calculated using 10,000 replicates. The resulting tree was visualized using the ITOL online platform version 5.6.3 (Letunic and Bork, 2019).

Whole-genome assembly

Sequencing results were checked for quality using FastQC version 0.11.8 (Andrews, 2010) and MultiQC version 1.7 (Ewels et al., 2016). Kmer depth was calculated using Jellyfish version 2.2.10 (Marçais and Kingsford, 2011) and Kmer-based estimates of genome size, heterozygosity, sequencing error rate and repeat content generated using GenomeScope (Vurture et al., 2017). A hybrid de novo genome assembly, combining BGISeq and PacBio data, was performed using SPAdes version 3.13.0 (Bankevich et al., 2012) for three strains: Ethiopia, Kenya and PREM 57485. All other genomes, for which only shotgun sequencing data were generated, were assembled de novo using SPAdes version 3.13.0 (Bankevich et al., 2012) (Table S2).

As the initial analysis by GenomeScope had revealed some samples were not homokaryons (heterozygosity-sequencing error >0; Table S2) and potentially contained multiple genetically-different nuclei per cell (Kooij and Pellicer, 2020), we used redundans version 0.14a (Pryszcz and Gabaldón, 2016) to identify duplicate contigs and reduce them into haploid assemblies. For strains without PacBio data, we also used redundans to scaffold them against the closest genome with PacBio data. PacBio reference selection was based on the ITS phylogeny (Figure S1) (Conlon et al., 2016) and was as follows: LINN Cape – PREM Cape; India – Kenya; Madagascar – Kenya; Kuranda – Ethiopia; Moombra – Ethiopia; Machattie – Ethiopia and Mexico – Ethiopia.

Scaffolds were then aligned to the chromosome-level, assemblies for Agaricus bisporus H39 (NZ_AEOK00000000.1) (Sonnenberg et al., 2016) and Coprinopsis cinerea CC3 (AACS00000000.2) (Stajich et al., 2010) using Last version 1060 (Kiełbasa et al., 2011). Scaffolds which aligned to neither genome were then BLASTed against Fungi in the NCBI RefSeq genomic database (Camacho et al., 2009; O'Leary et al., 2015) using the discontiguous megablast algorithm with no filtering of species-specific repeats. Any scaffolds which did not produce results, or had a bit score lower than 50, were removed from the assembly using Filter_fasta in BBTools version 38.86 (Bushnell, 2014). Assembly quality scores and completeness were calculated using QUAST version 5.0.2 (Gurevich et al., 2013) and BUSCO version 4.0.2 with the agaricales_odb10 dataset (Seppey et al., 2019) and Augustus version 3.3.3 (Keller et al., 2011) trained using Coprinopsis cinerea CC3 (AACS00000000.2) (Stajich et al., 2010) (Table S2). Comparison of genome assembly size between clade types was performed in R version 4.0.2 (R Core Team, 2020) using a phylogenetic ANOVA (Results & Discussion) in the package phytools version 0.7-47 (Revell, 2012) with the ortholog-based phylogenomic tree read in using version 5.4-1 of the ape package (Paradis and Schliep, 2018) and significance defined as p ≤ 0.05. The link between N50 and heterozygosity was tested in R version 4.0.2 (R Core Team, 2020) using a linear model (Results & Discussion) with significance defined as p ≤ 0.05.

Mitochondrial genome assembly

An initial mitochondrial genome assembly was performed using Norgal version 1.0.0 (Al-Nakeeb et al., 2017) with the best candidate being selected based on a BLAST search. We then selected the specimen with the longest circular contig from Norgal (PREM Cape) as a reference and the longest sequence from each strain as a seed in novoPLASTY version 4.0 (Dierckxsens et al., 2016). Mitochondrial genomes were annotated in Geneious Prime (2019) based on 70% similarity to the annotations of A. bisporus (NZ_AEOK00000000.1) (Sonnenberg et al., 2016), Macrolepiota fuliginosa (NC_045202.1) and C. cinerea (AACS00000000.2) (Stajich et al., 2010). Podaxis mitochondrial genomes were aligned using Mauve version 1.1.1 (Darling et al., 2004) to identify Locally Co-linear Blocks (LCBs) (Figure S3).

Genome annotation

Genomes were annotated using Augustus version 3.3.3 (Keller et al., 2011) with Coprinopsis cinerea CC3 (AACS00000000.2) (Stajich et al., 2010) (“coprinus_cinereus” in Augustus) as a training dataset. Amino acid and coding sequences were extracted from the resulting gff file using the getAnnoFasta.pl script in Augustus. Annotation quality was estimated using BUSCO 4.0.2 with the agaricales_odb10 dataset (Seppey et al., 2019). Intron sequences were extracted using BEDTools getfasta version 2.29.2 (Quinlan and Hall, 2010), coding sequence and intron lengths were calculated using stats.sh in BBTools version 38.86 (Bushnell, 2014) (Table S2). Comparison of the number and lengths of predicted genes and introns between clade types was performed in R version 4.0.2 (R Core Team, 2020) using a phylogenetic ANOVA (Results & Discussion) in the package phytools version 0.7-47 (Revell, 2012) with the ortholog-based phylogenomic tree imported using version 5.4-1 of the ape package (Paradis and Schliep, 2018) and significance defined as p ≤ 0.05.

Orthology and phylogenomic analyses

Predicted protein sequences from the Augustus annotation, combined with the A. bisporus H39 (NZ_AEOK00000000.1) (Sonnenberg et al., 2016), C. cinerea CC3 (AACS00000000.2) (Stajich et al., 2010) and Leucoagaricus leucothites (LU_LEL_1.0) (Floudas et al., 2020) genome assemblies, were used as inputs to Orthofinder version 2.3.12 (Emms and Kelly, 2019) using BLAST (Camacho et al., 2009), Mafft (Katoh and Standley, 2013) version 7.455, FastTree (Price et al., 2010) version 2.1.10 and combined into a single species tree using STRIDE (Emms and Kelly, 2017) and STAG (Emms and Kelly, 2018). All phylogenetically-weighted analyses used the phylogeny produced by Orthofinder.

Intron sequences were extracted from single-copy orthologous genes identified by Orthofinder (Emms and Kelly, 2019) using GenomeTools version 1.5.10 (Gremme et al., 2013) and BEDTools version 2.29.2 (Quinlan and Hall, 2010), aligned using Mafft version 7.455 (Katoh and Standley, 2013) and per-gene phylogenetic analyses performed using IQ-Tree version 1.6.12 (Nguyen et al., 2014). Model prediction was performed using ModelFinder (Kalyaanamoorthy et al., 2017) with 1000 bootstrap replicates. Individual gene trees, containing 6 or more strains, were combined into a species tree using ASTRAL-Pro version 1.1.2 (Zhang et al., 2020).

The trees (Figure 2A) were visualized in the ape and dendextend packages, versions 5.3 (Paradis and Schliep, 2018) and 1.14.0 (Galili, 2015) respectively, in R version 3.6.3 (R Core Team, 2020). Calibration of the phylogenomic tree was performed using the chronos function in ape (Paradis and Schliep, 2018) with the split of C. cinerea from Agaricaceae dated at 120-130 Mya (He et al., 2019), the maximum age of the termite specialist clades as 30 My (Bourguignon et al., 2016) and the maximum age of the Australian termite-specialist strain as 20 My (Bourguignon et al., 2016). The results (Figure 2A) were visualized using phytools version 0.7-47 (Revell, 2012) with colors based on the viridis package version 0.5.1 (Garnier, 2018) in R version 3.6.3 (R Core Team, 2020).

Functional gene annotation

Functional gene annotation was performed using InterProScan version 5.40-77.0 (Mitchell et al., 2018) with annotation of Gene Ontology (GO) terms, Panther families and KEGG pathways turned on using the “goterms,” “iprlookup” and “pathways” options respectively. The InterProScan results were then incorporated into the Orthofinder results using Kinfin version 1.0 (Laetsch and Blaxter, 2017). Per-orthogroup Gene Ontology (GO) enrichment analyses for each combination of clade types were performed using Pfam domains from InterPro (Mitchell et al., 2018) and Kinfin (Laetsch and Blaxter, 2017) in dcGO (Fang and Gough, 2012) with significance defined as padj ≤ 0.05.

CAZyme analysis

Functional annotation of carbohydrate-active enzymes was performed using HotPep version 1 (Busk et al., 2017) with default parameters. CAZymes for which specific EC numbers could be identified were manually annotated with their substrate using the ExPASy database (Artimo et al., 2012) and classified by the substrates presence in plant, bacterial or fungal cell walls. Comparisons of CAZymes and target substrates between clade types were performed using a phylogenetic ANOVA in the package phytools version 0.7-47 (Revell, 2012) in R version 3.6.3 (R Core Team, 2020). For the comparison of substrate types, adjusted p values were calculated using the false-discovery rate method in the p.adjust function with significance defined as p/padj ≤ 0.05.

The growth of Podaxis on elevated H2O2 concentrations was compared to the mean growth at 0 mM H2O2 (Table S7) using a linear model in R version 3.6.3 (R Core Team, 2020) with Tukey's post-hoc test (Results & Discussion). Significance was defined as p/padj ≤ 0.05.

Biosynthetic gene cluster analyses

Biosynthetic Gene Clusters (BGCs) were predicted using fungiSMASH version 5.1.1 (Blin et al., 2019). A potential link between the total number of predicted gene clusters and genome N50 values was tested using a linear model in R version 4.0.2 (R Core Team, 2020) with significance defined as p ≤ 0.05 (Results & Discussion). GenBank files from the fungiSMASH analysis were then further analyzed using the BiG-SCAPE BGC network prediction software (Navarro-Muñoz et al., 2020). Network data was illustrated using Cytoscape 3.8.0 (Shannon et al., 2003). Predictive phylogenetic frameworks were built from characterized fungal terpene cyclases and polyketide synthases (Key Resources). Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the JTT model, and then selecting the topology with superior log likelihood value. Evolutionary history was inferred by using the Maximum Likelihood method and Le_Gascuel_2008 model (Le and Gascuel, 2008). A discrete Gamma distribution was used to model evolutionary rate differences among sites: 5 categories (+G, parameter = 4.4685). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 1.81% sites). All positions containing gaps and missing data were eliminated (complete deletion option) and the analyses were conducted in MEGA X (Kumar et al., 2018).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
AQ553608	The Queensland Herbarium	AQ553608
AQ553622	The Queensland Herbarium	AQ553622
Kuranda	The Queensland Herbarium	AQ645840
AQ645843	The Queensland Herbarium	AQ645843
Moombra	The Queensland Herbarium	AQ795752
Machattie	The Queensland Herbarium	AQ799166
AQ799332	The Queensland Herbarium	AQ799332
AQ799334	The Queensland Herbarium	AQ799334
Ethiopia	The University of Copenhagen	C-F92630
Mexico	The University of Copenhagen	C-F101400
Kenya	The University of Copenhagen	C-F101401
India – Podaxis pistillaris	The Linnean Society of London	LINN1287.7
Linn Cape – Podaxis carcinomalis	The Linnean Society of London	LINN1287.8
PREM44075	The South African National Collection of Fungi	PREM44075
PREM44294	The South African National Collection of Fungi	PREM44294
PREM47477	The South African National Collection of Fungi	PREM47477
PREM47484	The South African National Collection of Fungi	PREM47484
PREM54389	The South African National Collection of Fungi	PREM54389
PREM Cape	The South African National Collection of Fungi	PREM57485
PREM58760	The South African National Collection of Fungi	PREM58760
PREM60320	The South African National Collection of Fungi	PREM60320
PREMN833	The South African National Collection of Fungi	PREMN833
PRUM3655	The University of Pretoria	PRUM3655
PRUM3913	The University of Pretoria	PRUM3913
PRUM4335	The University of Pretoria	PRUM4335
Madagascar	Direct collection (Bryn Dentinger)	N/A
Chemicals
Ampicillin	Sigma-Aldrich	CAS: 69-53-4
Streptomycin	Sigma-Aldrich	CAS: 3810-74-0
Chloramphenicol	Sigma-Aldrich	CAS: 56-75-7
Bovine Serum Albumin (BSA)	Sigma-Aldrich	CAS: 9048-46-8
Red Taq DNA Polymerase Master Mix	VWR International	Catalog number: 733-
Manganese chloride tetrahydrate	Sigma-Aldrich	CAS: 13446-34-9
Leucoberberlin blue	Sigma-Aldrich	CAS: 52748-86-4
Critical commercial assays
MSB Spin PCRapace	STRATEC Molecular	Catalog number: 1020220300
Exosap-IT	Affymetrix Inc.	Product number: 78200
Deposited data
ITS sequence data	NCBI – GenBank	MW430027-MW430047
Whole-genome sequence data	NCBI – BioProject	PRJNA681736
PCR primers
Primer ITS1F: 5’– CTTGGTCATTTAGAGGAAGTAA–3’	(Gardes and Bruns, 1993)	N/A
Primer ITS4B: 5’–CAGGAGACTTGTACACGGTCCAG–3’	(Gardes and Bruns, 1993)	N/A
Software
Geneious Prime	Geneious Inc.	2019.1.1
MUSCLE	(Edgar, 2004)	3.8.425
IQ-Tree with Modelfinder	(Hoang et al., 2017; Kalyaanamoorthy et al., 2017; Trifinopoulos et al., 2016; Nguyen et al., 2014)	1.6.12
ITOL online platform	(Letunic and Bork, 2019)	5.6.3
FastQC	(Andrews, 2010)	0.11.8
MultiQC	(Ewels et al., 2016)	1.7
Jellyfish	(Marçais and Kingsford, 2011)	2.2.10
GenomeScope	(Vurture et al., 2017)	N/A
SPAdes	(Bankevich et al., 2012)	3.13.0
Redundans	(Pryszcz and Gabaldón, 2016)	0.14a
Last	(Kiełbasa et al., 2011)	1060
BBTools	(Bushnell, 2014)	38.86
QUAST	(Gurevich et al., 2013)	5.0.2
BUSCO	(Seppey et al., 2019)	4.0.2
Augustus	(Keller et al., 2011)	3.3.3
R	(R Core Team, 2020)	3.6.3 and 4.0.2
Phytools	(Revell, 2012)	0.7-47
Ape	(Paradis and Schliep, 2018)	5.4-1
Norgal	(Al-Nakeeb et al., 2017)	1.0.0
novoPLASTY	(Dierckxsens et al., 2016)	4.0
Mauve	(Darling et al., 2004)	1.1.1
BEDtools	(Quinlan and Hall, 2010)	2.29.2
Orthofinder	(Emms and Kelly, 2019)	2.3.12
BLAST	(Camacho et al., 2009)	2.5.0+
Mafft	(Katoh and Standley, 2013)	7.455
FastTree	(Price et al., 2010)	2.1.10
STRIDE	(Emms and Kelly, 2017)	N/A
STAG	(Emms and Kelly, 2018)	N/A
GenomeTools	(Gremme et al., 2013)	1.5.10
ASTRAL-Pro	(Zhang et al., 2020)	1.1.2
dendextend	(Galili, 2015)	1.14.0
Viridis	(Garnier, 2018)	0.5.1
InterProScan	(Mitchell et al., 2018)	5.40-77.0
Kinfin	(Laetsch and Blaxter, 2017)	1.0
dcGO webserver	(Fang and Gough, 2012)	N/A
HotPep	(Busk et al., 2017)	1
ExPASy	(Artimo et al., 2012)	N/A
fungiSMASH	(Blin et al., 2019)	5.1.1
BiG-SCAPE BGC network prediction software	(Navarro-Muñoz et al., 2020)	N/A
Cytoscape	(Shannon et al., 2003)	3.8.0
MEGA X	(Kumar et al., 2018)	X
Other
Podaxis ITS sequences	Conlon et al., 2016	NCBI: KT844852.1 – KT844880.1
Agaricus bisporus H39 genome assembly	(Sonnenberg et al., 2016)	NCBI: NZ_AEOK00000000.1
Coprinopsis cinerea CC3 genome assembly	(Stajich et al., 2010)	NCBI: AACS00000000.2
Macrolepiota fuliginosa mitochondrial genome assembly	NCBI	NCBI: NC_045202.1
Leucoagaricus leucothites genome assembly	(Floudas et al., 2020)	NCBI: LU_LEL_1.0
Penicillium roqueforti	(Agger et al., 2009)	NCBI: AAA33694.1
Aspergillus terreus		NCBI: AAF13263.1
Fusarium sporotrichioides		NCBI: AAD13657.1
Coprinus cinereus		NCBI: XP 001832573
		NCBI: XP 001832925
		NCBI: XP001832925
		NCBI: XP 001836356
		NCBI: XP 001834007
		NCBI: XP 001832548
Botrytis cinerea		NCBI: ATZ56107.1
Omphalotus olearius	(Wawrzyn et al., 2012)	JGI model name: MUStwsD_GLEAN_10001317
		JGI model name: MUStwsD_GLEAN_10003938
		JGI model name: MUStwsD_GLEAN_10005581
		JGI model name: MUStwsD_GLEAN_10000810
		JGI model name: MUStwsD_GLEAN_10000811
		JGI model name: MUStwsD_GLEAN_10003820
		JGI model name: MUStwsD_GLEAN_10000831
		JGI model name: MUStwsD_GLEAN_10000543
		JGI model name: MUStwsD_GLEAN_10000292
Fomitopsis pinicola		NCBI: EPT00921.1
Stereum hirsutum	(Quin et al., 2013)	NCBI: XP_007306912.1
		NCBI: XP_007299839.1
		NCBI: XP_007301467.1
		NCBI: XP_007308318.1
		NCBI: XP_007299456.1
Armillaria gallica	(Engels et al., 2011)	UniProt: P0DL13
Coniophora puteana	(Mischko et al., 2018)	NCBI: XP_007771895.1
		NCBI: XP_007765978.1
Boreostereum vibrans	(Zhou et al., 2016)	NCBI: AMW90891.1
Aspergillus flavus	(Lackner et al., 2012)	NCBI: AAS90093.1
Aspergillus nidulans		NCBI: Q5B0D0.1
		NCBI: Q5AUX1.1
Aspergillus ochraceoroseus		NCBI: ACH72912.1
Aspergillus terreus		NCBI: Q9Y8A5.1
		NCBI: Q9Y7D5.1
Chaetomium chiversii		NCBI: C5H882.1
Cercospora nicotianae		NCBI: Q6DQW3.1
Aspergillus ruber	(Nies et al., 2020)	NCBI: EYE95336
Neurospora crassa	(Zhao et al., 2019)	NCBI: XP_965600
Trichoderma virens	(Liu et al., 2019)	NCBI: XP_013952638
Agaricomycetes sp. BY1	(Braesel et al., 2017)	NCBI: APH07629.1
		NCBI: APH07628.1
	(Brandt et al., 2017)	NCBI: ART89046.1
		NCBI: ART89047.1
Antrodia cinnamonea	(Yu et al., 2016)	NCBI: AST08390.1
Armillaria mellea	(Lackner et al., 2012)	NCBI: AFL91703.1
Coprinopsis cinerea		NCBI: XP_001835415.2
Strobilurus tenacellus	(Nofiani et al., 2018)	NCBI: ATV82110.1