Abstract

Background

Climatic and edaphic conditions over geological timescales have generated enormous diversity of adaptive traits and high speciation within the genus Eucalyptus (L. Hér.). Eucalypt species occur from high rainfall to semi-arid zones and from the tropics to latitudes as high as 43°S. Despite several morphological and metabolomic characterizations, little is known regarding gene expression differences that underpin differences in tolerance to environmental change. Using species of contrasting taxonomy, morphology and physiology (E. globulus and E. cladocalyx), this study combines physiological characterizations with ‘second-generation’ sequencing to identify key genes involved in eucalypt responses to medium-term water limitation.

Results

One hundred twenty Million high-quality HiSeq reads were created from 14 tissue samples in plants that had been successfully subjected to a water deficit treatment or a well-watered control. Alignment to the E. grandis genome saw 23,623 genes of which 468 exhibited differential expression (FDR < 0.01) in one or both ecotypes in response to the treatment. Further analysis identified 80 genes that demonstrated a significant species-specific response of which 74 were linked to the ‘dry’ species E. cladocalyx where 23 of these genes were uncharacterised. The majority (approximately 80%) of these differentially expressed genes, were expressed in stem tissue. Key genes that differentiated species responses were linked to photoprotection/redox balance, phytohormone/signalling, primary photosynthesis/cellular metabolism and secondary metabolism based on plant metabolic pathway network analysis.

Conclusion

These results highlight a more definitive response to water deficit by a ‘dry’ climate eucalypt, particularly in stem tissue, identifying key pathways and associated genes that are responsible for the differences between ‘wet’ and ‘dry’ climate eucalypts. This knowledge provides the opportunity to further investigate and understand the mechanisms and genetic variation linked to this important environmental response that will assist with genomic efforts in managing native populations as well as in tree improvement programs under future climate scenarios.

Keywords

EucalyptusTranscriptomicWater deficitHomeostasisEcotype

Background

Changes in climate over geological timescales have played a major role in the evolution of the Australian vegetation. Northward continental drift, beginning in the late Oligocence (38–23 Ma), led to aridity cycles of increasing severity and duration during the Miocene (25–10 Ma) which have promoted a range of adaptive traits characteristic of the Australian vegetation. Isolation of gene pools via developmental, climatic and edaphic factors, coupled with high degrees of outcrossing has promoted speciation among many Australian plant families, most notably the genus Eucalyptus L’Hér genus giving rise to >900 species. The Genus Euclayptus is widely distributed across Australia in tropical, arid, temperate and Mediterranean climatic zones occupying an extensive ecological range from regions with annual rainfall as low as 250 mm [14], to the wet sclerophyll forests of eastern and south-western Australia where rainfall exceeds 1500 mm year−1.

Among a range of edaphic and climatic factors governing species distributions, water availability is a significant determinant of eucalypt growth, survival and productivity (see: [1]). A range of adaptive and acclimation traits to enhance growth and survival under the effects of water deficit are prevalent including anatomical changes in leaf shape and area [55], changes in biomass allocation [47], variable stomatal control [27, 56], differences in cell wall reinforcement [20, 22], cell wall water storage [51] and cellular osmolarity [9, 21, 33, 35, 40]. These properties act in concert in many eucalypt species, undoubtedly combining to produce high levels of site-specific variation in performance and survival.

The primary function of physiological and chemical adaptations, either by tolerance or avoidance, ultimately serves to maintain cellular function. The scope of environmental conditions tolerated by the genus Eucalyptus, and the diversity of physiological responses to the availability of water detected among contrasting taxonomic groups [5] suggests both isohydric and anisohydric approaches to plant growth and survival are prevalent within the Genus. The concomitant effects of water deficit, such as redox imbalances within photosynthetic reactions or the capacity to dissipate excessive leaf temperatures require the enhancement of homeostatic mechanisms to maintain cellular function.

The complexity and interrelatedness of biochemical and physiological pathways of plants presents a major challenge to our ability to study plant scale responses. More recent technological advances to collect network scale molecular data [13] have led to progress in the adoption of a ‘systems biology’ approach (e.g. [26]) to the study of plants that is still being defined [58]. In addition to the scope of challenges faced by ‘omic’ approaches to comprehensively encompass spatio-temporal variation, interpretation of such datasets is challenging as the development of suitable bioinformatic tools lags behind our ability to acquire data, especially for non-model species. Emerging technologies in bioinformatics to enable a systems approach to multidisciplinary investigations will undoubtedly prove vital in the development of ‘plant systems biology’ and the investigation of dynamic systems from the gene to the landscape scale.

