Trees, fungi and bacteria: tripartite metatranscriptomics of a root microbiome responding to soil contamination

Direct detoxification/stress responses

The expectations for plant gene expression within a high hydrocarbon-contaminated environment would include an increase in cytochrome P450 monooxygenases. Out of 26 DE transcripts identified as encoding cytochrome P450 family proteins, 18 were in higher abundance within control roots (cultivated in non-contaminated soil) (Additional file 5). Only two were annotated as putative monooxygenases (c585325_g1_i1 and c601406_g5_i2), with one upregulated and one downregulated with respect to contamination. Another expected set of stress-induced detoxification equipment, known to be transcript abundance dependent, are glutathione S-transferases (GSTs), which catalyse glutathione conjugation to a broad range of xenobiotics in order to facilitate vacuolar loading (compartmentalisation of cellular pollutants). Only a single GST was downregulated in roots from contaminated soil, with 13 upregulated (8 tau class, including GSTU19, 37, 48 and 50) at consistently high relative abundance levels (tpm) (Fig. 2 and Additional file 5).

UDP-glucose:flavonoid 7-O-glucosyltransferase was the most abundant upregulated Salix gene in contaminated roots (c600230_g5_i1; 216.94 tpm), and two other putative isoforms were also highly upregulated (c600230_g5_i2, 90.71 tpm; c600230_g6_i1, 52.53 tpm; Fig. 2).

Such increases in flavonoid glucosyltransferase expression represent a common conjugation mechanism in plant detoxification metabolism [61]. Alternatively, UDP-glucose:flavonoid 7-O-glucosyltransferase, which allows symbiotic interaction with microorganisms in Glycine max [62], can be limiting for vacuolar loading of the isoflavone conjugates, and the corresponding conjugate-hydrolysing β-glucosidase is localised to the plant root apoplast, intriguing in terms of the extensive interaction with microbes evidently ongoing within the roots.

Apart from the UDP-glucose:flavonoid 7-O-glucosyltransferase and glutathione transferases, the most highly expressed DE transcripts in contaminated roots were aquaporins and dehydrin. This could be expected given petroleum hydrocarbon contamination is likely to reduce access to water and thus induce osmotic stress in roots. Three PIP1.1 transcripts (c483320_g1_i2, 125.11 tpm; c589371_g5_i4, 88.10 tpm; c589371_g5_i8, 30.63 tpm), three PIP2.1 (c599604_g6_i1, 79.73 tpm; c576656_g7_i3, 70.73 tpm; c599604_g8_i6, 4.23) and PIP2.7 (c602203_g2_i9, 21.70 tpm) were all highly abundant in contaminated conditions. Beyond a direct plant response to soil conditions, it has long been known that root-associated fungi, particularly ECM fungi, can improve not only nutrient uptake in trees, but also water [63–65]. It is interesting to note, in light of subsequent fungal gene expression, that extreme PIP upregulation has previously been observed as induced by the ECM fungi Amanita muscaria in fine ectomycorrhized roots of poplar (a close relative of Salix, sharing macrosynteny) [66] as well as in Salix by the AMF Rhizophagus irregularis [22].

Nitrogen and carbohydrate transport

Five amino acid transporters were DE, two in higher abundance in non-contaminated trees and three in contaminated trees. Two were distinct cationic amino acid transporters, one in higher abundance in non-contaminated conditions (c598534_g2_i3, 6.65 tpm) and one higher in contaminated conditions (c584679_g1, 18.25 tpm). Another of the transcripts higher in contaminated conditions (c594518_g1_i1) was most similar to AAP6 in poplar (87% identity, blastn e-value = 7e^-29), and the two remaining transporters (one in higher abundance in each condition) were both similar to vacuolar amino acid transporter 1. AAP6, an acidic and neutral amino acid transporter known to be expressed in roots, is distinctive in having high substrate affinity (so could potentially be relevant to low amino acid concentration uptake) as well as having an affinity for aspartate [67]. Seven nitrate/peptide transporters where DE with six in higher abundance in contaminated conditions, sharing close homology with poplar NPF (NTR/PTR family) 5.4, 3.1 and 5.2 proteins. Alongside PIP expression, these are hallmarks of root symbiosis [68] in terms of differential expression and are thought to represent oligopeptide import of AM- and ECM-packaged nitrogen, the principal fungi to plant facet of resource exchange.

In terms of potential concomitant plant to fungi exchange of resources, differential expression of a putative sucrose transporter (c582330_g6_i4) with high abundance in roots from contaminated soil was identified, which would be expected if increased ECM interaction were underway in these conditions. Three distinct trehalose-6-phosphate synthase transcripts were DE (c598745_g5_i1, c536396_g3_i1, c599008_g3_i2), all in higher abundance in roots cultivated in contaminate soil. While the role of trehalose synthesis in plants is somewhat obscure [69], trehalose is the dominant storage carbohydrate in AMF hyphae and implicated in ECM abiotic stress tolerance [70, 71], so expression here is intriguing in light of concomitant fungal DE. Additionally, the sugar transporter SWEET2 (c555872_g1_i5) was in higher abundance in roots from contaminated soil. SWEETs (Sugars Will Eventually Be Exported Transporter) have only been identified relatively recently [72] but have already been strongly implicated in plant-fungal interactions [73], including Solanum tuberosum roots colonised by AMF [74]. In Arabidopsis, SWEET2 has been shown to accumulate in root hairs, cap and epidermis, and in cells in close contact with the rhizosphere [75]. Interestingly, the authors suggest that SWEET2 functions as a bidirectional glucose transporter that could control glucose secretion through limitation/vacuolar loading as part of highly complex and co-ordinated interactions with rhizospheric microbes.

Community association

The most abundant downregulated contigs in Salix roots from non-contaminated soil encoded two distinct organ-specific S2 proteins (Fig. 2).

