Using next generation transcriptome sequencing to predict an ectomycorrhizal metabolome

Identification of expressed and enriched membrane transporters

In the ectomycorrhizal interaction, the individual enzymes and transporters of plant and fungal metabolism remain largely compartmentalized, prevented from directly interacting by the apoplastic space. To exchange nutrients, plant and fungus must have the capacity for importing and exporting nutrients across their cell walls. For the ectomycorrhizal metabolome, nitrogen, phosphorous and carbohydrate compounds are essential elements of nutrient exchange and Gene Ontology (GO) annotations [20] were used as a systematic approach to identify expressed genes in the mycorrhizal transcriptome annotated as members of these transporter classes. GO annotations for predicted proteins coded by best gene models from sequenced genomes were generated by automated BLAST homology to proteins of previously identified function [13]. These GO annotations are corroborated by individual publications on specific transporters for amino acids [21], ammonium [21, 22], sugar [19], and inorganic phosphate [13] which have been identified as encoded in the genome.

The lack of enrichment for ammonium transporter genes suggests that ammonium is not one of the primary mediums of interaction between aspen and L. bicolor in spite of the fact that ammonium salts are present in the growth medium. Amino acid transport, enriched in both aspen and L. bicolor suggests that amino acids are a major metabolite that is shared during ECM interaction. A number of nitrogen containing compound such as ammonia and amino acids have been proposed as candidate compounds for nitrogen exchange in ECM systems [26–28]. The correlation of the present model with amino acid transport is consistent with experimental data for arbuscular fungal and plant systems [29]. The transported amino acids are likely originated from the action of fungal enzymes as the growth medium lacks amino acids. Aspen transcriptome is not enriched for sugar transporters, but aspen does not require them. Sucrose molecules enter the apoplast from aspen roots via diffusion, where it is hydrolyzed into hexose sugars by a sugar invertase produced by aspen [30]. The glucose and fructose products of the invertase reaction are substrates for specific fungal hexose transporters [19]. The observed enrichment of sugar transporters in the L. bicolor transcriptome is consistent with the experimental data demonstrating transfer of carbohydrates across the apoplast interface in ectomycorrhizal roots of aspen [31]. The hexose transporter subset of enriched sugar transporters identified in this analysis overlaps with an expression profiling and characterization analysis of hexose transporters proteins encoded in the genome of L. bicolor[19].

Model of the Mycorrhizal Metabolome

The mycorrhizal metabolome map depicts all the metabolic compounds inferred by expressed genes annotated with enzyme functions (Figure 1). The respective metabolomes of these organisms are connected by a suite of transporters that shuttle essential nutrients across the apoplastic space. Potential gateways for linkage of these compartmentalized L. bicolor and aspen metabolisms were identified using the amino acid and sugar transporters enriched in the set of mycorrhizal expressed genes. The predicted mycorrhizal metabolome is comprised of 3472 metabolic compounds connected by 5156 enzyme reactions. Reactions are mediated by 948 specific enzyme functions, of which 168 enzyme functions are uniquely expressed by aspen and 85 enzyme functions are uniquely expressed by L. bicolor.

Mapping aspen expressed genes to metabolic pathways in KEGG

To better visualize the predicted Mycorrhizal Metabolome, reactions in the complete mycorrhizal transcriptome were mapped onto KEGG Global metabolism (map01100). Global Metabolism map is comprised of 145 metabolic sub-networks [32]. This map is an expert-curated metabolic pathway that highlights an organism's core metabolic pathways and is a useful tool for visualizing key features out of the much more inclusive set of all possible metabolic interactions in the complete mycorrhizal metabolome. The significantly expressed genes of the mycorrhizal transcriptome that code for proteins annotated with EC activities were mapped to the global KEGG Metabolism Map and generated a sub-set of enzymatic reactions predicted to take place in aspen- L. bicolor mycorrhizal system (Figure 2). In this illustration, the nodes correspond to metabolic compounds while the edges represent enzymatic transformations of the compounds. Although many unique enzyme functions are shared in both organisms, expressed enzyme functions unique to the aspen or L. bicolor transcriptome are indicated by green or red edge colors, respectively. Of all enzyme reactions represented in KEGG global metabolism map01100, 30% are expressed in complete mycorrhizal transcriptome, 25.2% in aspen transcriptome, and 21.5% in the L. bicolor transcriptome. Of the 145 KEGG subnetworks that make up global metabolism KEGG map01100, 58.6% in the complete mycorrhizal transcriptome, 31.7% in the L. bicolor transcriptome, and 46.0% in the aspen root transcriptome are detected with at least 30% of enzyme activities significantly expressed.

