Analysis of the transcriptome of the needles and bark of Pinus radiata induced by bark stripping and methyl jasmonate

Background Plants are attacked by diverse insect and mammalian herbivores and respond with different physical and chemical defences. Transcriptional changes underlie these phenotypic changes. Simulated herbivory has been used to study the transcriptional and other early regulation events of these plant responses. In this study, constitutive and induced transcriptional responses to artificial bark stripping are compared in the needles and the bark of Pinus radiata to the responses from application of the plant stressor, methyl jasmonate. The time progression of the responses was assessed over a 4-week period. Results Of the 6312 unique transcripts studied, 86.6% were differentially expressed between the needles and the bark prior to treatment. The most abundant constitutive transcripts were related to defence and photosynthesis and their expression did not differ between the needles and the bark. While no differential expression of transcripts were detected in the needles following bark stripping, in the bark this treatment caused an up-regulation and down-regulation of genes associated with primary and secondary metabolism. Methyl jasmonate treatment caused differential expression of transcripts in both the bark and the needles, with individual genes related to primary metabolism more responsive than those associated with secondary metabolism. The up-regulation of genes related to sugar break-down and the repression of genes related with photosynthesis, following both treatments was consistent with the strong down-regulation of sugars that has been observed in the same population. Relative to the control, the treatments caused a differential expression of genes involved in signalling, photosynthesis, carbohydrate and lipid metabolism as well as defence and water stress. However, non-overlapping transcripts were detected between the needles and the bark, between treatments and at different times of assessment. Methyl jasmonate induced more transcriptional responses in the bark than bark stripping, although the peak of expression following both treatments was detected 7 days post treatment application. The effects of bark stripping were localised, and no systemic changes were detected in the needles. Conclusion There are constitutive and induced differences in the needle and bark transcriptome of Pinus radiata. Some expression responses to bark stripping may differ from other biotic and abiotic stresses, which contributes to the understanding of plant molecular responses to diverse stresses. Whether the gene expression changes are heritable and how they differ between resistant and susceptible families identified in earlier studies needs further investigation. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08231-8.


Introduction
Plants have evolved a variety of constitutive and inducible defences to resist and tolerate herbivory. An assessment of the genetic mechanisms that influence these defences will enhance our understanding of their evolution [1]. Although structural changes in DNA are the major source of genetic variation [2,3], the phenotypic outcomes of several traits can be linked to gene expression [4][5][6][7][8]. However, the genes and genetic pathways that underlie most phenotypes are still unknown [2]. To date, most gene expression studies have focussed on identifying transcripts (different RNA products a single gene) or genes showing differential expression, or pathways associated with a phenotype (case/control) or condition (treated/untreated). In conifers, for example, transcript abundance has been examined with respect to biotic and abiotic environmental factors such as herbivory [9][10][11], pathogens [12], artificial wounding [13], drought [14], light intensity [15], seasonal changes [16], chemical stressors like methyl jasmonate [17], as well as associated phenotypic traits such as resistance and chemical composition [9,10]. Studies in conifer and non-conifer species that have simultaneously compared the expression from different stressors, such as mechanical wounding and methyl jasmonate, indicate both overlapping and non-overlapping gene expression and suggest that molecular mechanisms associated with varying stressors may differ [18][19][20].
In conifer-herbivory studies, most gene expression studies have focused on understanding induced defence responses, with a premise that these may be more important than constitutive defences as they are metabolically cost effective and expressed only when required [21,22]. Global transcriptome responses have been studied in both needles and bark, monitoring the expression of a wide range of genes related to the biosynthesis of primary and secondary compounds, and structural components [13,[23][24][25][26][27][28]. Most of these genes are expressed at basal levels in plants but some are only expressed in the presence of an appropriate stimulus. Some of the genes significantly respond to herbivory cues, by increasing or reducing their expression either locally at the site of the perceived effect or systemically throughout the plant [23,29,30]. Studies also show a high overlap in the genes that are differentially expressed when plants are subjected to different biotic and abiotic stresses [31,32]. However, the genes that show differential expression differ within and between target plant species [10,26], between plant tissues [23,33], as well as between biotic agents [34] and applied treatments [35]. Intra-specific differences in the timing of transcript expression have also been observed, where plants may respond to injury within hours or days, with short, or long, lasting effects [17,23,25,33]. Plant responses to different classes of herbivores may differ due to differences in herbivore oral secretions or mode of feeding and the amount of plant tissue damage [34,36,37]. While available conifer studies have documented changes in gene expression in response to insect herbivory [13,32], there are no studies from the perspective of mammalian herbivory, and none that link changes in gene expression to changing chemistry. Mammalian bark herbivory is fundamentally different from insect herbivory in the mode of feeding [22] and possibly the oral secretions. This particularly applies to mammalian bark stripping, which is of increasing concern to managers of conifer forests world-wide, including Pinus radiata plantations in Australia [38][39][40].
Pinus radiata is native to California [41], but is now a major plantation species in Australia (ABARES 2018) where it is subject to bark stripping, mainly by native marsupials (wallabies and kangaroos) [42]. The bark is stripped from the base of the trees during the early stages of growth [43][44][45], reducing tree growth rate, distorting stems and, in extreme cases, causing death [38,42]. The levels of bark stripping within plantations may be highly variable and progeny trials have shown a genetic, physical and chemical basis to this variation [42,46,47]. Further, chemical profiling in P. radiata shows that needles and bark respond differently to bark stripping and other forms of real and simulated herbivory, mostly by increasing levels of secondary compounds, especially terpenes and phenolics [48,49], and reducing levels of sugars and fatty acids [46,50]. This suggests changes in the expression of underlying genes that subsequently transforms the chemical phenotype. Indeed, the differences in timing of the induced changes in terpenes, phenolics and sugars [50][51][52] suggest corresponding differences in the expression of the underlying genes. However, while transcriptomic changes have been studied in P. radiata associated with ontogeny, wood formation [53][54][55] and fungal infections [56], those underlying the induced chemical changes to bark stripping have not been characterised.
The present study aims to quantify and compare the transcriptome changes that occur in response to artificial bark stripping of P. radiata and whole plant stress induced by application of the chemical stressor, methyl jasmonate. The longer-term goal is to identify genes that specifically mediate the previously shown induced Keywords: Transcriptome, Chemical phenotypes, Bark, Needles, Pinus radiata chemical responses to bark stripping in P. radiata, which may help develop strategies to reduce bark stripping. The specific aims of the study are to: 1) characterise and compare the constitutive transcriptome of P. radiata needles and bark; 2) identify genes which are differentially expressed following artificial bark stripping (aimed at mimicking mammalian bark stripping); and 3) identify genes which are differentially expressed following whole plant application of methyl jasmonate and compare these induced responses with those of bark stripping. The results are discussed in view of the holistic chemistry that has been characterised on the same individuals with the same treatments [50].