For Eucalyptus, understanding physiological and chemical responses of species is presently based on stochastic investigations of pathways that are part of a larger dynamic network. The molecular basis of acclimation offers a logical reductionist approach to investigating network complexity and is now receiving significant attention (e.g. [15, 39]). In plant systems, transcriptional regulation of pathways governing a cascade of chemical and physiological response mechanisms underpinning homeostasis are typically considered in isolation with little context to more general adjustments in metabolism and growth. Investigations collecting ‘omic’ data (for review see: [13]), such as the ‘transcriptome’, therefore offers considerable advantages to guide physiological investigations and provide upstream information on regulatory processes in a broader context of plant scale homeostasis. Whilst it is well known that quantitative and qualitative regulation of plant metabolism exists throughout the framework, transcript libraries offer a logical ‘first step’ in developing higher-level understanding of whole plant responses to environmental change.

Despite a broad understanding of physiological and chemical traits in Eucalyptus species conferring acclimation to the effects of aridity, less is known regarding the molecular basis of such traits. More recent investigations at both the genomic (e.g. [15]) and transcriptome level [54] are proving highly informative for our understanding of ‘omic’ responses by eucalypts in adapting to and mitigating against changing growth conditions. Studies into water stress response in model and crop species, which are generally short lived annuals with low tolerance to water deficit, tend to focus on identifying early gene expression changes that trigger the initial responses to a change in water status (e.g. early response transcription factors). In contrast, eucalypt species are long lived perennials thus subjected to seasonal changes in water availability that can vary in severity from year to year depending on climatic zones and prevailing weather conditions [31]. Tolerating such temporal variability and amplitude of osmotic potentials [4, 31, 57] is crucial for survival and reproduction.

Adopting a comparative transcriptomics approach (e.g. [3]), though a combination of experimental design, next generation sequencing and novel bioinformatics approaches, we investigate the gene expression profiles of two Eucalyptus species of contrasting taxonomy under controlled and quantified medium term water deficit treatment in order to understand the gene sets that become stably established in responses to prolonged water stress in each species. The species used in this study, Eucalyptus globulus Labill. and Eucalyptus cladocalyx F. Muell., have contrasting properties for tolerating water deficit with E. globulus occupying more ‘mesic’ and E. cladocalyx more ‘xeric’ environments [35] with differences in isohydric/anisohydric behaviour observed in traits such as stomatal control [35] and both qualitative and quantitative chemical traits for homeostatic regulation of the cellular environment [32, 34]. Despite these differences, both species are anatomically similar in leaf and plant structure. On this basis, this comparative transcriptome study facilitates an investigation of transcriptome wide differential gene expression on a background of Eucalyptus responses to the effects of medium term water deficit to provide insight into the molecular mechanisms that differentiate the ‘mesic’ and ‘xeric’ ecotypes.

Results

Leaf gas exchange in control (WW) plants was consistent for both species throughout the study compared to severe stress (SS) treated plants, with average net carbon assimilation values higher for E. globulus than for E. cladocalyx (Table 1). Both net carbon assimilation and stomatal conductance reduced with treatment intensity. Sub-stomatal carbon concentration (ci), similarly reduced with treatment indicating that limitations to carbon assimilation are largely attributable to stomatal limitations to diffusion rather than reductions in carboxylation arising to biochemical factors (Table 1). From these results it is clear that the main limitation is water availability in line with the treatments applied.

Sequencing data

Of the 272,167,757 raw reads 240,491,307 passed filtering as high quality (88%; Paired HQ 208,115,303). Of the HQ reads 120,719,749 aligned to 23,623 of the 33,916 annotated E. grandis v1.0 genes and these were used in the differential expression analysis. Most of the remaining reads either aligned to the plastid genomes (Stem Tissues Average 327,506 [3.7%]; Leaf Tissues Average 1,876,819 [24.4%]) or to ribosomal genes (Average 2,044,031 [24%]) with on average 85.5% (s.d. 4.1%) of HQ reads explained per sample (Additional file 1: Dataset S1). Investigation of the unmapped reads indicated these largely derived from unannotated ribosomal genes with ~94.7% of reads explained when mapped to the genome sequence as opposed to annotated genes. A measure of the inter-library variation between replicates (reflecting the biological variability in the data) is the common coefficient of biological variation [28] and was estimated as 0.5711 (Gene-wise estimates ranged from 0.321 to 1.803). A multi-dimensional scaling analysis of the transcriptome data (Fig. 1) shows clear discrimination between the two species and between WW and SS treatments with the tissues clustering by type. This high-level view clearly shows that the physiological differences between the species and treatments, as described above, are reflected in the transcriptome and that the isohydric strategy of E. globulus is reflected in a lesser transcriptome response in all tissues.

Fig. 1

Multi-dimensional scaling plot of transcriptome data distinguishes samples by tissue and species. Multi-dimensional Scaling plot showing the main separation to be between stem derived tissues (Brown) and leaf like tissues (Green) on the x-axis and species (EG – Eucalyptus globulus; EC – E. cladocalyx) on the y-axis