(An interactive version of Figure 2, allowing a user to zoom in on any of the 145 metabolic sub-networks of map01100, is available in Supplementary Data, Figure S1-a http://www.bio.anl.gov/molecular_and_systems_biology/MycorMetabolome/Supplementalfigures.htm.

Expressed mycorrhizal enzymes include core metabolic activities

The utility of the derived metabolic network (Figure 1) as a model for metabolism during mycorrhizal interaction was evaluated by analysis expressed sequences for core metabolic functions. The ability of the sequence data to fully represent expression of core metabolic functions would indicate the suitability of the transcriptomic data to investigate all predicted metabolic pathways. The approach and results are largely consistent with a previous transcript profiling study that delineated the major pathways of carbohydrate metabolism in L. bicolor[33]. To determine if genes that are expected to be expressed in aspen and L. bicolor are detected as expressed, we consider enrichment of core metabolism KEGG maps in the observed transcriptomics data. Expressed enzyme activities in the KEGG pathways for the core metabolic pathways pyruvate metabolism, glycolysis/gluconeogenesis (map00010), citrate cycle (map000230), Fatty acid metabolism (map00071), purine metabolism (map00230), and pyrimidine metabolism (map00240) are significantly enriched (p-value less than 0.0005) in both aspen and L. bicolor transcriptomes.

Identification of unique metabolic capabilities

Another set of unique activities in L. bicolor is the capability for the synthesis and utilization of allantoin via ureidoglycolate or urea (illustrated in Additional File 1,figure S1-b and S1-c). This capability is represented by the uric acid oxidase, allantoinase, allantoicase, and urease proteins encoded in the L. bicolor genome and expressed in the mycorrhizal metabolome. The assigned enzyme activities and annotations are supported by sequence alignments to conserved domain database profiles (CDD of NCBI; http://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml), for these enzymes [38]. The L. bicolor urate oxidase protein sequence (EC 1.7.3.3, Protein ID 190858) generates a position-specific scoring matrix (PSSM) alignment to the conserved domain for urate oxidase (d00445) with and E-value of 2.03e-61. The L. bicolor protein sequences for allantoinase (EC3.5.2.5), allantoicase (EC3.5.3.4), and urease (EC3.5.1.5) generate similar PSSM alignments consistent with the assigned activities (allantoinase,TIGR03178, E-value of 8.08e-121; allantoicase, TIGR02961, E-value of 4.71e-99; and urease TIGR01792, E-value of 0e+00. These observations suggest that the main nutrients L. bicolor exchanges with aspen roots for sugars are nitrogen-based compounds more complex than the ammonium salts already present in the media [26, 27].

The predicted mycorrhizal metabolome

Transcriptome data mapped at the level of KEGG pathways enabled identification of unique metabolic capabilities contributed by aspen and L. bicolor to the predicted mycorrhizal metabolome. To identify a priori the specific molecular mechanisms of those contributions it is necessary to consider the mycorrhizal metabolome model as a complete network of compartmentalized metabolisms interacting by transported metabolites through expressed membrane transporters. The process is enabled by identification of regions in the model represented in Figure 1 that exhibit unique interconnecting enzyme activities. Particular emphasis is assigned to those regions of unique activities that are also associated with compounds predicted to be transported by expressed, enriched amino acid transporters or sugar porters. Those identified regions underscore the most likely metabolic compounds that are synthesized by one of the ECM partners and transported across the apoplast to the partner organism. Several regions are identified using this strategy and lead to the identification of five compounds; fructose, sucrose, glycine, glutamate, and allantoin (highlighted in pink in Figure 1 and in interactive Additional File 1, Figure S1). These compounds meet the three criteria imposed for analysis of the predicted mycorrhizal metabolome; 1) the activities are uniquely expressed in one of the organisms, 2) the synthesized compound can be matched to appropriate exchange transporters, and 3) there is a reciprocal pathway in the partner organisms to enable utilization of the synthesized compound.

