Deciphering herbivory-induced gene-to-metabolite dynamics in Nicotiana attenuata tissues using a multifactorial approach

In response to biotic stresses, such as herbivore attack, plants reorganize their transcriptomes and reconfigure their physiologies not only in attacked tissues but throughout the plant. These whole-93 organismic reconfigurations are coordinated by a poorly understood network of signal 94 transduction cascades. To explore tissue-based interdependencies in Nicotiana attenuata’s 95 resistance to insect attack, we conducted time-series transcriptome and metabolome profiling of 96 herbivory-elicited source leaves and unelicited sink leaves and roots. To probe the 97 multidimensionality of these molecular responses, we designed a novel approach of combining an 98 extended Self-Organizing Maps (SOM) based dimensionality reduction method with bootstrap-99 based non-parametric ANOVA models to identify the onset and context of signaling and 100 metabolic pathway activations. We illustrate the value of this analysis by revisiting dynamic 101 changes in the expression of regulatory and structural genes of the oxylipin pathway and by studying nonlinearities in gene-metabolite associations involved in the acyclic diterpene glucoside pathway after selectively extracting modules based on their dynamic response patterns. 104 This novel dimensionality reduction approach is broadly applicable to capture the dynamic 105 rewiring of gene and metabolite networks in experimental design with multiple factors. 106 107 co-expression network representing first neighbors of identified HGL-DTG in N. attenuata. The cluster represents the strong connectivity between HGL-DTG metabolites and known biosynthetic genes in N. attenuata. Other previously characterized m/z features present in this cluster are enriched for O-acyl sugars associated signals.

based non-parametric ANOVA models to identify the onset and context of signaling and  This novel dimensionality reduction approach is broadly applicable to capture the dynamic  which involves intricate signaling pathways (Hahlbrock et al., 2003;Nakashima et al., 2009;1 2 5 Zeller et al., 2009;Walley and Dehesh, 2010). These transcriptional adjustments can be captured 1 2 6 by studying changes in the expression of genes in different tissues in order to elucidate the 1 2 7 influence of particular pathways as well as the relative contribution of a given tissue to the 1 2 8 whole-organism's response. Although poorly understood, signaling networks controlling these 1 2 9 transcriptional responses have been shown to be highly stress-condition specific, as clearly mechanical wounding and herbivory (Baldwin, 1990;Alborn et al., 1997;Halitschke et al., 2003; 1 3 2 Reymond et al., 2004;Wu et al., 2007). Experiments designed to study such intricate networks  numbers of false positives (Bittner et al., 1999;Getz et al., 2000). First, clustering algorithms subset of the experimental conditions (Swindell, 2006). Additionally, connections in GRN ( Priness et al., 2007) and Biclustering (Dharan and Nair, 2009) have been developed to address 1 5 0 this limitation. Moreover, interaction patterns deduced from these co-expression studies represent 1 5 1 the static wiring of the network, whereas networks will assemble dynamically as the organism 1 5 2 adapts to external stimuli. To best capture the temporal dimension as a variable affecting the structure of GRN,  performance of entire parts of the OS-elicited transcriptome, we isolated interactive motifs from 4 1 3 the SOM grids and analyzed their dynamic behaviors using network analysis.

1 4
We used this network-based approach to analyze the regulation of OS-elicited changes in 2008) as well as three other genes (NaHDR, NaHDS and NaISPD) predicted based on their 4 2 4 homology to Arabidopsis isoprenoid genes were mapped to the SOM grid, and found to be  Using a module-centric approach, we extracted this interactive motif to construct a gene-4 3 0 gene network using a statistically sound two-stage co-expression detection algorithm with an 4 3 1 FDR=0.05 and a minimal acceptable strength of 0.7 using the GeneNT package in R ( Figure 5A).

3 2
We found that this network shares properties of other biological networks. The scale-free 4 3 3 topology is of 0.96, which suggests that the network is comprised of many genes with few 4 3 4 connections, but a few genes with many connections. Additionally, the clustering coefficient, which provides a measure of cliquishness (0.516), and the measure of network heterogeneity, 4 3 6 which reflects the tendency of a network to contain hub nodes (1.40), are within ranges expected 4 3 7 for biological networks (Supplemental Figure 5). Interestingly, we observed the highest degree of  Network analysis also confirmed differential functional activation of the three NaGGPPS genes 4 4 2 present in N. attenuata, with NaGGPPS1 showing no significant interactive effect in TvS 4 4 3 comparison while the other two showed large differences in their degree of connectivity in the 4 4 4 motif 5c network (deg(NaGGPPS2)=224, deg(NaGGPPS3)=90).