While these proteins are of unknown function, they have been found to be downregulated in a Medicago truncatula 1-deoxy-d-xylulose 5-phosphate synthase 2 (catalysing the beginning of the MEP pathway) knockdown, which also reduced AMF colonisation of roots [76]. An acidic endochitinase (EC 3.2.1.14, c583738_g7_i2) was upregulated in contamination-treated roots. While this would intuitively be associated with plant defence against pathogenic fungi, a number of studies have found that this endochitinase can encourage symbiotic association with ECM fungi, such as Hebeloma, Suillus, Wilcoxina, Pisolithus, Paxillus and Amanita with spruce, birch and Eucalyptus; the suggested mode of action being that non-fungi-damaging apoplastic chitinase modifies ectomycorrhizal elicitors to facilitate symbiotic interaction [77–80]. When considered alongside the preceding data, Salix root DE reflects plant recognition of foreign molecular patterns [81], the non-recognition impact of effectors [82] as well as alterations in nutrient exchange and water availability known to be drastically altered by fungal interactions [65].

Division by taxonomy

Fungi represented the kingdom with the most diverse genetic response to contamination within the root samples, being the annotation source of the 8252 distinct DE contigs. The paradigm of direct up vs down-regulation collapses when multiple organisms and natural biological complexity are acknowledged in sequencing studies (when the dynamic nature of the metagenome is considered [2]), as the presence of an organism can vary as well as gene expression within organisms. Therefore, DE contigs can represent direct or indirect responses to treatment within a given organism of stable presence (with respect to treatment), including responses to highly diverse, changeable and potentially hostile biological environment, but also represent baseline transcription and metabolism of newly present, absent, growing or diminishing organisms within the system. The contigs annotated as fungal had very distinct patterning in relation to contamination; of the contigs present in both conditions, 88.8% (6184 contigs) had higher abundance in non-contaminated trees while, in contrast to this, of the contigs only present in trees of one condition, 96.3% (1239 contigs) were in contaminated conditions (Figs. 2b, c; 4a and 5). These expression patterns, even before investigating gene function, can be interpreted as a general downregulation of constitutive fungal expression due to contamination and the arrival of a distinct, contamination specific, set of expressing genes potentially representing a shift in community makeup. By comparison, only seven Salix DE contigs were expressed in only trees of a single condition.

The two Ascomycota species most represented within DE transcripts were the closely related Pyronema omphalodes (a saprophyte) and Tuber melanosporum (an ECM fungi) (Fig. 4b and Additional file 6). These species both dominated fungal differential expression and had high enough levels of primary annotated contigs (2303 contigs and 1077 contigs) to merit very confident identification of these species, or of close relatives, within the system. Transcript levels were almost all more abundant within roots of non-contaminated trees (94.1 and 98.4% for P. omphalodes and T. melanosporum, respectively). While an alteration in gene expression in response to hydrocarbon conditions is possible, the extremity of polar expression might also suggest a lower abundance for these organisms due to the impact of contamination stress and/or increased competition from more contamination-tolerant life forms. Tuber has been previously identified as associated with Salix [83, 84]; non-constitutively expressed ECM high-abundance markers were DE in Tuber and were indeed in higher abundance in non-contaminated conditions, such as RAS protein (c594475_g3_i1) [85], implying successful Tuber ECM interaction was underway in non-contaminated conditions but was suppressed by hydrocarbon contamination.

Unlike Ascomycota, Basidiomycota expression was dynamic with respect to contamination, as 1639 contigs annotated from three closely related species (from the Agaricoid family’s Hymenogastraceae and Strophariaceae [86]) Heboloma cylindrosporum (an ECM fungi), Galerina marginata (predominantly white rot) and Hypholoma sublateritium (white rot) were downregulated, while 1745 contigs, with a highly distinctive expression pattern, were upregulated but broadly annotated from 61 genera (Fig. 4c and Additional files 6 and 7). These upregulated contigs were distinctive in the extreme diversity of annotation origin and as the majority (71%) were present only in roots from contaminated soil. To investigate the origin of these phytoremediation responsive contigs, and provide additional evidence towards their species of origin, we explored the species diversity of equally good hits using ‘secondary annotation’ [2] (Fig. 4d; Additional file 1) as well as using nucleotide blastn (a strict blastn reduces observation of fungal data by 99.5%, Additional file 8).

Secondary annotation

By utilising secondary annotation, the retention of equally (statistically) ‘good’ annotation hits, all 12,576 DE contigs can be more confidently assigned to an organism of origin by elucidating any ambiguity in annotation which is often overlooked. The set of 1745 Basidiomycota genes most often had secondary annotation from known ECM (such as Scleroderma citrinum and Paxillus involutus) and saprophytes (such as Pleurotus ostreatus and Trametes versicolor), in agreement with the current awareness that the continuum of mutualism is complex [87]. Very little crossover was observed between Basidomycota and Ascomycota in terms of differential expression due to contamination response, as Basidomycota extensively dominated the increased expression in response to contamination conditions (Fig. 4a–d). Secondary annotation provides an additional benefit in terms of mining useful functional description from across the world’s major data repositories (Additional files 2, 3, 4, 5, 6, 7, 8 and 9), but it is also interesting to note that retention of secondary annotation independently reflects classical Basidiomycota phylogeny ([86]; Fig. 4d) (albeit with a crude thresholds limiting data retention to only the very closest homologues), particularly interesting is if this idea is considered in the context of the 189,839 total annotated contigs that include 491,505 secondary annotation (Additional file 4).

An important take home message here, when querying unknown sequences using BLAST, is the benefit of not taking the single top returns as fact. Pertsemlidis and Fondon [88] have detailed how differentiating proteins of high homology from BLAST scoring system should be performed with caution and others have outlined the extensive risk of confounding biological interpretation by not acknowledging the uncertainty of annotation [2]. Secondary annotation of Basidiomycota was spread widely across Agaricomycotina for upregulated contigs whereas downregulated contigs (> 95%) derived from a specific Agaricoid clade comprising the three closely related genera (Hypholoma, Galerina and Hebeloma) (Fig. 4c) which, understandably, shared the majority of secondary annotation with each other (Fig. 4d). For the distinct, upregulated Basidiomycota, secondary annotation revealed no dominant source of annotation (homology) to our contigs within Agaricomycotina (more detailed sequence origin investigation is presented in Additional files 1 and 8).

Upregulated Basidiomycota function