Laboratory mycorrhizal cultures of L. bicolor and aspen seedlings

Mycorrhizal aspen and L. bicolor were grown in biological duplicates. L. bicolor (Maire) Orton (strain S238N) culture was maintained on Modified Melin Norkan's (MMN) media as described [50] at 20°C. Populus tremuloides seeds were surface sterilized and germinated on McCown's woody plant media (WPM) in Petri dishes as described [51]. One week old germinated seedlings were transferred to Magenta vessels (Sigma, St. Louis, MO) containing the interaction medium (WPM with 1.5% of sucrose). The seedlings were grown under 16 hr light and 8 hr dark cycles at 24°C for 4-5 weeks until fine lateral roots were developed. L. bicolor mycelial plugs were transferred to Magenta vessels with aspen plants for ectomycorrhizae formation. Mycorrhizal roots were visualized approximately after 6 weeks. Mycorrhizal roots were collected, snap frozen in liquid nitrogen and stored at -80°C. Total RNA was extracted from mycorrhizal samples by CTAB method [52] and RNA quality was assessed by gel electrophoresis prior to library preparation. Total RNA was treated with RQ-DNase (Promega, Madison, WI) to remove any traces of DNA. The transcriptomic sequence was derived from the two independent biological samples but the data sets were combined for this analysis. Analysis of the independent set indicated 94% of aspen reads detected as expressed in one replicate were detected in the other. Significance of overlap for aspen is a p-value less than 1e-256. 95% of L. bicolor reads detected as expressed in one replicate were detected in the other. Significance for overlap for L. bicolor is a p-value less than 1e-256.

mRNA-Seq sample preparation

Procedures described for preparation of mRNA for the mouse transcriptome analysis [53] were used with some modifications. Ten ug of total RNA from each sample was hybridized to Sera-mag oligo (dT) beads (Thermo Scientific) for mRNA purification. Purified mRNA was fragmented by addition of 5X fragmentation buffer (Illumina, Hayward, CA) and was heated for 5 min at 94° C in a thermocycler. First strand cDNA was synthesized using random primers to eliminate the general bias towards 3' end of the transcript. Second strand cDNA synthesis was done by adding GEX second strand buffer (Illumina, Hayward, CA), dNTPs, RNaseH and DNA Polymerase I followed by incubation for 2.5 h at 16° C. Second strand cDNA was further subjected to end repair, A-tailing, and adapter ligation in accordance with the manufacturer supplied protocols. Purified cDNA templates were enriched by PCR amplification with Phusion DNA polymerase (Illumina, Hayward, CA) and the samples were cleaned using QIAquick PCR purification columns and eluted in 30 μl EB buffer as per manufacturer's instructions (QIAGEN, CA). Purified cDNA libraries were quantified using Nanodrop spectrophotometer and loaded onto Illumina flow cells with one of the flow cells reserved for a recurrent internal standard as control for sequencing efficiency. Over 25 million sequence reads were generated for analysis.

Detecting significant gene model expression

