Vertical and horizontal gene transfer shaped plant colonization and biomass degradation in the fungal genus Armillaria

Abstract

The fungal genus Armillaria contains necrotrophic pathogens and some of the largest terrestrial organisms that cause tremendous losses in diverse ecosystems, yet how they evolved pathogenicity in a clade of dominantly non-pathogenic wood degraders remains elusive. Here we show that Armillaria species, in addition to gene duplications and de novo gene origins, acquired at least 1,025 genes via 124 horizontal gene transfer events, primarily from Ascomycota. Horizontal gene transfer might have affected plant biomass degrading and virulence abilities of Armillaria, and provides an explanation for their unusual, soft rot-like wood decay strategy. Combined multi-species expression data revealed extensive regulation of horizontally acquired and wood-decay related genes, putative virulence factors and two novel conserved pathogenicity-induced small secreted proteins, which induced necrosis in planta. Overall, this study details how evolution knitted together horizontally and vertically inherited genes in complex adaptive traits of plant biomass degradation and pathogenicity in important fungal pathogens.

Figure 1. Genome statistics and reconstruction of ancestral genome sizes for 15 Armillaria species and 5 Physalacriaceae outgroups.

Numbers at nodes represent ancestral proteome sizes in the Physalacriaceae tree. Purple circles correspond to sizes for each node in the tree (for gene gains and losses at each node, and for the complete species tree, see Extended Data Fig. 1A). Blue circles represent BUSCO scores for each species. For genome sizes, darker regions show TE sizes (Mbp) from the primary scaffolds, and lighter color shows the genome assembly sizes in Mbp (for TE categories, see Extended Data Fig. 1B). For proteome sizes, darker color shows proteins with no known functional domains (unannotated proteins). Carbohydrate active enzymes darker color shows plant cell wall degrading enzymes, lighter color shows other CAZymes. Secretomes and small secreted proteins - darker color shows unannotated proteins, lighter color shows proteins with known functional domains (for general genome statistics of all 66 species in Dataset 1, see Supplementary Table 1). Clade names are based on Koch et al, 2017.

Figure 2. Plant biomass degradation related genes in Armillaria.

a) Phylogenetic PCAs for PCWDE gene families. Species abbreviations are shown only for Physalacriaceae species (all species names and PCA loadings are given in Extended Data Fig. 3 and Supplementary Table 3). b) Boxplot of copy numbers of 16 CAZy orthogroups co-enriched in Physalacriaceae (n=20) and in Ascomycota (n=21) with respect to white rot (n=23) and litter decomposer fungi (n=24). The box plots shows the median and interquartile range, with the upper whiskers extending to the largest value from the 75th percentile, and lower whiskers extending to the smallest value from the 25th percentile. Scale limits for the boxplot were set to 14, losing one sample point (Conioc1 in OG0000781 with 18 genes). Lifestyles of the species used are denoted by color.

Figure 3. Horizontal gene transfers into Armillaria and the Physalacriaceae family.

a) Barplot showing the number of phylogenetically validated HT genes identified in each genome. b) Schematic representation of donor and recipient relationships for HT events after phylogenetic validation. Size of the arrows are proportional to the number of HGT events inferred from all the nodes belonging to a specific donor group. The number of events for each donor group are listed along the arrows. c) Violin plot showing expression dynamics of phylogenetically validated HT genes and vertically transferred (VT) genes. X-axis shows the sample names and the y-axis shows log₂ FPKMs (For expression dynamics of HT and VT genes in developmental dataset, see Extended Data Fig. 4).

Figure 4. Enrichment of differentially expressed genes of wood-decay, pathogenicity, stress-response and other gene families in 6 RNA-Seq datasets.

The heatmap shows enrichment ratios for 24 gene groups (x-axis) from aggregated differential gene expression data across 6 experiments (a - upregulated, b - downregulated genes). Enrichment ratios were calculated by considering the number of DEGs in a given gene group in the given sample, all DEGs in that sample as well as all genes in the genome and all genes in that gene group. The total number of DEGs in each experiment is shown as a barplot at right. In the heatmaps, warmer colors mean higher enrichment ratios (for a complete list of enrichment ratios Supplementary Table 5).

Figure 5. Pathogenicity-induced SSPs (PiSSPs) of A. luteobubalina induce cell death in host plants.

a) Heatmap shows log₂ fold changes for SSPs significantly upregulated in at least one time point. Red shows higher and blue depicts lower logFC, followed by presence/absence matrix of homologs of unannotated SSPs in 131 species (Dataset 2). The two experimentally validated PiSSPs are shown by red arrows. b) Transient transformation of the non-host N. benthamiana and the host E. grandis with negative controls (empty vector; Pisolithus microsporus proteins Pismi_683611 and Pismi_689601), with a positive control (BAX) or with an in planta expression vector encoding Armlut_1165297 or Armlut_1348401. As further controls, frameshift versions of the A. luteobubalina sequences were included (Armlut_1165297_fs or Armlut_1348401_fs).