Given that the intricacy, or perhaps futility, in differentiating between saprophytes and ECM fungi has been well discussed [89, 90], determining the fungal mode of action from expression study is problematic [91] and is further confounded by the generally accepted belief of multiple ECM evolutionary events in Basidiomycota ([92]). The scale of mycorrhizospheric complexity is very high; once DE genes have been sub-selected as best annotated from fungi (keeping in mind that 2575 DE contigs were unknown, having no known sequence similarity), further sub-selected based on Basidiomycota annotation and even further sub-selected based on response to contamination (just those contigs in higher abundance in contaminated conditions), 1745 DE-annotated transcripts remain to describe potential functionality (driven by positive expression) within the system. Of these 1745 DE transcripts, 70% (1227 contigs) had unknown function and were annotated as hypothetical, predicted, putative or uncharacterised proteins from Basidiomycota species. These poor functional description terms are selected against during annotation if an equally good hit is available in secondary annotation, in practice re-mining any confident functional descriptions available within the major protein data repositories.

In terms of recognisable gene function, the most abundant fungal contig expressed in contaminated conditions was the relatively cryptic plasma membrane proteolipid 3 (Pmp3), best annotated from Moniliophthora, with 301.42 tpm (c553133_g5_i1, 2128.20 FC) (Additional file 7). The small hydrophobic Pmp3 (dissimilar to recognised hydrophobin structure) has been characterised as highly conserved in fungi, is environmental stress induced (cryptic in yeast under standard laboratory conditions), involved in cytotoxic cation tolerance [93] and sphingolipid synthesis [94] so could be related to membrane integrity maintenance given the challenge of contamination conditions. Generally, the most highly expressed contigs within the group of distinctive Basidiomycota (Figs. 4c and 5) were not putatively related to hydrocarbon degradation but, instead, were related to classical association with plant roots and/or the highly upregulated bacteria (Enterobacteriaceae sp.), namely carbohydrate import, nitrogen (nitrate, ammonia and amino acid) metabolism and transport, and bacterial interaction. This could be expected in light of the hydrocarbon degradation mechanisms DE in the Enterobacteriaceae sp. (outlined below) if the interaction did indeed involve some degree of tripartite mutualism.

Putative hydrocarbon degradation

Surprisingly none of the expected hydrocarbon degrading monooxygenases were identified within the 1745 DE Basidiomycota contigs which increased in abundance due to contamination (Fig. 5, Additional file 7). Only three dioxygenase encoding contigs were identified as DE, all of which had relatively low abundance levels. Two were poorly characterised dioxygenase family proteins (c596412_g1_i1, 8.78 tpm and c593762_g1_i1, 5.93 tpm), but one was a putative extradiol aromatic ring-opening dioxygenase (c598058_g2_i2, 4.23 tpm), functionality well recognised in bacterial PAH degradation studies [95–97] but less familiar in fungi. No lignin peroxidases, manganase-dependent peroxidases (class II peroxidases) or laccases (with potential degradation functionality [26]) were identified as DE, expected if a ligninolytic (white rot) mode of action was underway, particularly as the associated enzyme secretion is highly dependent on transcript levels [91]. Two DE contigs, however, did encode glutathione peroxidase-like proteins (c591991_g1_i5, 8.78 tpm; c591991_g1_i10, 2.14 tpm), and one encoded a thioredoxin-dependent peroxidase (c565655_g1_i1, 48.94 tpm) were upregulated, although this is a common response to contamination-induced oxidative stress [98].

Carbohydrate transport and CAZy

One of most intuitive approaches to distinguish between saprophitic and ECM ecological strategies (the most likely in Basidiomycota here) would be to compare the expression of genes for plant cell wall degradation and carbohydrate import mechanisms, although this is a non-trivial task [91]. For instance, the molecular mechanisms underpinning much of the carbon transfer to ECM fungi from plants is unclear; while there is strong evidence that up to 30% of the plant’s total photoassimilates can be transferred to ECM fungi [99], few hexose transporters have been experimentally validated. The monosaccharide transporter MST1 in both Amanita muscaria (fourfold upregulation during symbiosis [100]) and Laccaria bicolor [101] being the exception to this.

Basidomycota carbohydrate active enzymes were identified as upregulated in response to contamination conditions, with the majority only expressed in contaminated trees (only 9 contigs were also expressed in control trees). Of these, 20 contigs belonged to one of nine glycosyl hydrolase (GH) families, eight contigs to one of six glycosyl transferase (GT) families and one contig was a pectin/pectate lyase (PL; pectate lyase 1). The most highly represented group was GH13, with seven contigs (5 distinct enzymes including 3 likely isoforms) including a secreted alpha-amylase. It is possible that these may represent extracellular invertase activity (sucrose hydrolysis) as opposed to native fungal α-glucan metabolism.

The most abundant CAZy contig (104.28 tpm) was a GH45 Barwin-like endoglucanase (expansin family protein) (Table 1). Expansin domain containing proteins have previously been recognised in Laccaria bicolor as expressed only in ECM tissue [106], which is fascinating in light of the functional similarity to cell wall loosening in plants [107] and if considered in the context of non-necrosis inducing Salix root tissue remodelling and potential interaction with middle lamella and primary cell walls. In agreement with this, a carbohydrate esterase family 12 protein (putative rhamnogalacturonan acetylesterase, pectin-related) [108] annotated only from the orchid symbiote Tulasnella calospora was DE. Further to this, a pectin/pectate lyase was also DE, recently identified as one of the few cell wall-degrading enzymes retained by the ECM fungus Tuber melanosporum for cell wall remodelling of Corylus avellana as part of their symbiotic interaction [109]. The second most abundant CAZy was GH131 (best annotated from Hebeloma but with equally good hits in Plicaturopsis, Laccaria, Jaapia, Tulasnella and Gelatoporia), including a cellulose binding module (CBM1) with broad β-glucanase activity. While this activity allows for the potential for saprophytic action, remodelling of the plant cell wall and progression through lamella is in agreement with ECM colonisation of plant roots towards the effective formation of the hartig net [110] and, interestingly, GH131 has previously been identified as expressed in Jaapia argillacea as a potential adaption from saprotrophic ancestors [92]. Three GH5 contigs (exo-beta-1,3-glucanase) were also highly expressed including an ectomycorrhiza-upregulated exo-beta-1,3-glucanase [106, 111]. In disagreement with this is the lack of identifiable lignases or maganase-dependent peroxidases necessary for deconstruction of heavily lignified middle lamella as well as the differential expression of galactose oxidase (AA5), a secreted extracellular catalyse which has twice as much activity on galactomannan (mannan backbone with galactose sidechains) compared to galactose [112], important given this very common hemicellulose component is more indicative of plant secondary cell walls [113].