To detect significant gene model expression in collected transcriptomics data, we used the application 'BowStrap'. 'BowStrap' is a bootstrap style application of the available short sequence-aligning program Bowtie. As not every generated sequence read that does align can be uniquely aligned to a single location in a gene model, BowStrap takes into account those short sequence reads that align to more than one gene model location and calculated gene expression as number of aligned RNAseq Reads per Kilobase per Million reads (RRKM) with an associated standard error. The ultra-fast sequence alignment program "Bowtie-0.9.9" [54] was used to generate alignments of sequence reads to genomes and gene models. Bowtie uses a Burrows-Wheeler index to rapidly align sequences to a pre-processed indexed set of sequences. When Bowtie encounters a short sequence that does not align to a single location, Bowtie can return all the possible alignments for that sequence read which is required by 'Bowstrap'. 'BowStrap' reads the Bowtie output, and for some selected number of iterations, calculates an RPKM value for each gene model. At each iteration, if a RNA seq read was aligned to more than one possible location, that read is randomly assigned to one of those potential locations. From the set of RPKM values calculated across all Bowstrap iterations, an averages and standard deviations of RPKM are determined for each gene model. A statistical measurement to determine if expression of a gene models has been detected at significant levels is calculated as a Cumulative Normal Distribution (CND) and expressed as a pValue (Eq. 1)

(1)

Where μ is the average and σ is the standard deviation of a re-sampled gene model's RPKM. As this equation will fail to return a value if the average or standard deviation is equal to 0, if the μ is equal to 0, then the CND pValue is set equal to 1, else if σ is equal to 0, then the CND pValue is set equal to 0. A CND pValue close to zero indicates statistically significant levels of detected gene expression.