5
To infer OS-elicitation specific gene associations, we constructed a partial co-expression 4 4 6 network using the first neighbors of NaGGPPS2 for the W+W treatment from the set of genes  To further support this conclusion, we conducted a single time point DEG (differentially leaves ( Figure 5B). As expected, major processes highlighted by this analysis included secondary and genes from photosynthetic pathways as down-regulated processes in treated leaf tissue. The 4 6 1 biosynthesis of plastid isoprenoids by the non-mevalonate is essential for photosynthesis and 4 6 2 chloroplast function (Vranova et al., 2012). Dense connectivities were shared between the non-4 6 3 mevalonate/17-HGL-DTG pathways and photosynthetic genes suggesting yet-unknown 4 6 4 coordination mechanisms between these two processes and consistent with the previously  To facilitate the interpretation of the large transcriptomic differences observed between 4 7 1 elicited and unelicited tissues, we profiled downstream metabolite responses using a broadly-4 7 2 targeted UHPLC-qTOFMS metabolomic approach (Gaquerel et al., 2010). Often in such large 1 6 scale profiling, several compounds are not completely chromatographically resolved or are 4 7 4 observed with shift in their retention times, so mathematical procedures involving deconvolution 4 7 5 and retention time corrections need to be applied to extract accurate mass spectra with resolved modes using the XCMS package in R.

8 0
The first stage in the identification of differences between local and systemic responses is tissue type (treated and untreated leaves). As with the transcriptomic analysis and using the same 4 8 4 parameters (FDR=0.05, #bootstraps=200), the sampling time was not taken as another factor but 4 8 5 was used to find the best ANOVA structure along the optimal direction resulting in response metrics. In addition to the induced changes in responses to OS (treatment effect), a set of 4 8 7 metabolites of interest also showed small differences between the two leaf tissues (tissue effect); 4 8 8 hence they were classified in an additive bin. Therefore, in order to increase metabolite coverage,  from the negative mode showed major differences in induced responses between leaf tissues for 4 9 2 at least one time point in the series (Supplemental Files 2 and 3). This is many folds fewer than 4 9 3 transcriptomic differences observed between responses of the two leaf tissues. One reason for this 4 9 4 lower coverage could be the use of an targeted approach for profiling major classes of secondary 4 9 5 metabolites; moreover, some metabolites showed large constitutive biological variation which 4 9 6 thwarted the detection of significant differences between control and W+OS condition with an 4 5 3 1 power of this approach to cluster biochemically-connected metabolites, we additionally observed 5 3 2 highly coordinated dynamic responses between shikimate pathway-derived amino acids and  Time-lag corrected correlation on the interactive response metric supports the 5 3 6 reconstruction of gene-to-metabolite networks.

8
The nature of the co-regulation of functionally related genes and metabolites could vary 5 3 8 depending on biological activity of the studied metabolic class and the experimental conditions TvS comparisons using the interactive responses of genes and metabolites as associative metrics and used this approach to study previously reviewed oxylipin and 17-HGL-DTG pathways.

4 7
Although a tight correlation pattern is usually not expected for metabolites that are rapidly consumed by subsequent reactions, as it is the case for jasmonate production, we observed that genes of the LOX3 pathway correlated well with those in JA levels as well as associated underlying genes. Therefore to address this time lag, we applied lagged Pearson Correlation (PC), 5 5 7 estimated as following: We constructed a gene-to-metabolite network using lagged-time-specific data with PC>=0.98. signals derived from 17-HGL-DTGs. We observed strong coordination between17-HGL-DTGs, genes. Unknown m/z signals reported in this network represent potential molecular fragment of 5 6 7 the aforementioned compounds that will need to be confirmed by additional mass spectrometry- The comprehensive classification of leaf and root herbivory-activated gene expression 5 7 2 levels presented here significantly expands our knowledge on the multi-dimensionality of metabolite associations with a high degree of predictive power. Spatial-temporal maps produced in this study not only underscore the high plasticity of OS-elicitation responses but also provide a involved in antiherbivory processes. The two metabolic routes investigated in greater detail in shared patterns of regulation with other physiological processes.

8 3
Elucidating groups of genes involved in the sequential reorganization of biological 5 8 4 networks is extremely challenging using the available bioinformatics approaches, but we considering the ordering of time points before deriving the optimal direction for best ANOVA 5 9 1 structure. Here, we demonstrate that SOM-based herbivory-elicited interactive motifs are highly 5 9 2 tissue-specific, highly versatile and decrease in size for later time points reflecting specific 5 9 3 sequential responses deployed in different tissues. We speculate that this increased modularity in several groups indicating that insect herbivory results in major shifts in almost every aspect of a 5 9 7 plant's physiology (Hermsmeier et al, 2001;Schittko et al, 2001;Schmidt et al, 2001). The fact 5 9 8 that large sized interactive motifs are more pronounced in systemic tissues at 5, 9 and 13h after 5 9 9 elicitation indicates the involvement of organismic-level responses specifically triggered by OS-6 0 0 based signaling, with a large proportion of these also attributable to unexplored root-specific Analysis of interactive effects and response timings for jasmonates and oxylipin genes 6 0 8 illustrates that the extraction of the multifactorial statistical information in terms of response 6 0 9 patterns not only facilitates the clustering of genes involved in similar biochemical pathways --as NaAOC and β -oxidation components), associated upstream --MAPK (NaSIPK and NaWIPK), contrast with the pulse effects are however consistent with the existence, as suggested by loops modulating transition points in the pathway activation. To analyze the metabolic output of these large transcriptional changes, we also profiled 6 3 1 metabolites for the same experimental conditions and tissue types ( Figure 6). Even though 6 3 2 instrumental in highlighting the breadth of the complete metabolic profile being affected by the 2 1 treatment, the mining of these data is extremely challenging, notably due to compound-specific  provider for the production of structurally diverse plastidic isoprenoids such as carotenoids, phytol (a side-chain of chlorophylls), isoprene, mono-and diterpenes (Vranova et al., 2012). In N.  Organization of molecular processes within this OS-specific gene network supports known this study relies on the biogenesis of the GGPP backbone by a specific GGPPS, NaGGPPS2 hub region of a gene network, while the two other NaGGPPS copies (NaGGPPS1 and 6 6 0 NaGGPPS3) exhibit reduced connectivity with this hub region. We predicted the subcellular 2 2 mTP, and sTP. These predictions further support our findings on the differential activation of 6 6 6 these three GGPPS genes in response to OS-elicitation.

7
The ability to reconstruct tissue-specific gene-to-metabolite dynamics with high 6 6 8 confidence underscores the importance of extracted time response patterns for features of interest, 6 6 9 which single point analysis or pooled data using time series as third factor would fail to identify 6 7 0 (Supplemental Figure 6). This novel method of combining multifactorial analysis with the 6 7 1 information extracted from time series data is broadly applicable to investigate 6 7 2 signaling/metabolic pathways in other biological systems with time series data to deduce the processes. This strategy can further be used for improving clustering analyses by applying a Germany) with Klasmann plug soil (Klasmann-Deilmann GmbH, Geesten, Germany) and after 6 8 6 10 days, seedlings were transferred to 1 L pots with sand to facilitate sampling of roots. Plants research and interpreted results, wrote the paper. References        Figure 1: Work-flow for analyzing changes in the transcriptome/metabolome landscape of Nicotiana attenuata elicited by insect herbivory. Leaf and root tissues of wild type Nicotiana attenuata plants were harvested 1, 5, 9, 13, 17 and 21 h after leaves were wounded with a fabric pattern wheel and immediately treated with Manduca sexta oral secretions (W+OS) to study herbivory induced responses or water (W+W) to study OS specific induction. Replicated transcriptomic and metabolomic data were analyzed using multifactorial analysis with both factors (tissue and treatment) taken together across the time series to identify modules showing differential OS-elicitation. These modules were spatial-temporally resolved using batch-learning self organizing maps (BL-SOM). Herbivory-elicited interactive effects among transcript and metabolite levels were analyzed with a special emphasis on identifying short-and mediumterm changes in metabolism. To this end, interactive motifs from resolved maps were extracted and analyzed using network properties.     (NaGLA1) is oxygenated by NaLOX3 enzymes (LOX) and converted by a multi-enzymatic cascade into JA. JA-Ile formed after enzyme-dependent conjugation of JA to isoleucine and is the bioactive jasmonate interacting with the F-box receptor protein NaCOI1 and NaJAZ transcriptional repressor proteins. C12 and C6 aliphatic chains are produced by cleavage of NaLOX2-dependent hydroperoxides. Divinylether (DVE) oxylipins are produced by enzymatic rearrangement of NaLOX1-dependent hydroperoxides. (B) Hierarchical cluster analysis of ANOVA directions for interactive effects for genes involved in the two different branches of the lipoxygenase pathway. Elements of the NaLOX3 branch (red) exhibit differential regulation in leaf tissues 1h after OS elicitation. NaLOX2 and NaLOX1 branches (green) exhibit high interactive effects 5h after elicitation. (C) Temporal profiles for genes from NaLOX3 pathway showing a high degree of coordination at 1h while temporal profiles for NaLOX2 and NaLOX1 pathway show strong coordination at 5h.