In terms of the potential direct interface between the plant and putative ECM Basidiomycota, two Stereum hirsutum ricin B-like lectins of the carbohydrate-binding module family 13 (CMB13) contigs were very highly abundant in treated conditions at 144.48 tpm (14th most abundant, 1350 FC) and 98.05 tpm (22nd most abundant, 350 FC). Interchangeable homologous hits were also found using secondary annotation within Heterobasidion, Fomitiporia and Jaapia. Lectin carbohydrate-binding is well documented as one of the means of intense binding at the interface between the ECM hartig net and the plant host cell wall [114, 115]. CMB13 has been shown to have specificity for backbone xylan (such as the majority of hemicellulose within willow root 2° root cell wall), although is also found in a diverse array of non-xylanase glycosyl hydrolases [116].

The metabolic machinery necessary as a downstream consequence of monosaccharide import (in light of very highly expressed MST) includes the Embden-Meyerhof (EM) pathway, pentose phosphate pathway and tricarboxylic acid cycle (Additional file 7). Trehalose and glycogen metabolism within the EM pathway, clearly represented within DE contigs, are glycolytic pathways expected in ECM fungi [117]. Alpha,alpha-trehalose-phosphate synthase (TPS1), two trehalose-6-phosphate phosphatases (TPS2), two glycogen synthases, glycogen phosphorylase, glucokinase regulator, hexokinase, phosphoglucomutase, glucose-6-phosphate isomerase, two phosphoglycerate kinases, two phosphoglycerate mutase-like proteins, enolase, pyruvate kinase and two fructose-1,6-bisphosphatases were identified as DE and in higher abundance in root from contaminated soil. Within the pentose phosphate pathway: glucose-6-phosphate 1-dehydrogenase, 6-phosphogluconolactonase, two 6-phosphogluconate dehydrogenases, transketolase and transaldolase were also DE. While carbon was fated towards ethanol fermentation (pyruvate carboxylase was DE), the citric acid cycle was also comprehensively represented as upregulated: pyruvate dehydrogenase (E1 component subunit alpha; acetyltransferase component), four phosphoenolpyruvate carboxykinases, pyruvate carboxylase, three succinate dehydrogenases, two fumarate reductases, succinyl-coa synthetase (alpha chain and beta chain), three malate dehydrogenases, 2-oxoglutarate dehydrogenase, two homocitrate synthases, two aconitate hydratases and isocitrate dehydrogenase, supporting carbon scaffold assembly towards amino acid production.

Nitrogen management

A broad suite of proteases was also DE including four endopeptidases (including subtilisin and bleomycin), two Zn-dependent exopeptidases, three carboxypeptidases (glutamate, serine and zinc), four metallopeptidases, a methionine aminopeptidase, two putative aminopeptidase isoforms and two distinct aspartic peptidase A1 (Table 2). Secreted suites of proteases are thought to be similar between Basidiomycota saprotrophs and EMC [90]. However, the extracellular protein degradation pathways identified here distinctly match those expressed in the EMC P. involutus [120]. One of the markers identified from Laccaria biocolor but recognised across ECM are Ras/Ras-like proteins [85, 122, 123], small GTPases involved in cargo sorting in coated vesicules [124] (vesicular turnover is thought to be high in ECM due nutrient and signal exchange with hosts [125]) that have been shown to increase in transcript abundance after successful establishment of symbiosis [126]. Eleven Ras-like proteins were DE and in high abundance in the phytoremediation responsive Basidiomycota, importantly including likely Ras, or Ras interacting, proteins known as host interaction-specific ECM markers identified in Laccaria [85]. These comprised three highly abundant uncharacterised small GTPase-binding proteins (47.83 tpm, 31.86 tpm, 9.42 tpm), Rho1 (20.67 tpm), Ras-related protein Rab-5B, Ras protein, Rab-type small GTPase, Ras GTPase-activating protein, Rab GDP-dissociation inhibitor, GTPase foz1, Sar1-like protein member of Ras-family as well as an extremely highly expressed Rab5 (ypt5) GTPase family protein (57.28 tpm), involved in endocytotic vascular trafficking [127].

While a clear depth of functional detail can be revealed by metatranscriptomics, the majority of putative fungal proteins were uncharacterised (having no functionally characterised homologue), a useful reminder that the vast majority of the natural world remains to be explored.

Polyadenylation in bacteria

Polyadenylation of RNA has long been known to occur in bacteria [128, 129]. Some studies have established the potential involvement of polyadenylation in mRNA degradation [130, 131] and RNA quality control [132, 133] in a limited number of bacterial species, although the contemporary picture of polyadenylation in bacteria suggests complexity beyond this, as recently described by Kushner [134]. Using pulse-labelling, levels of polyadenylated RNA have been measured at up to 15% of RNA in Escherichia coli [135] (in contrast to those of up to 25% in gram-positive Bacillus brevis [136]) and more recent research using microarray analysis revealed, in an E. coli K-12 wild-type transcriptome, that 90% of transcribed ORFs underwent some degree of polyadenlyation [137].

Within the plant transcriptomic research community, bacterial mRNA is routinely discarded during early quality control of common bioinformatics pipelines as distinct from the target organism of study or discarded intrinsically when mapping to a reference genome. Gonzalez et al. [46] recently reported how the sub-selection of only transcriptomic sequences expected a priori can confound biological results, the example leading to S. purpurea biotic stress genes being misidentified as abiotic stress responsive genes due to a strong treatment-specific interaction of a plant herbivore (Tetranycus urticae). Similarly, Brereton et al. [2] demonstrated the potential for mistaken mapping (mis-mapping) RNA-seq reads from unexpected foreign organisms to technically confound results by mapping known foreign sequences to the S. purpurea genome. Given these potential pitfalls, it would seem prudent to acknowledge the complexity of extra-laboratory biological systems involving higher eukaryotes by investigating all mRNA molecules present within any sample even when strong expectations exist within the experimental design. The approach is limited here as an unknown absolute proportion (and community) of bacterial mRNA could be lost during polyA enrichment.