Experimental design
In 2015, 6-month-old seedlings from 18 full-sib families (each with 4 seedlings; total number of seedlings = 72) of P. radiata (D. Don) originating from the Radiata Pine Breeding Company deployment population, were obtained from a commercial nursery. Seedlings were transferred into 145 mm × 220 mm pots containing 4 L of basic potting mix (composted pine bark 80% by volume, coarse sand 20%, lime 3 kg/m 3 and dolomite 3 kg/ m 3 ) and raised outdoors in a common fenced area (to protect against animal damage) at the University of Tasmania, Hobart. At 2 years of age, plants were moved to a shade house and an experimental design established by randomly allocating the 18 families to three treatment groups (methyl jasmonate [MJ], artificial bark strippingstrip [strip] and control), each with 6 families. The three treatment groups were arranged in a randomized block design of 3 blocks, each block comprised a treatment plot of two families, with the treatment plots separated within each block to minimise any interference among treatments. Each family was represented by four plants arranged linearly, and randomly allocated to four sampling times (T0-T21). T0 represents the time immediately before treatment applications. T7, T14 and T21 represent respective sampling times at 7, 14 and 21 days after treatment (MJ and strip) application. All T0 seedlings (n = 18), irrespective of group allocation, were not treated and were used to compare the constitutive transcriptome of the needles and bark (i.e. plant parts). Additionally, all seedlings allocated to the control were not treated throughout the experimental period. One seedling from each family in the control and treated groups was destructively sampled at each sampling time to estimate differential expression (n = 18; Table 1). For each plant part, comparisons were made between the control (n = 6) and methyl jasmonate (MJ, n = 6) and between the control (n = 6) and bark stripping (strip, n = 6) treatments at each sampling time (T7, T14, T21) ( Table 1). Methyl jasmonate (MJ) was applied in a 25 mM solution by spraying the whole plant with a fine mist from a hand sprayer until 'just before run-off ' . The treated seedlings were sprayed in a well-ventilated area away from untreated seedlings to avoid cross contamination [57]. For bark stripping (strip), 18 plants were artificially stripped by removing a 30 cm vertical strip of bark, beginning 2 cm from the ground and covering 50% of the stem circumference, which is the average upper threshold of browsing observed in natural field conditions. Up to 20 young needles were randomly collected per seedling from different parts of the crown. The bark was sampled from different points of the stem, above and besides the area where the bark stripping treatment was applied, carefully avoiding the wood, following Nantongo et al. [50]. Individual samples were kept separate providing 144 samples for sequencing (2 plant Table 1 The treatments, sample size and pairwise comparisons that were made for each time and for the two treatments -bark stripping (strip) and methyl jasmonate (MJ). The seedlings of each family were grown in a line-plot and one was chosen at random for destructive harvesting at each time (T7 to T21). At T0, the sampled seedlings were destructively harvested just before treatment applications. At 7 (T7), 14 (T14) and 21 (T21) days after treatment, one seedling from each family (total number of seedlings per sampling time = 18, equivalent to the number of families and n = 6 are seedlings selected from each treatment) was destructively harvested  [58]. Quality trimming and filtering of data was performed using Trimmomatic v 0.39 [59]. On average, 99.9% of the sequences were retained at phred33 [60]. A de novo assembly of the pooled transcriptome was attempted using TRINITY v2.9.0 using default parameters [61], however due to the excessive computation requirements, it could not be completed with the available resources in the required timeframe. Accordingly, the filtered reads were aligned to the P. radiata reference transcriptome that is harboured at Scion (the New Zealand Forest Research Institute trading as Scion, Rotorua New Zealand) [54] with SALMON v0.14.1 using default parameters [62]. This reference transcriptome (www. ncbi. nlm. nih. gov/ biopr oject/ 482145) was assembled from a range of P. radiata genotypes and tissue types that were collected at different developmental and temporal stages. Most of the samples were from healthy seedlings under normal growth conditions but also included some pathogen infected seedlings [54]. The reference transcriptome has a total of 279,510 unique transcripts.

Differential transcripts expression analysis
Statistical analysis of differential expression was performed using the edgeR v3.24.3 package in R (v3.6.0) [63] using default parameters [64], except for the cut-off false discovery rate (FDR) in treated samples that was modified as described below. EdgeR uses the Poisson distribution model to examine differential expression of replicated count data, which makes it simpler than methods that use other statistical distributions [65]. Transcripts were first filtered retaining only those with a minimum expression change of 2 fold and with a minimum of 100 counts per million of a single transcript in at least two part x treatment x time groups. To adjust for library sizes and skewed expression of transcripts, the estimated abundance values were normalized using the trimmed mean of M-values normalization method included in edgeR. To detect differential transcript expression between the needles and the bark, the samples taken at T0 were used as these comprised a single plant from each of the 18 families (as treatments were not applied at this stage) and an FDR value of 0.05 was used. However, to establish transcript expression after treatment, instead of using an FDR of 0.05, a more conservative sample-specific approach was used [66], where transcript expression was initially compared between the samples collected from the control plants (n = 6), MJ-allocated (n = 6) or strip-allocated (n = 6) groups at T0 (before treatment) to check the inherent (potentially random) differences between sample groups. The p-values at which no differential expression was detected between these groups was set as the FDR for downstream pairwise comparisons. Accordingly, the p-value for detecting differentially expressed transcripts (DET) in the treated needles following both MJ and bark stripping was set at 1.0 × 10 − 11 . A p-value of 1.0 × 10 − 18 was set to detect DET in MJ treated bark and 1.0 × 10 − 10 to detect DET in the bark stripped samples. Twelve pairwise comparisons were performed. An upset diagram was generated using the UpSetR function in R to summarise the transcripts that were identified as significantly differentially expressed across different comparisons.
Principal component and unsupervised cluster analyses were performed to detect the dominant, relative expression patterns across the needles, bark and treatments. Following Ralph et al. [13], a subset of 500 transcripts with the highest variability and highest expression across the 143 libraries were selected in edgeR for this analysis. Principal components analysis (PCA), using Facto-MinerR version 1.41 [67] was based on the correlation matrix among all identified transcripts. Clustering and heat maps were generated using the heatmap.2 function from the gplots package in R, with a matrix of Euclidean distances from the log2 counts of normalised transcripts.

Sequence similarity search
For sequence similarity search and functional analysis of differentially expressed transcripts (DETs) the transcripts were blasted against the nucleotide BLAST database using BLASTn (https:// blast. ncbi. nlm. nih. gov/ Blast. cgi). BLAST analysis revealed that P. radiata transcripts were most similar to those predicted from genome sequences of P. taeda (BLASTn with e-value < 0.0001). Other species, mostly P. sylvestris, P. monticola, Picea stichensis and Pseudotsuga menziesii, showed high similarity with the P. radiata transcripts. Annotations of selected transcripts were done by comparing P. radiata transcripts to the sequences in the SwissProt database of annotated genes [68] using cut-off values ≤ 1. To gain clear patterns of the responses, only transcripts associated with genes of known function were included. However, there were many uncharacterised transcripts and proteins of unknown functions.

GO classification
Gene ontology (GO) classification was undertaken to understand the biological process, cellular component and molecular function categories represented in the genes exhibiting differential expression. These assignments were done for selected transcripts identified above using protein analysis through evolutionary relationships (PANTHER) version 14.1 [69]. This was first undertaken using transcripts that were differentially up-regulated in the needles over the bark and vice versa, with the aim of understanding the constitutive differences of the GO processes between the transcriptome of the needles and the bark. Secondly, the GO classification was performed on selected T1 transcripts to understand the differences in the up-regulated and down-regulated transcripts after treatment, as well as differences in the induced transcriptome of the strip and MJ treated samples. Due to the limited annotation resources available for conifers, gene family annotations were obtained using genomes of 10 species: Arabidopsis thaliana, Citrus sinensis, Cucumis sativus, Oryza sativa, Populus trichocarpa, Prunus persica, Saccharomyces cerevisiae, Theobroma cacao, Vitis vinifera and Zea mays. GO term classification was done for the top differentially expressed transcripts in the different conditions (time × treatment × part).

Results
The Pinus radiata reference transcriptome and read mapping RNA-seq of P. radiata generated a total of 2860 million 100-bp PE reads with a minimum of 20 million reads from each of the 143 samples. 87.6% of the reference transcriptome was represented among the study transcripts.
However, after the filtration criteria described above, only 6312 unique transcripts (2.6% of the reference transcriptome) were retained as the expression of the other transcripts was too low. The analysis was constrained to individual transcripts, which may not be unigenes.

Differential expression of the transcriptome
The overall relationships between the transcriptome from the different samples were visualised using a principal component analysis (PCA) plot ( Fig. 1) and the unsupervised hierarchical clustering (Fig. 2) of the top 500 variable transcripts in the transcriptome. Both figures show that the major differences in expression were due to plant parts (differences along the x-axis of Fig. 1 and the top x-axis of Fig. 2). Within plant parts, we noted genes that were: (i) up-regulated in the needles relative to the bark and generally non-responsive to treatment; (ii) up-regulated in the bark relative to the needles and generally non-responsive to treatment; (iii) up-regulated in either the needles or the bark and responsive to treatment; and (iv) not differentially expressed between the needles and the bark but responded to treatment by up-or down-regulation.

Differences in the constitutive needle and bark transcriptome
Of all 6312 transcripts considered for analysis, 5 transcripts were detected only in the needles and 13 transcripts were detected only in the bark. Most of these part-specific transcripts were uncharacterised ( Table 2). Gene level annotation of the top 10 transcripts expressed in each plant part are listed in Table 3 (superscript refers to ID number in Table 3). The type 2 light-harvesting chlorophyll a/b-binding polypeptide [1] that is possibly involved in photosynthesis, was the most expressed gene in both the needles and the bark and was represented by different copies of transcripts (isoforms). The needles had other photosynthesis-related genes expressed such as ribulose bisphosphate carboxylase/oxygenase (RuBisCO) [12] and PSI-D1 precursor [17] possibly due to its major role in photosynthesis. Genes related to secondary metabolism were also detected among these top 10 genes, suggesting that constitutive defence is important in P. radiata. These included dehydrin [2] , metallothionein [3] , chalcone synthase [4] , defensin [5] and pathogenesis-related proteins [8] and were represented by more transcripts in the bark than in the needles but their relative expression was not statistically significantly different between the needles and the bark.
At T0, 5469 out of the 6312 transcripts (86.6%) were differentially expressed between the needles and the bark. Of these, 3123 were up-regulated in the bark compared to the needles, while 2346 transcripts were up-regulated in the needles. The top 10 most strongly up-regulated transcripts in each of the bark and needles are shown in Table 4 (superscripts are identifiers to help locate the needle (N) or bark (B) transcripts in the ID column of the table). Besides the general function genes and those related with photosynthesis, there was an up-regulation of genes related to terpene [B9] and lipids biosynthesis [B7] in the bark and those related to sugars [N4] and phenolics biosynthesis [N1] in the needles. Of note is the up-regulation of genes involved in sugar transport in both the needles [N3] and the bark [B2] , but these are different genes.
To assess the overall constitutive functional differences in transcripts differentially upregulated in the needles and the bark, the GO annotation of the top 100 differentially upregulated genes in both plant parts was obtained. There were quantitative differences for all the molecular but not biological or cellular GO categories. In the molecular GO category, a greater proportion of the top upregulated genes in the needles were ascribed to catalytic activity in the needles than in the bark (Fig. 3).

Overall transcript expression in the needles and the bark after treatment
After treatment, considering all time points, a total of 1479 (23.4%) transcripts were differentially expressed at one time or another. More transcripts responded to treatment in the needles than in the bark and more transcripts were up-regulated than down-regulated (Fig. 4). For both treatments, most differential expression was detected 7 days (T7) after treatment and declined thereafter, although differential expressed transcripts were still evident in both treatments 21 days later (Fig. 4). MJ was applied to both bark and needles and caused more transcript expression than bark stripping in both the needles and the bark (Fig. 4). Indeed, no differential expression of transcripts was detected in the needles following bark stripping. Of the transcripts that were differentially expressed between the bark and needles at T0, only 20% and 1% of those respectively responded following either of the treatments in the bark and needles suggesting that the transcripts that did not differ constitutively (i.e. at T0) between the needles and the bark were more responsive to treatment. One uncharacterised transcript (NZPradTrx091980_C05) that was not present in the transcriptome of untreated samples was present after treatment. One isoform of ribulose bisphosphate carboxylase preprotein (NZPradTrx098233_C06) that is involved in photosynthesis was present before treatment but was missing in all the samples in the bark and the needles after treatment, including the untreated control samples.
Annotations of the top ten genes that were up-regulated or down-regulated for each condition (time × treatment × part) are presented in Table 5. Based on these genes, various functions were detected, indicating that multiple genes are involved in coordinating plant responses to stress. Most of the genes were up-regulated, for example genes associated with primary metabolism, secondary metabolism, digestive inhibitors, pathogenesis-related (PR) protein families, genes involved with physical strengthening of the cell-wall, transcription factors, phytohormones and signalling molecules as well as molecules involved in broad biotic and abiotic stress responses and broad function genes. In contrast, the general catalysts as well as molecules involved in transcription were down-regulated. A subset (968 out of 1479 = 64.7%) of the differentially expressed transcriptome studied was differentially expressed in only one treatment (strip or MJ) (Fig. 5, Table 5). Similarly, nonoverlapping differentially expressed transcripts, occurring in only one condition, were detected at different times in the needles and bark (Fig. 5, Table 5).

Gene expression after MJ treatment
A stronger response to the MJ treatment was detected in the needles than the bark, where 2206 versus 683 out of 6312 transcripts studied were differentially expressed, respectively (Fig. 4). Annotations of the non-overlapping, differentially expressed transcripts showed that MJ caused the unique differential expression of more genes that are directly involved in the metabolism of sugars, and time/part/treatment categories (columns) were clustered using Euclidean distance. The Z-score is calculated by subtracting the trimmed mean of the M-value of the individual from the grand mean of all the individuals and then dividing by the standard deviation. Trimmed Means of M values are estimated in edgeR by where highly expressed genes and those that have a large variation of expression are excluded, whereupon a weighted average of the subset of genes is used to calculate a normalization factor. Colouration; yellow = mean expression, blue = expression below the mean and red = expression above the mean. The categories on the x-axis were re-arranged based on similarity fatty acids and amino acids in both the bark and the needles compared with the bark stripping (Table 6).
Six transcripts were consistently differentially expressed from T7 -T21 ( Fig. 5) in the methyl jasmonate-induced transcriptome of the bark (B-MJ) and these were mostly up-regulated. Annotations of these transcripts showed that the genes were mostly involved in generating energy from various substrates, particularly glucose and fatty acids. In the needles treated with methyl jasmonate (N-MJ), 114 transcripts were consistently differentially expressed from T7 -T21 (Fig. 5). These genes were mostly directly associated with defence as well as chemical and physical structures, for example those involved in phenolic biosynthesis and structural components of the cell wall (Table 5).

Gene expression after bark stripping
Bark stripping did not cause any systemic response in the needles at any time point (Fig. 4). The strip induced bark transcriptome had, among the top genes, those involved in defence against pathogens, such as chitinases [U17] , PR10 [U39] and defensins [U18] . Bark stripping also caused differential expression of water-stress responsive genes [U12,U39] as well as genes related to replacement of tissues [U34] (Table 6). The difference in the representation of genes is likely related to the kind of damage incurred by the two stressors.
Both stressors caused differential expression of genes related to secondary metabolism (Table 5), including metabolism of monoterpenes (e.g. geranyl diphosphate synthase), phenolics (e.g. laccases) and alkaloids (e.g. phenylalanine ammonia-lyase). The differential expression of genes associated with lignification of cell walls were also identified for both treatments in the needles and the bark, emphasising the role of cell wall physical properties in stress responses. For some genes, the same gene was represented by different isomorphs in the different conditions such as geranyl diphosphate synthase in B-strip and N-MJ treatment/part combinations shown in Table 5. Only 6 differentially expressed genes were consistently differentially expressed following both treatments across all times and plant parts, except that no differential expression occurred in the needles following the strip treatment. Annotations of these transcripts mostly showed genes related to amino acid synthesis.

Time progression of genes
Not only did the treatments differ in the magnitude of their general response through time (Figs. 1, 4 and 5), but the pattern of response of individual genes differed between treatments. For the top ten expressed transcripts in the constitutive transcriptome (assessed at T0) of the bark and the needles (ID numbers 1 to 10 in Table 3), Fig. 6 shows the time progression of differential expression following stripping and methyl jasmonate application. There was a tendency for genes to be up-regulated or down-regulated following both treatments. Of the three genes (dehydrin, light-harvesting chlorophyll a/bbinding polypeptide and metallothionein) that showed marked down-regulation, only dehydrin showed   Table 5 Top 10 genes differentially expressed in each of the time periods from T7 to T21 in the bark (B) and needles (N) following bark stripping (S) or methyl jasmonate (MJ) treatment of two-year old Pinus radiata plants The Scion transcript code, predicted gene name and predicted functions of the known genes are indicated. Some genes were represented by more than one transcript (isoforms-different Scion P. radiata transcript codes that represent one gene in column 1) and multiple copies of an isoform as indicated by the numbers in the parentheses, for example +(2) = two copies of an isoforms relating to the gene were identified, where The superscript followingnumbers in the parentheses following the gene names represent the core function of the gene among the 11 broad categories listed in the table footnote. For example for the Peptide transporter PTR3-A-like, a the superscript a denotes that this gene was associated with primary metabolism (see footnote). However, it is recognised that some genes may fall in more than one category. Gene functions are mostly from UniProt [77] Scion transcript code

Functional classification of differentially expressed transcripts
To assess the overall effect of the treatments across different gene families and molecular processes, the GO terms were determined for the up-regulated and downregulated transcripts for each condition (time × treatment × plant part). There was an overall similarity in the GO terms for genes that were up-and down-regulated in the strip and methyl jasmonate treatments. For example, in the GO-molecular processes, differentially expressed genes were associated with catalytic activity both in the needles and the bark (Fig. 7, Supplementary Fig. 1). However, the proportion of the 100 top differentially expressed genes in the catalytic activity category varied markedly. For example in the bark, a great percentage of top down-regulated genes following bark stripping were in the catalytic activity category (72%) compared with the up-regulated genes (28%). Comparing GO terms for the top differentially expressed genes in the constitutive (needles versus bark) and induced transcriptome, indicated that some gene functions that were not strongly expressed in the constitutive state (T0) were notably up-regulated or downregulated after treatment, and this differential expression appears to be treatment specific (Fig. 7). For example, genes related to response to stimulus (GO:0050896), plasmodesma (GO:0009506) and cell junction (GO:0030054) were strongly up-regulated at T7 in the transcriptome of the bark stripped samples but not the methyl jasmonate samples. Accordingly, transcripts of many of the other GO categories were under expressed in the transcriptome of the bark stripped samples.

Discussion
We aimed to understand the differences in the constitutive needle and bark transcriptomes, the changes that occur following bark stripping and how they compare with those of methyl jasmonate that have been most commonly reported for conifer species [17,24,35,80]. While the results are based on a partial transcriptome, comparing the needle and bark transcriptome as assessed prior to treatment (T0) showed that there were minimal qualitative differences in terms of the transcripts found  in the plant parts. However, after treatment there was strong transcriptional response of the basal transcripts in both the needles and the bark, with the expression being different and with sometimes non-overlapping transcripts between plant parts, treatments and at each sampling timepoint. While the effects of methyl jasmonate have been previously reported in various pine species [17,24], this is the first study to illustrate transcriptional responses to bark stripping. The response to bark stripping was less than that to methyl jasmonate and was localised, as no systemic response extending to the needles was detected at any time point. Differences in responsiveness to both treatments were also detected between the classes of genes, where genes related to primary metabolism responded to treatments with a greater magnititude of up-regulation or down-regulation compared to genes associated with secondary metabolism. Among the genes that were homogeneously expressed between the bark and the needles were those related to basic life functions especially those related to primary and secondary metabolism. For example, ribulose bisphosphate carboxylase/oxygenase (RuBisCO) and a chlorophyll a/b binding protein were dominant both in the transcriptome of the needles and the bark. Similar observations were made in the needles of other P. radiata populations [81] and Pinus monticola [70], although these studies did not analyse how the transcriptomes change with treatment and the observations were limited to one plant part. Genes directly related to secondary metabolism, for example chalcone synthases, dehydrins and defensins, were among the basal genes, highlighting the importance of constitutive defences in P. radiata. Chalcone synthase has been identified in other conifers [82,83] and plays crucial role in phenolic biosynthesis [74]. Defensins have also been detected in various conifers where they inhibit the growth of a broad range of pathogens, including bacteria, fungi and viruses [75,76]. Dehydrins that represent a family of genes for drought tolerance have been detected in spruces and in other Pinaceae [72]. Metallothioneins that were strongly expressed both in the bark and the needles are important in protection against heavy metal toxicity [73] and have been documented mainly in Pseudotsuga menziesii [84,85]. They could reflect an adaptation to leached, heavy metal enriched soils in the coastal sites of California where P. radiata originates [86]. However, while the above genes are expressed at high amounts equally in the bark and needles, some transcripts were up-regulated in the needles or the bark. More up-regulation was detected in the bark, which contrasted with the higher expression  Table 3 and their differential expression in bark is shown following a bark strip and b methyl jasmonate treatments. The average change in expression was estimated at each time point by comparing the raw counts for the bark strip or methyl jasmonate induced transcripts and the transcripts from control treatments (mean of treatment -mean of control) for a specific time and were adjusted according to the differences in basal expression of the treatment groups at T0. T0 is before treatment applications and T7, T14 and T21 correspond to 7, 14 and 21 days after treatment application, respectively of transcripts in the needles than the bark reported in other P. radiata populations [81]. In both plant parts upregulated genes were predominantly related to the synthesis and transfer of macro-and micro-molecules, as well as transcription factors which are the key molecular switches orchestrating the regulation of plant responses to a variety of stresses.
After treatment with methyl jasmonate and bark stripping, there was an up-regulation and down-regulation of several genes involved in both primary and secondary metabolism both in the bark and needles, consistent with other studies that have characterised responses to other stressors in conifers [24,79]. The top genes that were up-or down-regulated in the present study overlap with those observed in similar studies with contrasting sources of stress in conifers [13,70,79,80,87], suggesting that changes in gene expression following stress are relatively conserved. Among the top expressed genes, results showed a down-regulation of hexokinases, granule-bound starch synthase and sodium-bile acid cotransporter as well as genes related with photosynthesis, suggesting reduction in sugar metabolism in the treated plants. However, cell wall invertase that mediates export of sucrose or enhanced import of hexoses at the site of damage was up-regulated in both methyl jasmonate and strip treated plants. Cell wall invertase (CWI) is an enzyme that cleaves sucrose, the major transport sugar in plants, irreversibly yielding glucose and fructose, which can be taken up by plant cells [78,88]. An increase in CWI should ideally lead to a reduction in sucrose, which is consistent with the drastic reduction in the amounts of sucrose that has been observed following methyl jasmonate and strip treatments in P. radiata. The up-regulation of CWI would also suggest an increase of glucose and fructose, but this was not the case as a strong reduction in the amounts of glucose and fructose was observed in treated samples [50]. This suggests that although fructose and glucose may be potentially enhanced by an increased break down of sucrose, their utilisation for energy and carbon skeletons for other organic compounds or for tissue recovery exceeds their production, supporting the concept that defence is costly in terms of energy [89]. Gould, Reglinski [90] detected a repression of photosynthesis in P. radiata as a response to stress that could lead to a reduction of sugars. Sugars have also been shown to function as signalling molecules, in a manner similar to hormones [88,91], but their down-regulation contrasts to the up-regulation of other signalling molecules. However, according to Eveland and Jackson [92] sugar signals are generated either by relative ratios to other metabolites, such as C:N, not necessarily carbohydrate concentration.
In addition to the sugar-related genes, the other primary metabolism genes that were responsive to the treatment included those genes related to fatty acid metabolism such as the medium-chain-fatty-acid-CoA ligase and UDP-rhamnose:rhamnosyltransferase that were up-regulated and those related to fatty acid hydrolysis, such as carboxylesterase, that were down-regulated. Observations on the same population showed a reduction in fatty acids following treatment, consistent with their potential use as precursors to the formation of secondary compounds [93]. Accumulating evidence has suggested lipids and lipid metabolites as important regulators of plant defence [94]. Genes related to amino acid synthesis were also among the top expressed genes. Increase in amino acid levels have been detected in plants under stress and is hypothesized to protect plant cells against dehydration [95,96]. Amino acid accumulation has been observed to be strongly related to abscisic acid signalling [95]. Molecules related to abscisic acid signalling were also strongly up-regulated similar with pathogenicity response in the Pinus pinaster -Fusarium circinatum pathosystem [97]. This study contributes to the body of literature demonstrating the crucial role of phytohormones in host defense response [98].
Genes related directly to secondary metabolism were not detected among the top differentially expressed genes following treatment although they are abundant in the constitutive transcriptome of both the needles and the bark, consistent with the observations in spruce [10]. However, the relatively weak transcriptional response to treatment of individual genes related to secondary metabolism in this study contrasts with other studies [13,17] and could possibly be due to the timing of the sampling, which was done 7 days after treatment application. In various studies, maximum expression of genes is shown to be attained within 5 days after treatment application [13,17]. On the same population, a weak response of terpenes and phenolics was observed following similar treatments [50], which probably suggests an inherently weak response of secondary compounds and associated genes to stress in P. radiata. Defence genes being strongly expressed in the constitutive but not in the induced transcriptome may suggest existence of trade-offs in induced gene expression [99], analogous to the trade-offs in constitutive versus induced chemical responses that have been detected in P. radiata [21]. Although alkaloids have not been well researched as important defence compounds in conifers, genes related to alkaloid biosynthesis such as RS-norcoclaurine 6-O-methyltransferase were among the top expressed genes but were down-regulated after treatment. There were also many proteins of unknown functions that were up-regulated or down-regulated at various time points, which potentially explains the many unknown chemical compounds that were quantified on the same plants.
Considerable overlap was observed between the methyl jasmonate and the strip induced transcriptome. However, results also indicate that bark stripping can induce transcripts that are not induced with methyl jasmonate and vice versa. Defence responses for bark stripping may differ from methyl jasmonate since bark stripping causes tissue and water loss at the injured sites, and damaged plants are also easily infected by pathogens through these wounds. In this case both defence and repair responses are required. Hence the dominant genes in the strip-induced transcriptome involved pathogenesisrelated (PR) genes and those related to fibre synthesis. The expression of PR genes could also be related to the historical relationship between P. radiata and various pathogens [100]. No systemic transcript responses were observed in the needles to bark stripping. Coupled with the chemical changes that were observed in the needles following bark stripping on the same population, for example the reduction of glucose and fructose at T7 and T14 [50], this observation suggests that some chemical stress responses, possibly those involving sugars, may not involve on-site gene expression changes and may result from passive reallocation of chemistry within the plant. For other compounds like terpenes, it has been indicated that passive changes normally occur only in the constitutive environment and that stress-induced changes in terpenes are entirely of a de novo nature [101].
A key finding from this study is that the main transcriptome change associated with either treatment was clearly earlier than the main chemical changes observed on the same population [50]. The maximum differential expression of the transcripts was observed 7 days after treatment whereas most chemical change were detected 14 and 21 days after treatment, consistent with a time-lag between gene and phenotypic expression. This discrepancy may be associated with trade-offs between gene expression and other cellular resources, including the nutritional quality of the plant [99]. One GO-term that was significantly enriched after treatment was response to stimuli and, consistently, genes related to signalling were among the top expressed genes. For example, 1-aminocyclopropane-1-carboxylate oxidase, which is related to production of ethylene; lanC-like protein 2-like for abscissic acid and Tify domain containing protein for jasmonates were strongly responsive. Ethylene is one of the major signalling molecules in plant defences in addition to others, such as jasmonic acid, salicylic acid and abscisic acid [102]. Ethylene can act synergistically or antagonistically with jasmonic acid in the regulation of both stress and developmental responses. The connection between these two signalling pathways has been demonstrated genetically to be the transcription factor for the ethylene response [103], that was also strongly expressed. This suggests that jasmonates, abscisic acid and ethylene are involved in induced responses of P. radiata under different stresses. The involvement of jasmonates and ethylene in induced defence responses has been shown in other pine species [20]. In other species, abscisic acid has been shown to be involved in defence responses and has been reported to play a negative role in the regulation of the major photosynthesis genetype 2 light-harvesting chlorophyll a/b-binding polypeptide [71] -which was reduced after treatment in this current study.

Conclusion
There are marked quantitative differences in the needle and bark transcriptome of Pinus radiata both in the constitutive and induced states. The transcriptome triggered by bark stripping substantially differed from methyl jasmonate triggered responses suggesting that some molecular aspects of bark stripping may differ from other biotic and abiotic responses. Gene annotation revealed that some of the differentially expressed transcripts have putative functions in plant defence signalling, transcription regulation, biosyntheses of primary and secondary metabolites and other biological processes. The diversity of these genes reflects the complexity of stress responses. The expressed genes provide a basis for further identification of candidate genes that affect bark stripping through variation in their expression levels while the uncharacterized genes that responded to simulated herbivory and methyl jasmonate provide a rich resource for future studies. Gene expression can be used by breeders to exploit phenotype variability among individuals within or between populations. It also remains to be tested whether variations in the transcript levels, particularly the differentially expressed components in reponse to the artificial stress treatments can be linked to the susceptibility classes identified in the field [46] and whether they are heritable.