Transcripts show differential ecotypic responses to water deficit

Differential expression analysis (Additional file 2: Dataset S2) aimed to identify genes showing a general response to treatment, and secondly, to identify genes with ecotype specific responses to treatment. In all cases a stringent FDR of 0.01 was used to determine significance. In total 460 genes showed a significant differential expression response to treatment, with 188 having a higher expression in the SS treatment and 272 having higher expression in the WW control. Of the 460 genes, 141 displayed a significant tissue specific response with 69% in E. cladocalyx and 71% in E. globulus exhibiting a stem tissue response.

By way of validation we compared the set of differentially expressed genes discovered in this study to those discovered in an independent water deficit experiment where similar stress and well watered treatments were applied to biological replicate plants from the same two species. As only leaf tissues were assayed in the validation study we would only expect to see validation of genes that showed a leaf tissue specific response and a global tissue response, approximately half the genes detected. Also, as the two studies differed in the level of replication we would expect a difference in power to detect differential expression. However, if the genes reported in this study were largely false discoveries it would be extremely unlikely for us to observe a significant overlap between the detected gene sets and further, the direction and magnitude of changes would also be expected to be uncorrelated. Of the 460 genes detected as significantly differentially expressed in stem and/or leaf tissue in this study (FDR <0.01) 232 genes in the validation study (leaf only) showed a significant differential expression and this change is in the same direction and of a similar magnitude in more than 99% of genes (allowing for experimental scaling; Additional file 3: Figure S1). This validation result strongly supports that the signals reported in this study are true signals and not false discoveries.

Gene ontology functional classification of differentially expressed genes by treatment and species. Gene ontology functional classifications undertaken using Panther for genes significantly more highly expressed in response to the water deficit (SS) treatment (a) and control (WW) (b) as well as for genes showing a significant species specific response to water deficit (SS) treatment (c) and control (WW) (d)

Genes showing significantly higher species specific expression (FDR < 0.01) in either the water deficit (SS) treatment or control (WW) in both E. cladocalyx and E. globulus from most to least significant

Of the 33 genes that showed a significant response to the SS treatment, 12 showed an opposing response by species and 21 showed a similar but enhanced response relative to E. globulus with 12 of these genes, including the top 3, being unknown. The expansin-like B1 (2) was the only gene family to be represented more than once. In the case of the 41 genes that showed a significant response in E. cladocalyx to the WW control, 23 showed an opposite response and 18 showed a similar but enhanced response. Ten genes in this group where unknown while other highly represented gene families included; disease resistance protein (TIR-NBS-LRR class) putative (3), NB-ARC domain-containing disease resistance protein (3) and soybean gene regulated by cold-2 (2). For the genes in E. globulus that showed a significant response to the WW control, all six exhibited an opposite response when compared to E. cladocalyx and comprised of only one unknown gene. The 20S proteasome beta subunit G1(2) and NB-ARC domain-containing disease resistance protein (2) gene families were represented more than once.