Despite polyA enrichment, an extremely broad diversity of bacteria was observed within the system (including contigs that were not DE; Fig. 6a)). The majority of these contigs derived from Proteobacteria, including Alphaproteobacteria (14%), Betaproteobacteria (34%), Deltaproteobacteria (13%) and Gammaproteobacteria (23%). Actinobacteria (14%), Bacteroidetes (4%) and Firmicutes (4%) were also highly represented. This community makeup is similar to that previously reported in metagenomic studies of contaminated soils [25].

A different picture was discerned in DE contigs (Fig. 6b); of 638 DE contigs, 86% were in higher abundance in contaminated trees. Enterobacteriaceae species were the most represented, annotating 72% of all the bacteria contigs. Importantly, 100% of these Enterobacteriaceae contigs were of higher abundance in contaminated trees (Fig. 2b, c), suggesting strong biological association with contamination and/or the other organisms. Therefore, their functionality was explored, albeit with caution. Of potential relevance to this, growth under more challenging conditions than standard rich media, or repression of growth using transformation, can increase the level of polyadenylated transcripts within E. coli [138–140].

Decoding polycistronic transcriptional units

Determining transcriptional units (putative operons) in bacteria is computationally problematic in terms of prediction from genomic sequence [141, 142]. While prediction of E. coli K12 operon structure is as advanced as in any bacteria, little confidence can be established in the transferability of such predictions given operon variation across organisms and very high responsiveness to environmental change. The unique nature of any extra-laboratory environment, such as that explored here, as well as the high likelihood of numerous novel organisms present (not the least being the bacteria themselves present within the system) further complicates transcriptional unit prediction. That being said, top-down de novo assembly provides very high confidence in the presence of observed sequences within the system, of particular technical value is the high requirement for equal mapping coverage applied within Trinity’s contig assembly [3]. Of the 639 DE and potentially polycistronic contigs annotated as bacterial in origin, 3134 protein coding regions were identified using Trinity’s transdecoder. In total, 489 of these contigs were identified as deriving from species within the family Enterobacteriaceae (100% in higher abundance in contaminated trees) comprised 2834 proteins, although a substantial proportion were duplicates with only minor sequence variation (Additional file 10). These could represent common expression by distinct organisms or multiple operons from a single organism (being somewhat reminiscent of eukaryotic splice variants at a sequence level). The majority of the 2834 Enterobacteriaceae putative proteins were annotated from E. coli (68%, across a broad spectrum of strains, Additional file 10) or the genus Shigella (29%). A recent identification and genome sequence of a poplar growth promoting endophyte (Enterobacter sp. 638) by Taghavi et al. [31] shared similar levels of homology as those found here but at a genome level, with 69% of predicted CDS being similar to E. coli.

Upregulated bacterial gene function

A broad spectrum of contigs encoding proteins putatively involved in survival within a rhizospheric environment were identified as DE due to contamination conditions (Additional file 10), including those putatively interacting extracellularly and/or intracellularly with ECM hyphae and/or plant roots. As the vast majority of DE bacterial genes were best annotated as Enterobacteriaceae (either E. coli or Shilgella strains), E. coli K-12 MG1655 gene nomenclature was used where possible.

Hydrocarbon degradation and biosurfactant production

For PAH degradation, bacterial dioxygenase enzymes are considered the principal means of ring fission [26, 143] (although cytochrome p450 monooxygenase are common), whereas alpha-ketoglutarate-dependent dioxygenases (alkB) are one of the most recognised enzymes driving degradation of aliphatic hydrocarbons [25] (not to be confused with alkane monooxygenase, alkB). The contig c596278_g1_i5 was upregulated here (containing yojIalkBadaapbEmqo), functionally comprising an ABC transporter ATP-binding protein (yojI), alpha-ketoglutarate-dependent dioxygenase (alkB), regulatory protein (ada), thiamine biosynthesis lipoprotein (apbE) and malate:quinone oxidoreductase (mqo) (Fig. 7 and Additional file 10). Additionally to this, a number of genes well characterised as accompanying toluene tolerance [144, 145] were DE, including a putative operon which contained the toluene transporter subunits mlaCDEF and kdsCDlptACkdsCD encoding a lipopolysaccharide export system. Degradation of alkanesulphonates, previously identified as involved in crude oil degradation by metagenomic study [25], occurs through desulfonation and relies on expression of the operon ssuABCDE [146, 147] comprising an alkanesulfonate transporter subunit, a putative alkanesulfonate transporter subunit, an FMNH(2)-dependent alkanesulfonate monooxygenase, putative aliphatic sulfonate binding protein and NAD(P)H-dependent FMN reductase. Here, the entirety of the ssuABCDE operon was assembled and DE in higher abundance in contaminated conditions (two similar contigs c595976_g1_i2 and c595976_g1_i9). As would be expected alongside the ssu system, the associated tau system (comprising tauABCD genes) was also upregulated within two separate putative operons containing tauAB proteins (taurine transporter subunits A and B; c596422_g1_i2) and tauCD (taurine transporter subunit C and alpha-ketoglutarate-dependent taurine dioxygenase; c596422_g2_i1). These oxidising mechanisms, such as the well-known alkanesulfonate monooxygenase ssuD [148], suggest the bacterial species present in association with willow roots is actively expressing an enzyme suite capable of degrading the hydrocarbons present within contaminated soil and so may play an important role within any tripartite interaction of mutual benefit.