In our analysis, all sequence reads were trimmed to 46 bp. Optimal read length for expression analysis is impacted by system and sequence characteristics [55] but the 46 bp read length was empirically determined to yield good sequence alignments to L. bicolor and aspen gene models. Bowtie indexes were generated from sets of published JGI gene models for L. bicolor and aspen. The default Bowtie conditions were used to generate alignments for all sets of sequence reads to gene models, except for setting Bowtie to return all possible sequence alignments. We used 1000 iterations in 'BowStrap' for the calculation of average and standard deviations of RPKM values. The source code for a PERL implementation of 'BowStrap' is available at [http://www.bio.anl.gov/molecular_and_systems_biology/bowstrap/Bowstrap_download.html]. 'Bowtie' and 'BowStrap' output files are available at the Gene Expression Omnibus (GEO) as series record GSE28157 (http://www.ncbi.nlm.nih.gov/geo/ ).

Analysis of alignments generated with Bowtie indicated a total of 53% of the total sequence reads aligned to L. bicolor or aspen gene models. 55.2% of those aligned reads were to L. bicolor gene models and 45.8% to aspen gene models. A small fraction of the total aligned reads, 0.49%, aligned to both L. bicolor and aspen gene models. Of the entire set of aligned reads, 82% aligned to a single location with the balance aligned to two or more possible locations in the set of gene models. Alignment of reads to the genome scaffolds indicate 57% of the sequence reads aligned to L. bicolor or aspen genomic sequence. We observed that 3.2% of sequences aligned to gene models but not to genomic sequence indicating that these sequences align across gene model splice junctions. A fairly high fraction (11.1%) of sequence reads align to genomic sequence but not to gene models. These aligned sequences may indicate the presence of genes that are expressed but have not been identified in the set of gene models or errors in the published gene models. This is consistent with a recently reannotation of the L. bicolor genome which increased the number of gene models by ~6% (JGI gene portal). The revised annotation may also account for a fraction of the 36.1% of sequences that align to neither gene models nor genomic sequence. Some of these unaligned reads may come from errors or deviations in the assembled L. bicolor and aspen genomic sequences. Some fraction of the unaligned sequence reads likely derive from erroneously called or multiple splice sites in the set of gene models. Previous analysis of deep RNA sequencing in L. bicolor have identified that up to 60% of L. bicolor gene models contain at least one deviation from expressed gene sequences [56]. The last likely source of un-aligned sequence reads is sequencing errors and the error rates inherent in Illumina sequencing have been previously well characterized [57–59]. A single insertion or deletion in a sequence read will render it un-alignable under the 'Bowtie' criteria used here. More than two substitution errors in a sequence read will also prevent it from finding an alignment. aspen The observed percentage of reads aligned at a high stringency to the gene models is comparable to similar transcriptomic experiments [60–63]. To generate a measurement of the specificity of reads that were found to align to gene models, we generated a synthetic set of reads that were the reverse, but not the reverse compliment, of all observed sequences reads (entire dataset). This synthetic set of reads retains the frequency and relative distribution of nucleotides from the observed data. Aligning the synthetic read set to gene models identifies a total of 28 (0.0001%) aligned reads. Inspection of the aligned synthetic reads indicated that they are comprised of 2 or 3-mer repeats. Alignment with the synthetic read set demonstrates that observed alignments are highly specific and miss-alignments are not expected.

In our transcriptomic data, at a CND-pValue < 0.001, for aspen there are 27318 expressed gene models out of a total possible 45555 (60.0%) gene models. For L. bicolor, there are 12856 significantly expressed gene models out of 26014 (62.4%).

Gene Model Annotation

Annotation to gene model predicted proteins was provided by JGI publicly available downloads using the following files: "Lacc.ecAnnot.txt" for (EC) annotation for L. bicolor best gene models, "Lbicolor_GO.txt" for Gene Ontology (GO) annotation for L. bicolor best gene models, "mntPop.ec.Annotation.txt" for Enzyme Commission (EC) annotation for aspen gene models, and "Poptr1_1.goinfo.tab" for Gene Ontology (GO) annotation for aspen gene models. The genome annotation resources used in this study were derived by a combination of automated and manual approaches that are described in the Additional material associated with the genome sequence publications.

Using a strategy employed by other investigators [64–66], we used gene models from the JGI sequenced and closely related P. trichocarpa as surrogates for the P. tremulous genes. There are 12,813 published P. tremuloides ESTs in NCBI with an average length of 472 bp, for which the longest sequence is 848 bp and the shortest sequence is 101bp. The limited number of available EST sequences does not represent a complete catalogue of P. tremuloides genes and the average sequence length limits functional prediction approaches. To verify that this is a reasonable substitution, we BLAST the P. tremuloides ESTs against the JGI predicted best gene models for P. trichocarpa. 95% of P. tremuloides ESTs hit JGI P. trichocarpa gene models at an e-value threshold of 1e-6. The average percent identity is 93% with a minimum observed 73% identity. 422 (3.3%) P. tremuloides ESTs share 100% identity with P. trichocarpa gene models. While the species are not identical and the possibility of presence of unique genes in one species cannot be completely dismissed, the level of similarity observed is suitable for a system-scale analysis of the P. tremuloides transcriptome.

10% of L. bicolor gene models have EC number annotations and 37% have GO Molecular Function annotation. 13% of aspen gene models have EC number annotations and 41% have GO Molecular Function annotations. Gene models with GO annotations amino acid transport (GO:0006865), Ammonium transporter activity (GO:0008519), Sugar porter activity (GO:0005351), and Inorganic phosphate transporter activity (GO:0005315) are the relevant transporters in the mycorrhizal metabolome.

Predict mycorrhizal metabolome from transcriptomic data

In this manuscript, we use the term "enzyme function" to describe a specific annotation applied to an enzyme, i.e. "Phosphotransferases with an alcohol group as acceptor". We use "enzyme reactions" to refer to metabolite transformations catalyzed by an enzyme function, i.e "ATP + D-Glycerate ↔ ADP + 3-Phospho-D-glycerate". A unique enzyme function may catalyze more than one enzyme reaction and an enzyme reaction may be catalyzed by more than one unique enzyme function. A metabolite is a molecular compound that is a reactant or product in an enzyme reaction. A network, constructed from annotated unique enzyme functions, enzyme reactions and inferred metabolites, is used to generate the predicted mycorrhizal metabolome. The set of metabolic reactions in KEGG databases was used to represent the set of all possible enzyme reactions in the mycorrhizal metabolome and Enzyme Commission (EC) number annotations for enzyme function for genes models were used to assign unique enzyme functions to the predicted proteins encoded in the mycorrhizal transcriptome. A unique enzyme reaction was included in the mycorrhizal metabolome if a gene annotated to an enzyme function that catalyzed that reaction was detected as significantly expressed in the mycorrhizal transcriptome. All metabolites associated with included enzyme reactions were included in the predicted mycorrhizal metabolome.

Map genomic and expressed gene models to KEGG metabolic pathways

Annotation of proteins encoded by predicted gene models were obtained from publically available data resource at the Joint Genome Institute (http://www.jgi.doe.gov/) ftp archive. Lists of reactions in the form of reactants, products, and mediating EC enzyme activities were acquired from KEGG (ftp://ftp.genome.jp/pub/kegg/pathway/). Sets of enzyme activities associated with individual KEGG pathways were also collected from the KEGG website. The KEGG database of metabolic interactions includes 1329 organism-specific metabolomes, including 14 complete plant and 51 fungal metabolomes. Of the 378 pathways in KEGG, 27 are plant-specific pathways. The set of plant and fungal metabolomes includes those specific to P. trichocarpa and L. bicolor. KEGG includes over 8000 individual metabolic reactions, between nearly 15 thousand metabolites and mediated by over 2500 specific EC activities. Every gene model in aspen or L. bicolor annotated with a specific EC activity was mapped to KEGG enzymatic reactions as described below.

Identify enriched annotations and pathways in sets of expressed enzymes

In this analysis, it is required to determine if a particular annotation appears in a subset of genes, relative to the frequency that annotations appear in the superset, at a frequency greater than expected by chance. For this, the Cumulative Binomial Distribution (CBD) was used (Eq. 2). The CBD as written here returns the probability that there are at most number x successes in n trials where each trial has a probability of success p.

(2)

Where x is equal to the total number of successes in a sample, n is equal to the sample size, p is equal to the probability of success of each sample, and Annot is a specific annotation. For the application here, the number of successes x is equal to the number of times a particular annotation Annot is present in a subset of gene models and the sample size n is the number of gene models in the subset. The probability p is the proportion of gene models in the superset attributed with annotation Annot. A CBD pValue close to zero indicates an enrichment of that annotation in the subset while a CBD pValue close to one indicates a depletion of that annotation relative to the distribution of the annotation in a larger set.

Significant enrichment of specific annotations in the sets of expressed mycorrhizal genes, relative to the total number of those specific annotations in the annotated genomes and aspen and L. bicolor, was calculated. For the calculation, the number of successes is equal to the number of genes with a specific annotation in the set of expressed gene models. The sample size is the total number of expressed gene models. The probability of success is the frequency of that specific annotation in the complete genome. Specific annotations considered were gene models with an EC activity, sugar transporters, nitrogen transporters, and inorganic phosphate transporters.

KEGG pathways were identified if they were enriched for expressed EC activities in common between the sets of aspen and L. bicolor expressed gene models. For each KEGG pathway, the number of successes is equal to the number of EC activities in the set of expressed aspen and L. bicolor gene models that map to that pathway. The number of trial is the total number of EC activities expressed in mycorrhizae by aspen and L. bicolor. The probability of success is equal to the number of EC activities associated with that KEGG pathway divided by the total number of EC activities expressed by aspen and L. bicolor. Specific KEGG pathways were also identified if they were statistically significantly enriched for mapped genes with EC activities unique to expression in either aspen or L. bicolor. In this analysis, for each KEGG pathway, the number of successes is equal to the number of gene models in the set of expressed aspen or L. bicolor gene models that map to that pathway. The number of samples is equal to the total number of expressed gene models in mycorrhizae by aspen or L. bicolor. The probability of success is equal to the number of aspen or L. bicolor gene models that map to that KEGG pathway divided by the total number of expressed gene models by aspen or L. bicolor. Enriched KEGG pathways were reported at a significance of p < 0.01, and if at least 25% of the genes or activities associated with that KEGG pathway were also in the complete genomic set of aspen or L. bicolor gene models.

The complete set of KEGG pathway EC activities was collected from ftp://ftp.genome.jp/pub/kegg/pathway/ec/. Identifying expressed activities and transported small molecules in KEGG global metabolism was performed using the KEGG PATHWAY Query tool.