Functional potential of the AOD lesion microbiome and identification of two abundant Gram-positive bacteria through metagenome binning
We have previously described the taxonomic composition of the AOD microbiome and the shift from the microbiome of healthy trees [4]. However, in order to identify the specific mechanistic processes that mediate lesion formation, a deeper functional analysis comparing the AOD microbiome and healthy microbiome was required. Thus, we aimed to identify categories of genes associated with AOD virulence and highlight the dramatic shifts taking place on a functional level. Samples were taken from three non-symptomatic Q. robur trees in Attingham Park, along with two symptomatic AOD Q. robur trees with drier lesions, and two symptomatic Q. petraea trees at Hill court park with fresher lesions where fluid was actively seeping from the lesions. As host (oak) DNA would dominate the sequencing libraries obtained using DNA extracted from host tissue, significantly reducing the microbiome signal, the host DNA was depleted using an NEB Next Microbiome enrichment kit in order to focus on the microbiota for this stage of the analysis. The functional potential of the microbiome was initially identified using MG-RAST [19] (Additional file 1: Table S1, Additional file 2: Table S2, Additional file 3: Table S3, Additional file 4: Table S4, Additional file 5: Table S5, Additional file 6: Table S6, Additional file 7: Table S7 and Additional file 8: Table S8), revealing that approximately 67–95% of the predicted genes in AOD lesion microbiomes were bacterial, whereas only 0.6–6% of genes in non-symptomatic samples were derived from bacteria. Functional microbiome analysis of SEED subsystem gene category abundances (functional gene categories) demonstrated a clear distinction between symptomatic and non-symptomatic samples, corroborating previous investigations of the AOD microbiome shift, and providing a genetic database for our holistic multi-omic analysis (Fig. 1a) [2,3,4, 6]. Specifically, the MG-RAST analysis demonstrated that functional SEED subsystem gene category abundances, associated with normal plant activity were significantly more abundant in non-symptomatic tree samples, while those associated with bacterial activity, bacteriophage activity and plant defence were more abundant in symptomatic trees (Additional file 9: Table S9, Fig. 1b).
The metagenome databases of annotated genes were combined revealing 627 distinct genes identified in all symptomatic samples, but undetected and thus considered absent in the non-symptomatic samples (Additional file 10: Table S10). These genes were 99% bacterial and included proteins involved in virulence (membrane transport, PCWDEs, iron scavenging, motility and chemotaxis, stress responses and regulation). In contrast, 223 distinct genes (80% plant associated, 20% fungal), involved mainly in general metabolism, were identified in all non-symptomatic samples, but in none of the symptomatic samples (Additional file 11: Table S11).
Although several new species of bacteria have been isolated from AOD lesions, it was important to attempt to identify any as-yet uncultivated microorganisms in the AOD lesion microbiome that may play an important role in lesion formation. Using metagenomic binning, we generated two draft genomes belonging to the genera Clostridioides and Carnobacterium, and other bacteria (Additional file 12: Table S12), from metagenomic reads, thus extending the AOD-associated microbiome of putatively unculturable organisms for further analysis. The predicted genome sizes were 5.45 and 2.42 Mbp and the genome completions were 97.5 and 66.7%, respectively for the Clostridioides and Carnobacterium bacteria, and both exhibiting < 5% contamination. However, as no ribosomal RNAs were identified in the Clostridioides genome, both genomes were determined to be medium quality drafts as outlined in the guidelines by Bowers et al. [20]. These genomes were annotated and found to contain PCWDE, secretion systems, catalases and ROS-defence associated genes, flagella and other virulence-associated genes (Additional file 12: Table S12 and Additional file 13). As a contrast to large public databases such as Swissprot where our newly identified bacterial species are not yet represented, the annotated genomes were combined into a narrowed-down database for specified bioinformatic analysis, containing the annotated oak transcriptome from NCBI (Additional file 14: Table S13), the annotated genomes of AOD-associated bacteria B. goodwinii, G. quercinecans, R. victoriana and other bacteria whose genomes had homology to AOD lesion metagenome coding domains; Escherichia coli, Dickeya dadantii, Pectobacterium carotovorum, Erwinia billingae, Serratia marcescens and Clavibacter michiganensis (Additional file 13). This database allowed us to specifically clarify the activity of each bacterial organism of interest in the AOD microbiome, as well as specify which host genes and proteins are active.
Metatranscriptome analysis reveals active host defences and the primary active pathogens of the AOD microbiome
We performed metatranscriptome analysis of symptomatic and non-symptomatic inner bark tissue to identify the functional activity of the host and its microbiota and to specifically determine the function and relative role of B. goodwinii, G. quercinecans, R. victoriana and other bacteria in the AOD lesion microbiome. The metatranscriptomes were profiled and compared across samples (Additional file 15: Table S14, Additional file 16: Table S15, Additional file 17: Table S16, Additional file 18: Table S17, Additional file 19: Table S18, Additional file 20: Table S19 and Additional file 21: Table S20). In symptomatic tissue, 11–21% of annotated predicted transcripts were bacterial (2–3% in non-symptomatic tissue) (Additional file 15: Table S14, Additional file 16: Table S15, Additional file 17: Table S16, Additional file 18: Table S17, Additional file 19: Table S18, Additional file 20: Table S19 and Additional file 21: Table S20). These data demonstrate that the genetic shift in the metagenome also translates into an enhanced functional activity of bacteria.
We identified 216 transcripts present in all AOD symptomatic samples but absent in non-symptomatic trees (based on annotation, Additional file 22: Table S21). Of the symptomatic transcripts, 94% were associated with bacteria, partly involved in virulence (regulation, membrane transport, signalling, biofilm, sporulation and PCWDEs), the remaining 6% were associated partly with plant defence and cell wall synthesis, revealing increased bacterial phytopathogenic activity in symptomatic tissue. Conversely, 263 transcripts were identified across all non-symptomatic samples, but not in any symptomatic sample (97% eukaryotic, 3% bacterial, Additional file 23: Table S22).
In order to generally assess quantitative differences between non-symptomatic and symptomatic metatranscriptomes, we performed a gene expression analysis, identifying 6419 differentially expressed genes in symptomatic tissue (Additional file 24: Table S23). The 3064 downregulated genes in symptomatic tissue mainly displayed homology to plants (69%), fungi (16%), animals (9%) and bacteria (5%), while the 3355 upregulated (upregulated is used generally here to include bacterial genes that may for example be unrepresented in non-symptomatic tissue) displayed homology to bacteria (61%), plants (32%), animals (4%) and fungi (2%). Many of the upregulated bacterial genes were associated with virulence, for example, biofilm formation, chemotaxis and motility, effectors, efflux pumps, PCWDEs, virulence regulation, reactive oxygen species (ROS) defence, protein secretion systems and toxin production (Fig. 2a). Upregulated plant genes associated with symptomatic trees included calmodulin binding and production, disease resistance, hormonal signalling, PCWDEs, membrane receptors, ROS production and protection and the WRKY superfamily stress response transcription factors (Fig. 2b). Furthermore, 31 downregulated and 66 upregulated genes were identified as originating from viruses. The majority of downregulated viral genes belonged to plant viruses, while the majority of upregulated viral genes belonged to bacteriophages (Additional file 24: Table S23).
To further clarify the host response in the symptomatic tissue, we performed a geneset enrichment analysis (GSEA) on the expression data homologous to Arabidopsis genes (Fig. 2c) [21, 22]. GSEA demonstrated upregulation of 190 genesets and downregulation of 276 genesets (Additional file 25: Table S24). Upregulated genesets were partly associated with defence and wounding response, ROS burst, cell wall modification and cell death. Downregulated genesets were partly associated with organelle organisation, circadian rhythm, protein transport, cell-cell signalling, photoperiodism and regulation of development, demonstrating a redistribution of resources in the host from growth to defence.
We developed a narrowed-down database containing bacteria identified as abundant and important in AOD lesion formation in our previous study [4], along with bacterial draft genomes extracted from the metagenomic data, and the Q. robur transcriptome available in the NCBI database, which we annotated (Additional file 14: Table S13) [23]. The narrowed-down database allowed for specifying the representation of all bacteria of interest and the oak host in AOD lesions the metatranscriptomes and metaproteomes, and comparing symptomatic against non-symptomatic field samples (Additional file 26: Table S25, Additional file 27: Table S26, Additional file 28: Table S27, Additional file 29: Table S28, Additional file 30: Table S29, Additional file 31: Table S30 and Additional file 32: Table S31). We could significantly detect (FDR < 0.05, fold change > 2) 499 upregulated genes in symptomatic samples and 92 downregulated genes. Of the 499 upregulated genes, 295 were B. goodwinii genes, 198 were oak genes and 6 belonged to the predicted bacterium of the genus Clostridioides (Additional file 33: Table S32). However, studying the identified gene transcripts of each sample, we could detect PCWDEs (for example oligogalacturonate lyases, proteases, pectate lyases, pectate disaccharide lyases and cellulases), toxins (for example toxin A, toxin RTX, HigB-2 and Colicin V) and other virulence-associated genes (for example effectors HopM1, YopJ, flagella and pathogenicity factors) as active (Additional file 26: Table S25, Additional file 27: Table S26, Additional file 28: Table S27, Additional file 29: Table S28, Additional file 30: Table S29, Additional file 31: Table S30, Additional file 32: Table S31 and Additional file 33: Table S32). The genes of interest to virulence and survival were mainly identified as belonging to B. goodwinii (Additional file 34: Figure S1A) but also to the putative Clostridioides and Carnobacterium. Using our narrowed-down database, the percentage of bacterial reads in the metatranscriptome of AOD symptomatic samples ranged between 1.1 and 3.3%, with B. goodwinii corresponding to 33–82% of those reads. In non-symptomatic samples, the percentage of bacterial reads was only 0.02–0.07%, with no clear difference between different species. For the host, we could identify genes covering similar categories as for the Swissprot analysis (Additional file 34: Figure S1B), with several defence-associated regulators and genes activated, such as mitogen-activated protein kinase, kinase AtM2K9, suggesting activation of the mitogen-activated protein kinase stress response cascade. Furthermore, the defence-associated genes AtNDR/HIN1-like protein 3 and 13 were significantly upregulated suggesting enhanced defence, along with genes for toxin and ROS resistance such as AtGSTU8 and AtDTX29.
Metaproteome analysis suggests B. goodwinii is the dominant pathogen in the AOD lesion microbiome
In order to further confirm the results of the metagenomic and metatranscriptomic data, by identifying active proteins of the microbiome members of interest to the AOD disease aetiology, and the host defence, we performed metaproteomic analysis using mass spectrometry. Using genes identified via metagenomics and metatranscriptomics as a reference database, we identified 629 proteins, of which 59 were bacterial and 570 eukaryotic (Additional file 35: Table S33). As expected, bacterial proteins were detected in all tissue samples but were generally higher in abundance in symptomatic samples. However, several core proteins of phytopathogenic bacteria were detected only in symptomatic tree samples, such as FliC, PepB, PotF and OmpA, along with various plant-associated transketolases, alcohol dehydrogenases and aldolases (Fig. 3a). In a comparative analysis of symptomatic versus non-symptomatic samples, bacterial proteins PelA (plant cell wall-degrading enzyme) and OsmC (survival) were significantly increased in symptomatic tissue (Fig. 3b), along with host-associated proteins CHI1 (chitinase), CAT2 (catalase), LRX4 (cell wall formation regulation), HIR1 (hypersensitive reaction) and GRP1 and TLP1 (lignin formation). Metaproteomics clearly demonstrates an increased abundance of bacterial proteins, and a host defence response, in AOD lesions.
Using the narrowed-down database described for the metatranscriptomics, we re-performed the metaproteome analysis (Additional file 36: Table S34). We found a total of 122 proteins significantly (p < 0.05, fold change > 2) differently abundant in symptomatic samples compared to non-symptomatic samples using student’s T test, or undetected in all non-symptomatic samples and detected in all symptomatic samples (Additional file 37: Table S35). Additionally, a Benjamini-Hochberg multiple test correction was performed on the p values, which yielded significance (FDR < 0.05) for one gene, BG_04787, a 4-hydroxy-tetrahydrodipicolinate synthase. However, the identified proteins in general supported the metagenomic and metatranscriptomic results. Among the proteins with increased abundance in AOD symptomatic samples were those of B. goodwinii and G. quercinecans, as well as that of our predicted Clostridioides. Only one protein, predicted from oak transcript DN950644.1, homologous to a vesicle transport protein found in Arabidopsis thaliana, was increased in abundance in healthy samples; the other 121 proteins were more abundant in symptomatic samples, of which 23.5% were bacterial. Of those proteins, 18 proteins belonging to B. goodwinii were detected, including a superoxide dismutase, aldolases, a catalase and lipoproteins, suggesting survival in the oak host. Several host pectinesterases, proteases and chitinases were increased in abundance, suggesting PCWDE maceration taking place and anti-fungal activity. Furthermore, 3 proteins belonging to G. quercinecans were significantly increased in symptomatic tissue, including a Colicin V secretion protein, along with a signal recognition protein and a hypothetical protein of unknown functions, suggesting bacterial competition. The data of the narrowed-down database here provides a similar but re-focused picture of bacterial survival and virulence from the most prominent AOD microbiome member, B. goodwinii, along with host defence in the AOD lesions (Additional file 38: Figure S2).