Genes necessary for three biosurfactants were found to be most represented in metagenomic study of crude oil contaminated soil [25], (in descending order of abundance): trehalose lipids [149], polyol lipids and mono/di-rhamnolipids [150]. Interestingly, the principal enzymes involved in the trehalose degradation pathway were DE, with two putative operons, the first containing treAdhaLKMR (c602247_g2_i1), from E. coli and the second (c567422_g1_i2) containing treF (from shigella) alongside a transcriptional regulator yhjB (from E. coli). These potential operons comprise the periplasmic trehalase (treA) and PTS proteins (dhaLKM, although no transport activity is expected [151]) as well as the cytoplasmic trehalase (treF) expected for trehalase utilisation [152]. While this suggests trehalose may not be being employed as a biosurfactant, it is compelling in terms of the potential interaction with ECM Basidiomycota as the specific contamination-responsive DE contigs include trehalose biosynthesis genes. The entire rmlABCD (often termed rfb) operon (c433320_g1_i1) was in higher abundance in contaminated conditions, comprising the l-rhamnose synthesis genes essential for rhamnolipid production in E. coli [150, 153]. Additionally to this, rhamnulose-1-phosphate aldolase (rhaD), rhamnulokinase l (rhaB) and rhamnose-isomerase (rhaA) were all upregulated allowing provision of l-rhamnose from dihydroxyacetone-phosphate [154], present within two operons comprising rhaDyiiL with genes for PTS system IIA and IIB components (c599913_g4_i1, operon structure common in Firmicutes) and rhaABSRyifCLMaidB (c433320_g1_i1). While the l-rhamnose synthesis pathway is considered well characterised, the exact 3-hydroxy fatty acid precursors of rhamnolipids are less confidently known, although it has been hypothesised that both the type II fatty acid synthesis pathway [155] and β-oxidation [156] can provide lipid precursors. An operon was upregulated (c601578_g3_i4) containing numerous members of the type II fatty acid synthesis pathway (FAS-II): fabG (syn: beta-ketoacyl-ACP reductase, similar to rhlG of rhamnolipid synthesis in Pseudomonas aeruginosa [157]) fabDFHychHyceDFGthiK, the operon structure being unsurprising in E. coli. rhlAB, recognisable as driving rhamnolipid production in Pseudomonas aeruginosa [155], was not detected however, and an operon containing fabA, potentially a competitor for β-hydroxydecanoyl-ACP [158], was DE. Of β-oxidation, two fadIJsixA upregulated operons (c591107_g2_i1 and c591107_g2_i2), a common operon if E. coli K-12, were in higher abundance under contaminated conditions. fadIJ are recently discovered homologues of β-oxidation fadAB in E. coli with a suggested preference for short and medium chain fatty acid degradation [159].

Community interaction

These DE potential operons encompass a wide array of genes considered cryptic in E. coli (silent under standard laboratory conditions), for example, cryptic genes in E. coli K-12 substr. MG1655 [160] that were expressed here included ecpABCD (c589067_g3_i3), ebgAC (c601454_g4_i1), gspAB c602128_g2_i1), gspCD (c600465_g1_i1) and chiAgspO (c600465_g2_i2) (Fig. 7 and Additional file 10). The expression of genes which are silent in many E. coli species cultured under laboratory conditions is perhaps unsurprising given the mycorrhizal environment; this is further supported by DE of genes representing fungiphile metabolism [161] in the Enterobacteriaceae species responding to contamination conditions. The most convincing in this regard were genes involved in interaction with chitin (fungal cell wall) and the well-known fungal exudate oxalate, used for habitat manipulation through pH lowering and increasing availability of nutrients [162] as well as for lignocellulosic degradation in many Basidiomycota. The putative operon c591446_g1_i5 comprised the well-studied two-component (stimulus-response) regulatory system evgSA [163] with an Oxalyl-CoA decarboxylase, acetyl-CoA:oxalate CoA-transferase and transporter similar to the metabolic machinery essential for E. coli tolerance of oxalate [164]. The putative operon c599707_g2_i2 comprised 8 proteins including chiPQ (oprD family, a chitoporin for uptake of chitosugars and chitosugar-induced lipoprotein which is not constitutively expressed [165, 166]) and a phosphoglymutase, suggesting tight interaction with the chitin phosphotransferase system (PTS) system [167]. In line with this is the upregulation of a potential operon containing four proteins including chbBC (N,N′-diacetylchitobiose-specific enzyme IIB component and permease IIC component of the PTS system) alongside osmE and nadE (an uncharacterised osmotically induced lipoprotein and NH(3)-dependent NAD(+) synthase). Further supporting this, the operon chbACFGR (c599728_g5_i1) comprising the more classically recognisable chitobiose operon [168], with the exception that chbB is absent (expressed in the previous operon alongside a distinct chbC). The cryptic chitinase chiA (c600465_g2_i2) was also expressed in an operon including bfr and gspO, pertinent given the ECM fungal interaction suggested by gene expression (potentially explaining cryptic expression as the fungal environment is difficult to replicate in culture) and logical in light of associated [169] DE of (non-constitutively expressed) type II general secretion systems (gspAB, c602128_g2_i1 and gspCD, c600465_g1_i1 functionality is discussed below).

While the clear expression of chitin degrading mechanisms could suggest the interaction with ECM fungi is endocellular biotrophy or necrotrophy, the possibility of extracellular biotrophy cannot be discounted as any complex interaction could involve modification of, or macro-(biofilm)-adhesion with, the Hartig net involving chitinases. In line with this, Chittero et al. [170] report biotrophic interaction of Pseudomonas fluorescens and Bacillaceae through the use of chitinases with fruit bodies of the ECM Tuber borchii; while others found bacterial chitinase production did not inhibit ECM fungi [171, 172] and that when comparing rhizospheric and ectomycorrhizosphere bacteria, those capable of hydrolysing chitin were present only within the ectomycorrhizosphere [173].

Biofilm formation

There are five stages to biofilm formation which are generally accepted: (i) surface contact and reversible attachment, (ii) irreversible attachment, (iii) microcolony formation and early development of architecture, (iv) maturation and (v) dispersal. The accomplished mini-review by Van Houdt and Michiels [174] outlines the proposed surface determinants associated with each stage in E. coli. Flagella have been associated with the initial surface contact and reversible attachment as well as final dispersal and motility (although development and maturation stages are not dependent on flagella). Irreversible attachment coincides with production of type I fimbriae, extracellular polysaccharides and curli fimbriae. Development of biofilm architecture also involves curli and extracellular polysaccharides as well as colonic acid production and the self-recognising adhesion antigen 43 (also called flu). Maturation of the biofilm architecture also involves colonic acid production and curli as well as conjugative pili. Dispersal of bacterial cells back into the environment again involves flagella. Operons encoding proteins from each of these processes were DE, all being in higher abundance in contaminated conditions. The formation and adhesion of bacterial biofilms to fungal hyphae has been well characterised as facilitating the exchange of nutrients [28, 29] but has also been shown to allow the mobilisation of pollutant degrading bacteria through the environment for mutual benefit [175, 176].