A network analysis of treatment responsive genes (FDR < 0.2) and of known pathways using Cystoscape™ revealed two main linked clusters; Cluster A, comprising 42 genes linked to 42 pathways, and cluster B, comprising of 63 genes linked to 89 pathways, as well as approximately 70 small unlinked clusters containing between 1 and 12 genes (178 in total) with links to between 1 and 14 pathways (165 in total) (Fig. 3a). Due to the complexity of these networks, a brief description of the key observations for the main clusters is detailed below with the low resolution figures provided (Fig. 3 and Additional file 5: Figure S2) being for demonstrative purposes as it is not possible to present the extent of detail in this format. Those wanting to view additional details and/or undertake further analysis of the network presented here are encouraged to download the cytoscape V3.3.0, (http://www.cytoscape.org/) and view the Cytoscape file (Additional file 6: Dataset S4) and interpretation instructions (Additional file 7).

Fig. 3

Metabolic pathway network analysis of differentially expressed genes in response to water deficit treatment. Cytoscape graph depicting genes with an FDR of <0.2 with known links to pathways in plant metabolic networks. Genes and pathway interactions are found in two main clusters, cluster ﻿a (a, i), and cluster b (a, ii) as well as numerous small ‘unlinked’ clusters (a, iii and iv). Sub-clusters are also identified within cluster a (b i, ii, iii, iv) and cluster b (c i, ii, iii, iv, v, vi)

Cluster A

Cluster A (Fig. 3ai, b) contained 12 genes from 7 families that showed a response (FDR < 0.2) to the SS treatment e.g. more highly expressed, compared to 26 genes from 14 gene families that showed a response in the WW control. Eighteen of the 42 pathways showed differential gene responses to treatment, e.g. at least one linked gene showing higher expression in the SS treatment and WW control, of which 12 genes showed a higher response to the SS treatment and 23 genes in response to the WW control. Four pathways showed an exclusive response to the SS treatment only and were linked to one gene, whereas 23 pathways showed an exclusive response to the WW control and were linked to 19 genes.

Number of these gene family members more highly expressed in SS treatment

Number of these gene family members more highly expressed in the WW control

A

i

p450 cytochrome

6

-

6

A

i

2-oxoglutarate and Fe(II)- dependent oxygenase

2

2

-

A

i

Chalcone and stilbene synthase

2

-

2

A

ii

UDP-glucosyl transferase

6

-

6

A

iii

O-methyltransferase

4

1

3

A

iii

P-loop containing nucleoside triphosphate

2

-

2

A

iii

ATP-dependent caseinolytic (Clp) protease/crotonase

2

2

-

A

iv

Cinnamyl-alcohol dehydrogenase

4

4

-

A

iv

Caffeoyl-CoA 3-O-methyltransferase

2

-

2

B

i

Rhamnose biosynthesis

2

-

2

B

i

UDP-D-glucuronate 4-epimerase

2

-

2

B

ii

Phosphofructokinase

3

1

2

B

ii

ATP-citrate lyase

2

-

2

B

iii

AMP-dependent synthetase and ligases

2

1

1

B

iii

Phospholipases

2

-

2

B

iv

AMP-dependent synthetase and ligase

2

-

2

B

v

Aldehyde dehydrogenase

4

2

2

B

v

Acyl-activating enzyme

2

1

1

B

v

UDP-glucosyl transferase

2

-

2

B

vi

Tyrosine transaminase

3

2

1

B

vi

2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase

3

-

3

B

vi

Homogentisate 1,2-dioxygenase

2

2

-

Gene families more highly represented in sub-clusters within cluster A and B identified as part of the network pathway analysis in Cytoscape (Fig. 3)

Within this cluster, nine genes showed a different response to treatment at a species level, (solid lines indicating species treatment interaction FDR <0.2) (Fig. 3b). Four of these genes showed a differential response, e.g. a greater relative change in expression between treatments in one species, in the SS treatment and five in the WW control where in all cases the change was observed in E. cladocalyx. Two genes showed a significant differential response to the treatment (FDR < 0.05) in the leaf tissue of E. cladocalyx (one in the SS treatment and one in the WW control) which were linked to nine different pathways (Table 4).

Genes showing a significantly higher species specific expression (FDR < 0.05) linked to metabolic pathway networks in E. cladocalyx seedlings subject to either a water deficit (SS) treatment or control (WW) over an 8 week period. Table notes in which treatment (TR) and tissue (TiR) the significant response was observed in

Cluster B

Cluster B (Fig. 3aii and c) contained 31 genes from 29 gene families that showed highest expression in the SS treatment and 32 genes from 20 gene families that showed a highest expression in the WW control. Of the 89 pathways present, 38 showed differential responses to the treatment, e.g. at least one linked gene showing higher expression in the WW and SS treatment, where 26 genes showed the response in the SS treatment and 32 genes in the WW control. 32 pathways showed a response only in the SS treatment and were linked to 23 genes, whereas 13 pathways linked to 13 genes showed a response in the WW control only.

Within this cluster, 18 genes showed a differential response to treatment at a species level (solid lines indicating species treatment interaction FDR < 0.2) (Fig. 3c). Eight genes showed a differential response, e.g. greater relative change in expression between treatments in one species, to the SS treatment and ten to the WW control where, in all but three cases, this response was observed in E. cladocalyx. Nine genes showed a significant differential response to the treatment (FDR < 0.05) in mostly the stem tissue of E. cladocalyx (five in the SS treatment and four in the WW control) which were linked to 13 different pathways (Table 4).

Minor clusters

Of the 70 small clusters identified (Additional file 5: Figure S2) 20 showed a significant response e.g. higher expression, as a result of the SS treatment only, 30 as a result of the WW control only and 20 showed responses in both SS treatment and WW control. These clusters comprised 178 genes from 118 gene families where 78, representing a total of 55 gene families, showed a significant response to the SS treatment and 100, representing 63 gene families, showed a significant response to the WW control.

Within these 70 clusters, 43 genes showed a notable differential response to treatment at a species level (solid lines indicating species treatment interaction FDR < 0.2) (Additional file 5: Figure S2). Nineteen of these genes showed a differential response, e.g. greater relative change in expression between treatments in one species, in response to the SS treatment and 24 in response to the WW control where in all but five cases, this response was observed in E. cladocalyx. Thirteen genes showed a significant differential response to the treatment (FDR < 0.05) in often both the leaf and stem tissue of E. cladocalyx (seven in the SS treatment and six in the WW treatment) which were linked to 16 different pathways (Table 4).

Discussion

This comparative expression profile for Eucalyptus species of contrasting ecotype under the effects of medium term water deficit provides insight into the molecular mechanisms that underpin adaptive and acclimation responses. Plant responses to treatment induced a stomatal response with concurrent physiological changes indicating that the main limitation for water deficit treated plants compared to the well watered was water availability. The treatments triggered a cascade of gene expression responses associated with the concomitant effects of water deficit with several patterns in gene expression emerging and offering insight into the comparative stress response mechanisms employed by Eucalyptus species. Notably, this study highlights the enhanced transcriptional responses observed in key pathways of species adapted to more water limiting environments.

This investigation demonstrates the importance of a plant scale transcriptomic approach encompassing multiple tissue types. Over 78% of genes identified as responding to the stress treatment occurred in the stem, demonstrating the potential for significant regulatory processes being modulated from this tissue. Previous and current approaches investigating plant responses to resource availability are focused on leaves potentially overlooking significant molecular, chemical and physiological regulation of processes governing stress responses. Such results clearly demonstrate the importance of ‘whole plant’ studies adhering to the literal interpretation of the term.

Encompassing multiple tissue types exacerbates a major challenge to ‘omic approaches – the complexity of interpreting such large volumes of data. The bioinformatic approach outlined here, has successfully highlighted patterns in transcript data and produced interpretable representations of network structure. It is important to recall that our intention for this study is not to comprehensively summarise the entire transcriptome, nor the pathway analysis constructed, but to offer insight into the transcriptional response mechanisms to applied water deficit and provide a framework for further investigation into the molecular, chemical and physiological mechanisms underpinning tolerance to water deficit in the genus.

The known unknowns

Among the most highly induced transcripts observed in this study, many were of unknown function. Whilst repositories of Eucalyptus sequence data are improving rapidly in scope (see for example, [15, 39, 54]) much work remains to complete the functional annotation of the genome, especially given the broad ecotypic scope, chemical composition and diversity of physiological and growth responses within the genus. The targeted experimental and bioinformatic approaches employed here prioritises a small set of genes as candidates for further investigation. This use of contrasting ecotypes has enabled us to distinguish common mechanisms and species-specific responses from the ‘background’ of gene expression necessary for general cellular function. This approach offers great potential for more functional, mechanistic approaches to characterising transcriptome wide responses to the onset of water deficit with implications for both upstream (genomic) and downstream (proteomic, metabolomic) responses.

Photoprotection and redox balance

Amongst the most immediate challenges to photosynthetic tissues during resource imbalances is avoiding the consequences of excessive photon flux. In this study, several patterns emerge at the network scale pertaining to the maintenance of redox balance. For example, in cluster A, a range of photoprotective response mechanisms from the phosphoenolpyruvate pathways are observed to respond including transcripts governing redox balance, most notably those of several cytochrome P450 enzymes which play a central role in the oxidation of organic compounds [41]. The role of P450s in oxidising organic compounds, or as the terminal of electron transport chains (including that of NADPH) represents significant homeostatic control over redox balance. Such control aids in the preservation of photosynthetic infrastructure and the reactions of primary photosynthesis (by providing reductant for the photosynthetic reactions) thus mitigating excess photochemical energy experienced during times of stress or alternatively in times of excessive photoassimilate production. Similarly, antioxidant systems such as the glutathione and ascorbate pathways appear regulated in cluster A and in the unclustered networks supporting previous studies on the role of antioxidant pathways and redox balance in Eucalyptus leaves (e.g. [10, 17]). Sequencing and subsequent annotation of these pathways provides solid evidence for their inducible responses in Eucalyptus tissues and together with the induction of P450 enzymes provide a useful basis for further characterisation in maintaining redox balance under conditions requiring quenching of excess photochemical energy or as an indicator of photoassimilate supply to cellular metabolism.

Phytohormones and signalling molecules

Rapid elicitation of signalling molecules represents upstream regulation of plant processes governing plant scale resource allocation and utilization. A range of genes involved in the synthesis of phytohormones is observed to respond in the present network. Genes encoding phytohormones in the gibberellin, auxin and cytokinin groups are observed, most notably among cluster A and the associated gene analogues encoding photoprotective properties. The wide distribution of phytohormonal responses across the network likely reflects the broad chemical and physiological processes that rely on hormonal regulation in higher plants. Of particular note is the up regulation of genes governing the synthesis of quercitin in E. cladocalyx. Quercitin is a strong antioxidant compound [50] that is antagonistic to the functions of the auxin class of phytohormone [12]. Inhibiting the role of auxin may lead to broad changes in resource allocation to the scale of plant growth habit. Such influence may lead to changes in plant scale resource allocation and concomitant traits associated with resource limitation. It is tempting to suggest that the expression of quercitin and its antagonistic effects on a major hormonal signalling pathway may serve to explain at least some of the variation in growth habit and water use observed across the genus [35, 42, 49, 57]. Further investigations of this gene expression pattern and spatio-temporal assessments of transcript and compound abundance may serve to unravel this fundamental process associated with ecotypic variation within Eucalyptus.

Primary photosynthesis and cellular metabolism

A range of pathways and gene analogues governing the cycling of primary photosynthates and respiratory substrates were observed across the network. Combined, these indicate a broad reconfiguration of primary photosynthetic reactions in response to the water deficit treatment. For example, UDP-D galactose-4-epimerase-1 governs the inter-conversion of galactose to glucose primarily for input into glycolysis. This gene shows differential regulation in response to the treatment (Fig. 3ci). UDP-D galactose-4-epimerase-1, producing glucose-1-phosphate, may also provide substrate for alternative pathways, most notably to the synthesis of inositols via inositol 3 phosphate synthase and downstream pathways previously linked to stress tolerance – including that within the Eucalyptus genus [2, 4, 32]. Cycling of photoassimilates and subsequent removal from the primary photosynthetic reactions of the cell is thought to alleviate sugar-mediated repression of photosynthesis in plants (e.g. [16, 18, 46]) thus avoiding, at least in part, photo-oxidative damage. A major component of this is the galactosylation of sugars for symplastic loading mechanisms characteristic of many tree species [52] including that of Eucalyptus [36, 37, 45]. Gene analogues governing the galactosylation of carbohydrates to form di-, tri- and tetra-saccharides such as melibiose, raffinose and stachyose is observed (Fig. 3ci) with potential implications for rates of carbon export from leaves.

The interconnectedness of pathways within primary metabolism is a major contributor to the capacity for homeostatic regulation of the system to balance inputs such as water and light. It is therefore not surprising that response mechanisms such as redox balance, photoprotection and changes in the cycling of primary photosynthates dominate the response to the water deficit treatment and share commonalities with previous investigations in response to alternative stress types (e.g. [17, 54]). An example of this is the expression of phosphogluconolactonase (Fig 3ci) closely linked with the pentose phosphate pathway. The pentose phosphate pathway may replace glycolysis under suboptimal conditions producing NADPH and a range of 4 and 5 carbon sugars. Products of the pentose phosphate pathway may provide substrate for several biochemical processes such as input back into glycolysis or the shikimic acid pathways for subsequent synthesis of aromatic compounds. Within primary metabolism, the necessity for tight regulation of the network means that subtle differences in transcript abundance may have large consequences on the observed phenotype. It is therefore important to remember that transcript abundance – even that of differentially expressed genes among contrasting ecotypes – is not a sole indicator of the significance of the response.

Secondary metabolism

Previous investigations into secondary metabolites in Eucalyptus have illustrated a broad diversity and range of concentrations among leaf oils [23, 24], waxes [25], phenolics [7] and tannins [29, 44] with several efforts to characterise genetic associations that underpin this diversity [19, 43]. In this study, both common and differential response in the expression of candidate genes governing the synthesis of secondary metabolites was observed between the two species. For example, tyrosine transaminase (Fig. 3cvi) is differentially expressed in response to treatment. Tyrosine is a phenolic amino acid thus may act as a precursor to the synthesis of other phenolic compounds or alkaloids [53]. Linkages between the synthesis of tyrosine and the synthesis of compounds known to have homeostatic redox functions such as the vitamin E in close proximity to gene analogues for plastiquinol biosynthesis indicate that redox balance is a major factor in plant acclimation strategies to the effects of imbalances between the availability of water and excessive photon flux. The diversity of secondary metabolism in the genus is reflected within the network compiled and serves as a valuable sequence resource to investigate differentially expressed chemical traits within the genus.

This study reports a novel experimental approach that achieves a comparative investigation of transcriptome responses two Eucalyptus species of contrasting ecotype to medium term water deficit. The range of transcript responses among the species highlights the mechanisms employed to maintain cellular function under a contrasting environmental conditions. These patterns, in several cases having the potential to influence broader structural and chemical processes governing adaptation to water deficit, are likely to be reflected in gene expression among a broader range of species within the genus and thus are candidates underpinning inter-species variation. This study highlights the strength of a transcriptomic approach to characterise plant responses in the context of metabolic network scale expression profiles and follow up experiments on more species and a range of time points will enable discovery of the genes that regulate species specific responses to water deficit as well as to the gene sets that define those responses. Compilation of such repositories for transcriptome data will undoubtedly serve as a valuable resource for further investigation at a range of scales into the responses employed by Eucalyptus to tolerate such a wide scope of environmental conditions.

Conclusion

This study describes and annotates an enhanced transcriptional response to water deficit in key pathways and associated genes in a ‘dry’ (xeric) climate eucalypt species adapted to water limited environments compared to a species originating from less arid conditions (mesic). The observed transcriptional response was of a greater magnitude and diversity in stem tissue. Of the differentially expressed genes, a significant proportion were of unknown function, indicating both unique and poorly understood mechanisms underlying the success of members in this genus in water limiting environments. Genes with known function were found in pathways linked to photoprotection and redox balance, phytohormone and long-distance signalling, primary photosynthesis and cellular metabolism. This suggests that the ‘dry’ climate eucalypt species has tightly regulated mechanisms across the metabolic network for broad scale modification of their biochemistry to reduce water loss while coping with the resultant increases in heat and light energy associated with the inability regulate this via modified physiology. Identification of these pathways and associated genes provides the opportunity to further investigate and understand the mechanisms and genetic variation linked to this important environmental response and assist in genomic efforts linked to managing native populations or tree improvement programs.

Methods

Plant material and experimental design

Seeds of E. cladocalyx were obtained from the Australian Tree Seed Centre at the Commonwealth Scientific and Industrial Research Organisation (CSIRO) and were germinated and raised in a naturally illuminated glasshouse at the University of Melbourne, Creswick (l37.43° S, 143.89° E). Seedlings were germinated in seed raising mix (Osmocote™>). Day/night air temperature was maintained at 20/12 °C and relative humidity at 60/70%. Irradiance as measured at the canopy level was approximately 600 μmol m−2 s−1 (PPFD) for a 13 h photoperiod (6:00–17:00). After 8 weeks, 60 plants were re-potted into 9 L containers containing a mixed peat/sand 50/50 v/v fertilised with 4 g L−1 slow release fertiliser Nutricote-100, 13/13/13 N/P/K with oligo elements, and watered six times daily by drip irrigation to field capacity. Seedlings were grown for 20 weeks prior to experimentation. For E. globulus, 60 ramets of clone X46 were sourced Narromine Transplants (Narromine, Australia) and potted into 9 L containers. Seedlings were grown following the same protocols and conditions as E. cladocalyx and experimentation took place after 14 weeks of growth in the glasshouse.

For experimentation, 24 E. cladocalyx and 36 E. globulus seedlings of similar height and leaf area were selected and segregated into two treatments: well watered (WW) and severe water deficit (SS). Treatments were determined gravimetrically based on calculated water use of control plants [30]. SS plants were given 45% of water used by WW to exhaust the reservoir of water in the pot by week 8 of the treatment period.

Leaf gas exchange measurements

Patterns in leaf gas exchange across the light period were determined on six plants of each species in each treatment using a LICOR 6400 Infra-red gas analyser (LICOR, Lincoln, Nebraska, USA) with a blue-red diode illuminator for leaf-level photosynthesis from 9 am (two hours after the photoperiod began) to 5 pm. A single leaf was selected on each plant and measured sequentially with at least 10 independent measurements taken on each leaf on each day. Fully expanded leaves were chosen for measurement with a different leaf used in each diurnal cycle. Light conditions within the LICOR chamber were set to tracking mode to approximate the growth chamber conditions. CO2 mole fraction in the reference air stream was set to 400 μmol mol−1 and temperature and relative humidity in the measuring chamber was maintained within the LICOR chamber to approximate ambient conditions. Net CO2 assimilation rate (A, μmol m−2 s−1) as well as stomatal conductance to water vapour (gs, mmol m−2 s−1) were recorded and used to compute the ratio of intercellular to atmospheric CO2 concentration ci/ca.

Tissue collection

Plant tissues were harvested between 12 and 2 pm on a sunny day at the end of week 8 of the treatment period and placed immediately into liquid nitrogen before being stored at −80 °C prior to RNA extraction. E. globulus (EG) tissue sourced was pooled from four separate individuals from the WW and SS trails including Apical Tip (AT), comprising of tissue above the first discernable node of the apical stem, Fully Expanded Leaf (FL), large leaves located in the top 1/3 of the plant that had full exposure to light, Primary Stem (ST), node 2, 3 and 4 of the apical stem and exhibiting a juvenile state (squared shape), Secondary Xylem (XY), xylem tissue scrapped with significant pressure from the surface of the stem using single edge razor blade after phloem removal in stems >2 cm in diameter, and Secondary phloem (PH), the secondary phloem tissue removed from stems >2 cm in diameter when sampling the XY tissue. In total there were 10 different tissue/treatment tissue samples for this species (Additional file 8: Figure S3). E. cladocalyx (EC) tissue was sourced from one individual in the WW and SS trails as individuals were not clonally related and included the following tissues; AT, ST and FL, as per E. globulus description, except in the case of ST where the first node of a number of lateral shoots was sourced (E. cladocalyx did not show strong apical dominance and was highly branched), as well as Secondary Stem (S2), where all xylem (developing and developed) and phloem tissue from a stem >1 cm in diameter was sourced due to difficulties of separating phloem from the xylem in the SS trial. In total, eight different tissue/treatment tissue samples were present for this species (Additional file 8: Figure S3).

Information from the plant metabolic network database Plantcyc (http://www.plantcyc.org/ ) and the pathway/genome database Biocyc (http://websvc.biocyc.org/) as well as the Gydle mapping as part of this study were used to visualise gene pathway networks using Cytoscape (V3.1.1, http://www.cytoscape.org/). To create the cytoscape graph, data was converted to triples data, loaded into a triplestore (OpenLink Virtuoso) and then linked via the TAIR identifiers. Sparql queries where then used to connect identified genes from the Gydle mapping to pathways. For network display, identified genes with FDR value of <0.2 for the DroughtW_FDR (response to water deficit) and SpeciesEG_DroughtW_FDR (species specific response to water deficit), were used (Additional file 2: Dataset S2). This modest FDR was chosen to facilitate the development of a suitably representative gene pathway network to support elucidation of responses. To assist with interpretation, additional data on the identified genes was visualised through attributes including graph line colour, style and width as well as gene pie charts found at each node. Details of these attributes can be found in the Additional file 7. Gene ontology functional classification was undertaken using Panther [38] and Singular Enrichment Analysis was undertaken using AgriGO [11] on genes with and FDR < 0.01.

Validation of transcripts

Validation of differential gene expression was through an independent (unpublished) study of leaf gene expression responses in the same two eucalypt species, treated under a similar medium term water stress regime. Protocols for tissue collection, RNA extraction, sequencing, bioinformatics and analysis were as described above with the exception that library construction targeted total RNA rather than mRNA as in this current study. The analysis compared both stem and leaf tissue from this study to leaf tissue from the validation study where logFC and FDR values where used to compare expression patterns across common genes found in both studies that showed a response to the drought treatment.

Declarations

Acknowledgements

The authors would like to thank Dr Luibov Volkova, Mr Raymond Dempsey and Mr Julio Najera for their assistance with gas exchange and sample collection.

Funding

This work was supported by the ARC Future Fellowship program (FT 120100200, Merchant), and a University of Melbourne Early Career Researcher (ECR) grant (Spokevicius).

Availability of data and materials

Data supporting this research paper can be found in the additional files section.

Authors’ contributions

AS and AM conceived the study, undertook the glasshouse trails, and prepared the manuscript, AM oversaw the collection and analysis of the gas exchange data, AS and JT undertook the RNA extraction, mSEQ library preparation and sequencing, JT, PR and MN undertook the sequence alignment, annotation, analysis of the sequence data, creation of pathway network and some manuscript preparation. CM and AM undertook the validation study. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Seeds of Eucalyptus cladocalyx were sourced from Australian Tree Seed Centre at the Commonwealth Scientific and Industrial Research Organisation (CSIRO) while Eucalyptus globulus seedlings were sourced Narromine Transplants (Narromine, Australia). Both orgnisations comply with all appropriate laws and guidelines when sourcing plant material. No voucher specimens were made.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Additional file 3: Figure S1.Plot of log fold change in gene expression for leaf and stem tissues between water deficit (SS) treatment or control (WW) treatment as detected in this experiment (X-axis) versus log fold change in gene expression in leaf tissues between water deficit (SS) treatment or control (WW) treatment (Y-axis) as detected in a follow up water deficit experiment on seven eucalypt species (N = 56). Grey (+) symbols represent all genes while red (o) represent genes validated at an FDR of < 0.01 in both experiments. It is clear that validated genes also show responses in the same direction between experiments. (PNG 25 kb)

Additional file 5: Figure S2.Small unlinked clusters depicting genes with an FDR of <0.2 with known links to pathways in plant metabolic networks. Small unlinked clusters have been arranged according to whether they show a species specific water deficit response (solid lines), with or without a water deficit only treatment response (dashed lines), (a), or water deficit only treatment response (dashed lines) (b). Clusters have been arranged further into those with a higher expression in response to the water deficit treatment (SS) (i), a differential response to treatment where some genes in a cluster show higher expression in response to water deficit treatment (SS) and others to the well watered control (WW) (ii), or higher expression in response to the well watered control only (WW) (iii). (PNG 2743 kb)

Additional file 6: Dataset S4.The full cyctoscape file to allow for more detailed examination of genes/pathways and how they have responded to the water deficit treatment (SS) and well watered control (WW). To view file, please download the Cytoscape V3.3.0 program free from the following website: http://www.cytoscape.org/ . Interpretation instructions can be found in Additional file 7. Please note the node pie colours are not available in this file and readers should refer to Fig. 3 and Additional file 5: Figure S2. (CYS 639 kb)

Additional file 7:Description of details relating to line colour, style and width as well as node pie colour to assist with the interpretation of Fig. 3 and Additional file 5: Figure S2 and Additional file 6: Dataset S4. (DOCX 15 kb)