As substantial numbers of genes involved in biofilm formation were DE, we compared gene expression here to the list generated by Tenorio et al. [177], who systematically investigated genes affecting biofilm formation in E. coli through overexpression. While the study was performed using LB media, we expect extensive overlap in machinery involved in potential biofilm formation within an ECM environment. Ten genes were identified as fundamentally altering biofilm abundance: ccmf; fdrA (syn yahF); flgFGIL, fliH (c601578_g1_i1); gspA; secY and wcaK (Fig. 7 and Additional file 10). All were DE here with the exception of secY and wcaK; however, other colonic acid production genes were upregulated (wbz wcaJLF) as well as extensive DE of the secYEG translocon (peptide export complex proteins) including, secA, secE and SecDFyajC. A total of 11 of the 35 genes whose overexpression altered biofilm architecture in LB media were DE with increased abundance in roots of contaminated soil, while 11 of the 27 genes causing filmentous morphology of the biofilm were also upregulated. The flagellin synthesis genes FliA and FliC, and flagellar biosynthesis master regulator FlhC, were identified as in high abundance in treated samples as well as fliEFGHIJKLMNOP, fliDST, flgKL and flgGHIL. Flagella can allow increased fitness when growing/forming biofilm on fungal hyphae [178], although E. coli are thought to require flagella only for the initial stages of biofilm formation [177]. Expression of these genes has been shown to be positively regulated by polyadenylation [139]; this is particularly interesting given the suggestion that this mechanism could be specifically in place to allow the bacteria to adapt and survive in challenging environmental conditions (through PAPI regulation).

The chaperone-usher pathway is a delivery system for type I pili through the outer membrane. The periplasmic transport of proteins is essential for successful pili biogenesis and the chaperone-usher pathway utilises the SecYEG translocon to achieve this, also upregulated here in contaminated conditions. Concordantly, the fimDFGH operon (c600702_g6_i3, structurally common in E. coli K-12 strains), encoding the outer membrane usher protein for type 1 fimbrial synthesis and three minor component of F-type fimbriae, was upregulated in contaminated conditions, as was sfmCD (c575401_g1_i2; fimACI was not identified as DE). Interestingly, given the likely interaction with fungal cells, sfm is a chaperone-usher operon that is silent under laboratory conditions but was shown to express when the type 1 fimbriae complex is knocked out and could promote adhesion to eukaryotic cells [160]. Exploring the idea that a substantial amount of biofilm-related expression within this microbiome environment could indeed be cryptic within laboratory environments, Korea et al. [160] identified a wide range of these chaperone-usher fimbrae as associated with distinct surface specialties in E. coli K-12.

The contigs assembled here from root samples and upregulated in contaminated conditions did indeed represent the majority of these previously identified cryptic chaperone-usher fimbriae, including yfc (c601154_g1_i3; yfcRQOU(usher papC)SP(chaperone papD)B(prmB)aroC), sfm (c575401_g1_i2; sfmCD); fim (c600702_g6_i3; fimDFGH), yad (c601794_g1_i1; yadHG, c601794_g2_i3; yadCKLMhtrEecpD and c601794_g3_i1; yadECDI), yde (c591908_g2_i3; ydeTQSR), yqi (c601683_g14_i1; yqiGHI), yeh (c598114_g3_i1; yehABC) and yra (c600252_g5_i1; yraJ) (Fig. 7 and Additional file 10). These represent putative operons, predicted in E. coli K-12, which do not contribute to adhesion of E. coli under normal laboratory conditions but do promote biofilm formation on uncommon abiotic surfaces (‘unknown environmental niches’ [160]) and do contribute to adhesion to eukaryotic surfaces. The ELF operon (elfADCG) [179] was observed within a single contig as the seven member elfADCG-ycbUVF operon (c601320_g2_i5) predicted in E. coli K-12 with an additional protein, pyrD dihydroorotate dehydrogenase gene (quinone). Interestingly, a number of proteins expressed within these contigs are suggestive of potential functionality of association with a host: ydeTQSR included a serine/threonine-protein kinase HipA and the toxin component of a HipA family toxin/antitoxin system, and pantothenate synthetase, 3-methyl-2-oxobutanoate hydroxymethyltransferase and aspartate 1-decarboxylase (both of which participate in pantothenate and coA biosynthesis) were expressed along with yadECDI.

In addition to these extra-laboratory expressed operons associated with biofilm formation, a number of DE genes were associated with swarming in Pseudomonas aeruginosa, a process strongly related to biofilm formation. Four of the six genes explored as swarming negative mutants by Overhage et al. [180] were necessary for good biofilm formation; all were DE here and in higher abundance in contaminated conditions: rhlE (c601160_g1_i1, ATP-dependent RNA helicase), lptA (c602252_g5_i2, lipopolysaccharide transport periplasmic protein), gshB (c586692_g1_i1, glutathione synthetase) and acsA (c600318_g2_i1, acetyl-CoA synthetase). Additionally, two variants of traABCEKLNPUVWYtrbCL (c582511_g1_i6 and c582511_g1_i7) are comprised of conjugative apparatus (pili construction) common to F-like plasmids and necessary for mobility of genetic elements between bacteria (Fig. 7 and Additional file 10). None of the bacterial secretion systems are thought to be constitutively expressed but are instead triggered by recognition of host molecular pattern through adhesins [181, 182]. As discussed with relation to chitinase secretion, gspAB and gspCDO (with chiA) of the type II secretion system were upregulated in contaminated conditions. Given this pattern, it would be expected that xcpU (the general secretion pathway protein H, gspH in E. coli) would be DE; however, it was not identified as such.

The common pilus (ECP) operon is widely conserved throughout E. coli and required for early-stage biofilm development and host cell recognition [183], although usually silent in E. coli K-12 MG1655 under laboratory conditions [184]. As discussed above, an operon containing the regular common pilus machinery, ecpRABCDEykgJyagU (c589067_g3_i3; synonyms include but are not limited to matABCDEFykgJyagU) was in higher abundance under contamination conditions. Interestingly, the yagTSRQ operon (c589067_g1_i4; synonym paoABCD), common in E. coli K-12 strains, was also DE, encoding periplasmic detoxification enzymes (xanthine oxidase family [185]) including yagT, a twin arginine translocator (tat) pathway signal sequence. While the secYEG translocon can transport unfolded proteins, the tat can transport folded proteins to the periplasm (around 6% of all secreted proteins are thought to be tat dependent [186]). These secreted proteins include hydrogenases, dehydrogenase and nitrate reductases. The tat-dependent hydrogenase-2 operon (redox) (hybGFEDCBAO; c600034_g1_i12) was identified as DE, including the Tat pathway signal sequence domain protein hybA [187].

Curli fimbriae are outer membrane adhesins employed by E. coli to facilitate biofilm adhesion [177, 188, 189]. A master regulator of biofilm formation, csgD [189], is at the centre of the curli fimbriae regulatory network. Two variants of the csgD containing operon csgDEFGycdXYZ (c596842_g1_i2 and c596842_g1_i4) were DE and in higher abundance in contaminated roots (Fig. 7 and Additional file 10). csgDEFG encode four proteins necessary for curli assembly: csgG is involved in pore formation in the outer membrane, whereas csgE and csgF are periplasmic proteins which physically interact with csgG [188]. csgD is a positive regulator of csgAB (classically expressed on a distinct operon), the major curli extracellular structural subunit (csgA) and its nucleator (csgB). While csgAB was not identified as DE here, it has been recognised that this is often the case (csgAB not being expressed in cells forming curli) and that assembly can occur via interbacterial complementation, it may also be true here that subunit production is not transcript dependent or is constitutive and constant (and so not DE). Additionally to curli, the prominent surface protein flu was also upregulated (c601850_g2_i6), a self-recognising adhesin known to induce autoaggregation to promote biofilm formation [190].

Nitrogen management

Similar to the recently sequenced plant growth-promoting endophyte Enterobacter sp. 638, which was isolated from poplar stems, it seems likely that no nitrogen fixation was occurring within the Enterobacteriaceae species responding to contamination as no nif genes were identified as DE, with the exception of the pyruvate-flavodoxin oxidoreductase (nifJ/ydbK) (Fig. 7 and Additional file 10). Two variants of the fixABCX operon were also upregulated (c602208_g1_i1 and c602208_g1_i3; intricately linked with N fixation in diazotrophs [191, 192]), but in E. coli, these genes form an important reducing role in carnitine metabolism [193, 194].

Differential expression of well-characterised respiratory nitrate reduction (nar) operons were identified including narK (c598952_g5_i3) used for nitrate/nitrite extrusion in E. coli [195]; narUZ (c598132_g7_i1) where narU is also used for nitrate/nitrite extrusion while narZ is part of nitrate reduction, as well as narXLychONChaC (c598952_g3_i7; similar to 282397 (TU0650)), nitrate/nitrite sensor and regulator proteins. Two slightly divergent versions of the periplasmic nitrate reduction (nap) operon was also DE, napABCGHccmABCFE (c602069_g2_i2 and c602069_g2_i4; similar to 425132 (TU1418)) [196]. The important nitric oxide reduction (dissimilatory reduction) equipment in E. coli, via flavorubredoxin expression encoded by the NorVW operon [197], was also DE. The nitrogen assimilation control regulator nac was also present within an upregulated operon ((c598923_g2_i7; common to in structure to 580938 in E. coli BL21(DE3)) containing the cbl regulator, involved in aliphatic sulfonate utilisation, and yeeO MATE efflux transporter (involved in flavin secretion). In terms of ammonia generation, two versions of anaerobically expressed nitrite transport and reduction nir operons [31, 198] were DE: nirBtsgA (c598930_g2_i1) and nirBCDcysG (c598930_g1_i1). Amino acid (glutamine) synthesis under nitrogen demand has also been shown to be fine-tuned (ammonia/ammonium) by the nitrogen assimilatory protein glnK regulation of amtB (ammonium transporter) [199], present here within an upregulated operon containing the amtB and two multidrug efflux proteins, mdlAB, (c602197_g2_i5; all commonly expressed together). Additionally to this, nrfA, a cytochrome c-552 which catalyses nitrite to ammonia in a formate-dependent manner, was upregulated (c600318_g2_i1) as well as an operon containing yddG, an aromatic amino acid extrusion protein, in a operon (c601567_g3_i3) containing the nitrate-inducible formate dehydrogenase major subunit fdnG.

The conversion of inorganic nitrogen, through nitrate assimilation and reduction, to the biologically more useful nitrite and ammonia (for example, to construct amino acids) provides bacteria with nutrition for survival within the environment but also a potential currency for mutualism within the microbiome, or mycorrhizosphere, system. E. coli only assimilates nitrate under anaerobic conditions [200]; this broad spectrum of assimilatory and dissimilatory nitrate reduction equipment with higher abundance under contaminated conditions is perhaps unsurprising given that facultative anaerobic Enterobacteriaceae can outcompete obligate anaerobes in anoxic, organic compound-rich forest soils [201].

Bacterial functional role

So-called mycorrhizal helper bacteria MHB [30] are thought to promote the symbiotic association between plant roots and mycorrhizae, creating a tripartite relationship [29], such as the Bacillus sphaericus isolate EJP109 promoting Suillus luteus (ECM) growth in association with Pinus sylvestris as well as the Streptomyces sp. nov. 505 and S. anulatus (Beijerinck) Waksman 1003 promoting Amanita muscaria (ECM) growth in association with Picea abies [172, 202]. MHB and PGPR (plant growth-promoting rhizobacteria) strains of Pseudomonas have been visualised attaching to ECM hyphae (suggesting extracellular biotrophic mycophagy) of the genus Laccaria, in a bacterial and fungal strain-specific manner [203]. Differentiating between the potential roles of interacting bacterial species within the biological system is difficult; Leveau and Preston [204] suggested three potential modes of interaction, necrotrophy, extracellular biotrophy and endocellular biotrophy, as well as a spectrum of expected functionality associated to each mode. In general here, expression suggested extracellular biotrophy, however, not only are the other two modes of interaction possible but a more complex continuum of interaction would seem likely given the diverse environment. Before such functionality can be confidently elucidated, in addition to the necessary promotion of cross-disciplinary microbiome research, observation of gene expression without constraint to an organism/s expected within a biological system needs to become a standard requirement for transcriptomic investigations.