- Research article
- Open Access
Transcriptomic analysis of wound xylem formation in Pinus canariensis
BMC Plant Biologyvolume 17, Article number: 234 (2017)
Woody plants, especially trees, usually must face several injuries caused by different agents during their lives. Healing of injuries in stem and branches, affecting the vascular cambium and xylem can take several years. In conifers, healing takes place mainly from the remaining vascular cambium in the margin of the wound. The woundwood formed in conifers during healing usually presents malformed and disordered tracheids as well as abundant traumatic resin ducts. These characteristics affect its functionality as water conductor and its technological properties.
In this work we analyze for the first time the transcriptomic basis of the formation of traumatic wood in conifers, and reveal some differences with normal early- and late-wood. Microarray analysis of the differentiating traumatic wood, confirmed by quantitative RT-PCR, has revealed alterations in the transcription profile of up to 1408 genes during the first period of healing. We have grouped these genes in twelve clusters, according to their transcription profiles, and have distinguished accordingly two main phases during this first healing.
Wounding induces a complete rearrangement of the transcriptional program in the cambial zone close to the injuries. At the first instance, radial growth is stopped, and a complete set of defensive genes, mostly related to biotic stress, are induced. Later on, cambial activity is restored in the lateral borders of the wound, even at a high rate. During this second stage certain genes related to early-wood formation, including genes involved in cell wall formation and transcription factors, are significantly overexpressed, while certain late-wood related genes are repressed. Additionally, significant alterations in the transcription profile of abundant non annotated genes are reported.
Organisms usually suffer injuries throughout their life. In multicellular organisms these injuries can cause the damage or loss of differentiated tissues or organs, and ease the entry and spread of pathogens. The analysis of the similarities and differences in the regeneration process in animals and plants has been on the spotlight in recent years [1,2,3]. Animals can often regenerate these damaged tissues and even, in certain cases, the lost organs, and due to constant regeneration of certain tissues such as skin, eventually no signal of the injury remains after some time.
On their side, plants do not regenerate continuously their tissues. Proliferation is usually limited to certain niches: the root and shoot apical meristems (including axillary buds) and the vascular cambium and the phellogen in woody plants. If damage occurs, the plant generates new tissues and organs from these meristems or eventually develops new meristematic niches from living cells, usually parenchymatic ones [4, 5]. This is the case of the traumatic periderms developed, for instance, from cortical parenchyma or from mesophyll to seal a wound in a young stem or a leave.
When a woody branch or stem suffers a deep wound, affecting the secondary xylem, the vascular cambium must be restored. In certain angiosperms proliferation from xylem parenchyma or from immature xylem conducting elements has been described, as in Tilia , Eucommia  or Populus . These cells can reverse their differentiation pathway and divide profusely, giving rise to a parenchymatic callus. Later on, a new, traumatic vascular cambium differentiates within this callus, and new secondary xylem and phloem are produced.
On the contrary, this proliferation from (partially) differentiated cells is not usual in conifers. In these species healing takes place mainly from the remaining vascular cambium in the margins of the wound, as described recently for Pinus canariensis .
Anyway, the traumatic wood formed this way can be easily distinguished from normal wood. Traumatic wood usually presents malformed tracheary elements and fibers, with altered lignification patterns, and with a high proportion of parenchymatic cells. Orientation of these elements is also very often distorted [9,10,11,12], probably due to altered hormonal flux, and also to altered mechanical signals, as suggested by Chano et al. . This disorganized xylem implies an evident disadvantage for water and nutrient transport . In the case of conifers, especially Pinaceae, traumatic wood also presents a very high proportion of resin ducts, as described for Cedrus libani , Larix decidua [14, 15], Picea abies [15, 16], Pinus nigra  or Pinus pinaster . Actually, formation of traumatic resin ducts is the basis of traditional resin exploitation, very common in the past for several species of Mediterranean pines, and with increasing interest in the last years .
Since plants do not renew their secondary xylem, but generate new sheets of xylem, centrifugally, year after year, traumatic xylem also remains in the damaged branch or stem, leaving a “scar” in the wood. These scars have proven to be very useful, for instance, for dendrochronology studies [15, 19]. However, traumatic wood presents undesirable characteristics from a technological point of view. Although the higher density due to the increase in resin content can improve certain mechanical qualities of wood it also causes problems at machining and blunting . In addition, disordered and not properly formed traumatic tracheids contribute to alter the physico-mechanical properties of wood. Therefore, lumber dealers consider traumatic wood as a defect, lowering the price and reducing the applicability of wood pieces with important scars.
Several works have focused in the consequences of traumatisms on wood development in conifers: early-late wood ratio, ring width, formation of traumatic resin ducts (f.i., [10, 21,22,23,24]), and a few others in the description of the healing process from an anatomical point of view [9, 25,26,27]. However, although the molecular aspects of the response to traumatism has been analysed in different angiosperms (f. i., [4, 28,29,30,31]), the process has been less studied in gymnosperms, where most works have focused in the induction of traumatic resin ducts by traumatism, insect attack or fungal infections [32,33,34,35]. In this work we focus on the molecular basis of traumatic wood formation in a gymnosperm, P. canariensis, known for its extraordinary healing ability. For this purpose we have performed deep wounds in the stem of P. canariensis trees, affecting the vascular cambium, and have assessed the transcriptomic profile during the healing process and traumatic wood growth.
Results and discussion
Identification of genes induced and repressed in response to wounding
In order to analyse the transcriptomic response in the cambial zone and differentiating xylem in the borders of deep wounds performed in the stem of pine trees, we hybridized a 60K two-color cDNA microarray (Agilent, USA) which includes genes involved in P. canariensis xylogenesis . Samples were collected at three dates during wound response: i) H1 was collected seven days after wounding, ii) H2 after 75 days, when development of traumatic wood was evident and while the trees outside the wound area were still forming early wood, and finally iii) H3 92 days after wounding, when the trees were already forming late wood. Controls for each sample were collected at the same sampling dates from branches distant from the wound, in order to distinguish local effects caused by wound response from constitutive changes in gene expression during the vegetative season.
Figure 1 shows the distribution of genes selected as over- and underexpressed at each sampling point. We identified 1408 differentially expressed genes (DEG), genes significantly overexpressed or repressed compared to normal wood formation. Table 1 shows a selection of 91 DEGs with the strongest response (induction or repression), grouped following the functional processes they are presumably related to, according to their top BLASTx hit, as previously described . The complete table can be found in a supplementary table (Additional file 1).
Immediate response H1 included 837 DEGs; 619 DEGs were detected for H3, while only 336 DEGs were detected at H2. Just 69 genes were identified as DEG for the three sampling dates, H1, H2 and H3. Moreover, just 87 genes were identified as DEGs exclusively for H2, while 348 were exclusive DEGs for H3 and up to 658 for H1 (Fig. 2).
Enrichment analysis of DEGs pointed out an increase of mRNA levels for the categories “defense response”, “response to stress” and different forms of “response to stimulus”, as “response to abiotic stimulus” or “response to biotic stimulus”, among others, in the Biological Process (BP) category. As well, other enriched GO terms were “nucleic acid binding transcription factor activity” and “sequence-specific DNA binding transcription factor activity”, for the Molecular Function (MF) category, and “extracellular region”, “cell wall”, “cell periphery” and “external encapsulating structure” GO terms of the Cellular Component category.
Hierarchical clustering of DEGs
Twelve clusters were established according to the transcription patterns detected for DEGs throughout H1, H2 and H3 (Fig. 3).
Cluster A includes genes clearly induced at H1, which keep high transcription levels during H2 and H3, although at a minor degree. Genes included in clusters B, D and F were also overexpressed at H1, but later on their transcription levels decrease, being even repressed at H2 and/or H3. Genes in cluster C show a faint overexpression throughout the three phases, while clusters E and H show an increasing overexpression at H2 and H3. The opposite pattern is reflected by clusters G and K, with an increasing repression at H2 and H3. Genes in cluster I show a significant repression at H1, followed by a recovery of transcription to normal levels at H2 and H3. Finally, clusters J and L are characterized by a general repression during the three phases.
Clustering of samples revealed consistency among biological replicates, as shown in a supplementary figure (Additional file 2). Samples harvested at H1 clustered together and separated from the other sampling dates; on their side, H2 and H3 samples were included in the same group. Slight irregularities (for instance, sample Pc3H2 is closer to H3 samples than to the other H2 ones) can be due to the genetic variability among trees. This result supports the differentiation of two major phases (H1 and H2/H3) in the response to wounding, as discussed below.
To validate the reliability of the transcription profiles obtained from microarray hybridization, we selected 12 genes for qRT-PCR analysis, covering the main tendencies described above and the putative function of the genes. Thus, we selected three genes directly involved in cell growth and cell wall formation, as an expansin (Contig 03225, cluster H), a CeSA-like (Contig 00654, cluster L) and a CCoAOMT (Contig 06476, cluster G), transcription factors also involved in xylogenesis, as a MYB46-like (Contig 12050, cluster L), a WOX4-like (Contig 06813, cluster H), a bHLH35-like (Contig 05923, cluster C)an ATHB13-like (Contig 20304, cluster B), a NAC2-like (Contig 00787, cluster B), and a WRKY51-like (Contig 05551, cluster L). Finally, we also analysed a gene coding for a PAL protein (Contig 20555, cluster B), involved in salicylic acid biosynthesis and presumably related to defense, an EXORDIUM-like protein (Contig 09007, cluster G), presumably involved in cell proliferation, and a Major Allergen Pru AR1-like (Contig 22185, cluster A), putatively involved in defensive response.
Profiles obtained by qRT-PCR for these genes match the ones obtained from microarray hybridization, with high correlation coefficients, thus validating the general tendencies described above for microarray analysis (Fig. 4).
Detailed analysis of the transcription patterns of the differentially expressed genes (DEGs) leads to the identification of two major phases in the response to wounding.
H1. Immediate response
A complete rearrangement of the transcriptional program takes place as immediate response to wounding. At H1, a general repression of genes involved in the normal development of early wood is detected. In particular, genes related to meristematic activity, cell division or synthesis of cell wall show their transcription levels significantly lowered. This is consistent with anatomical observations: Chano et al.  described how cambial activity is stopped in response to a recent wound, and no growth in the cambial zone is further detected up to approximately 4 weeks after wounding. On the contrary, numerous genes putatively involved in the defense against stress (including biotic stress) are significantly induced at this point, serving as a defense against opportunistic pathogens infecting the wound. Noteworthy, many of these genes were reported to show their transcript maximum in normal xylogenesis during late wood formation . Expression of genes related to stress and defense processes in differentiating late wood has also been reported in other species. For instance, Mishima et al.  described the abundance of “defense mechanism genes” in the “cessation of growth clusters” obtained from cambial zone and differentiating xylem in Cryptomeria japonica In normal late wood formation, these genes could act as a preventive defense against putative pathogens infecting the tree just before dormancy, since dormancy could hamper the display of an induced response during winter. Later on, well differentiated late wood constitutes a barrier against eventual infections that could take place during winter, as described by the CODIT (Compartmentalization Of Decay In Trees) model .
Thus, among the genes involved in cell wall development typically overexpressed during early wood development and repressed at H1 we can find transcription factors such as the HD-ZIP class III family member ATHB15-like or a MYB46-like transcription factor, reported to be involved in cell wall biosynthesis in Arabidopsis , included in clusters I and L, respectively (Fig. 3). Other genes directly involved in the cell wall biosynthesis and repressed at this stage, were some CAZymes (f.i., Contigs 03231, 11436, 13611, 19457, 00766, or 08356), COBRA or KORRIGAN endoglucanase (Contigs 01405, or 10173, as well as a CCoAOMT (Contig 02447), involved in lignin biosynthesis and deposition. As well, a homologous of the rice NAC29 transcription factor (Contig 13895) has been found to be locally induced at H1 (FC=22.81). This gene was previously reported to participate in normal late-wood development in the Canary Island pine , coexpressed with a putative cellulose synthase-like protein, consistently with its reported role as activator of CesA in rice . However, no significant induction of CesA was found in H1. This fact can be due to the observed repression at this point of MYB546, another activator of CesA . Additionally, NAC29 could be also involved in other routes related to wound stress, and not only in the synthesis of cellulose, since no growth was detected at H1. Actually, this gene has been described to be involved in the response to stress caused by high salinity and drought in bread wheat .
Conversely, among the genes significantly overexpressed at H1 we can find genes related to oxidative stress, hydrolytic enzymes and hormonal signaling. Oxidative stress is one of the main effects of mechanical damage and infections. Cell lysis results in the production of hydrogen peroxide, which is toxic for pathogens, but also for plant cells, triggering the hypersensitive response . Peroxidases are then induced for ROS (reactive oxygen species) detoxification [43, 44]. Interestingly, in addition to their role in response to pathogens [45, 46], peroxidases are also involved in lignin biosynthesis and suberization [47,48,49]. Although some peroxidases were repressed at H1, several contigs coding for different isoforms of a peroxidase12-like protein were overexpressed at this time (f.i., Ppnisotig 01747 with FC value of 97.76, cluster D). In the same way, other genes involved in oxidative stress were induced at this stage, as Contig 03079 (cluster B), putatively coding for a lacoylglutathione lyase, involved in the glutathione-based detoxification  and in the response to drought and cold stress . In the same way, a glutathione-S-transferase (Ppnisotig 06171 found in cluster D) or a thioredoxin (Contig 09180, cluster B), also involved in the anti-oxidative plant defense , were also overexpressed at H1.
Another important group of genes induced at this stage are those coding for hydrolytic enzymes that attack pathogen cell wall. Among them we find Contig 18804 and Ppnisotig 13431 in cluster A, homologous to PI206, a disease resistance response protein firstly described in Pisum sativum, where is induced after inoculation with Fusarium solani [53, 54]. In the same way, and also in cluster A, the putative PR-4-like proteins Contig 22375, Ppnisotig 13133 and Contig 19053 are highly induced at H1 (FC values of 343.24, 206.36 and 96.55, respectively). PR-4 protein was first described in Solanum tuberosum , named also win-1 and win-2 for “wound-inducible genes”. In Capsicum chinense L., PR-4 was found to have both RNAse and DNAse activity in the extracellular space during stress conditions . Other putative PR4-like proteins with endochitinase activity , as Ppnisotig 08058 and Ppnisotig 00751, were also induced at H1 and found in cluster D. This cluster also includes other chitinases, such as endochitinase a-like proteins (Contig 21216, Contig 10307 and Contig 00126), and chitinase 1-like protein (Contig 23442). A major allergen pru ar1 homologous (Contig 22185) was also found in cluster A, strongly induced at H1 (FC 138.29). This protein was first described in Prunus armeniaca during ripening and annotated as a pathogenesis-related protein . Ppnisotig 12265, found in cluster G, corresponds to a putative antimicrobial peptide 1, which are widely present in living organisms, and possess antifungal and antibacterial properties . Finally, we can also mention Contig 00602 (cluster B) and Contig 17617 (cluster D), coding for two defensins, the most abundant antimicrobial peptides in plants, involved in defense-related processes, biotic stress response and plant development , which were also reported to be expressed during normal late wood differentiation in P. canariensis .
In the first stage after wounding the plant displays an extensive hormonal signaling. For instance, Contig 00715 and Ppnisotig 12073, included in clusters B and C, respectively, encode for putative 1-aminocyclopropane-1-carboxilate oxidase (ACO) proteins, involved in the synthesis of ethylene, known to be involved in different stress- and defense-related processes [61, 62].
Jasmonic acid (JA) is known to trigger a complex signaling network, both locally, activating the expression of wound-induced genes, and systemically, via the systemin peptide , mediated by ethylene . However, we have not detected any DEG related to JA biosynthesis. The restrictive criteria used in this work to select DEGs can account for this result. Additionally, a local repression of the JA-dependent pathway by ethylene production has been reported in Arabidopsis , where the existence of an additional JA-independent pathway has also been described. This could also be the case for Pinus canariensis.
Two genes coding for salicylic acid-binding protein 2-like (SABP2-like) proteins were overexpressed at H1 (Contig 03482, cluster B) and at H1 and H2 (Contig 14053, cluster D). These proteins are involved in the plant immune response, through their salicylic acid (SA)-stimulated lipase activity . SA is also involved in the expression of plant pathogenesis-related genes , and is thought to be an antagonistic of JA , blocking its synthesis [69, 70]. This would also be consistent with the lack of detection of DEGs related to the JA-dependent wound response pathway in this work.
Several non-annotated genes differentially overexpressed at H1 showed high levels of overexpression, specifically 63 DEGs in cluster B, 11 in cluster D, as many of the DEGs related to defense and stress mentioned above, and 2 in cluster F (f.i., Contigs 09209, 22448 or 24621 in cluster B, Contigs 23569, 22397 or 19474 in cluster D, with FC values over 20, or Contigs 03012 and 03111 in F). In addition, other non-annotated DEGs were repressed at H1, mainly grouped in clusters I (86 DEGs, with FC values below -10 for Contigs 02729, 13781, 19504 and 16419), and L (30 DEGs, with FC values close to -10 for Contigs 10360, 12514 or 14477). Additionally, other poorly annotated genes can be found among the H1-related DEGs. For instance, Contig 03506 was strongly induced at H1, with a FC value of 63.88, and kept overexpressed in H2 and H3. This sequence was annotated as homologous of the hypothetical protein SELMODRAFT_115352, predicted in the clubmoss Selaginella moellendorfii . As well, other remarkable contigs poorly annotated were Contig 20761 (cluster D), a predicted protein with FC value of 24.06 in H1, and Contig 34794 (cluster L), repressed to -10.94 at H1 and annotated as homologous to chickpea uncharacterized protein LOC101509257 .
H2/H3: Development of traumatic wood.
Chano et al.  described that noticeable formation of traumatic xylem begins 4 weeks after wounding. Accordingly, we collected traumatic wood samples 11 weeks after wounding (H2), when traumatic growth was visible at healing borders. At this date, early wood was still being formed . Two weeks later, when the trees were already differentiating late wood , another samples of traumatic wood were collected in independent wounds (H3).
As expected, after the first phase, characterised by the cessation of growth and by the expression of defensive genes, cambial activity resumes at the wound margins and development of traumatic wood is evident. Consistently, genes related to cell proliferation and cell wall biosynthesis are expressed. Thus, transcription patterns are more similar between H2 and H3, and differ more from H1.
While most of the genes involved in xylogenesis during latewood formation do not change their normal transcription patterns, and therefore are not detected as DEG at H2 or H3, several genes characteristic of early wood formation appear as overexpressed at these phases. This is the case for Contig 06813 (cluster H), encoding for a WOX4-like transcription factor. WOX4 belongs to the WUSCHEL-related HOMEOBOX (WOX) family, which is involved in the differentiation in the organizing center of the apical shoot , in procambial and cambial growth with function in vascular bundles development [74, 75] and in the regulation of proliferation from stem cells niches in root and shoot meristems after embryogenesis  together with CLAVATA (CLV; ). Moreover, a homologous of the clavata3-like protein (CLV3) was found to be induced at H3 as well (Contig 14178, cluster H), suggesting similar combined roles in response to wounding and meristematic activity during tissue regeneration and traumatic wood development. In the same way, homologues of two expansins (Contigs 03225 and 18811), two KORRIGAN endoglucanases (f.i., Contigs 18777 and 10173) or several CAZymes (f.i., Contigs 00603, 13281, 09907, 05424, or 05066), typically expressed during early wood formation in P. canariensis , are overexpressed at H3, when late wood is already differentiating in other parts of the stem. On the contrary, other CAZymes (Contigs 01916, 17013 or 21865) or a cellulose synthase (Contig 15857), typically expressed during late wood formation in P. canariensis , are repressed at H2 and H3, in the same way as an early wood induced CCoAOMT (Contig 06476), crucial in lignin biosynthesis. Repression and overexpression of putative early and late wood genes during H1, H2 and H3 are summarized in Fig. 5.
These results are consistent with anatomical observations. As shown in Fig. 6 not a clear difference between early and late wood is observed in the traumatic wood grown during 18 months after wounding. On the contrary, a high number of resin ducts appear in this traumatic wood, as already reported by Chano et al . Accordingly, several genes related to resin synthesis have been detected as overexpressed at H2 and H3. Oleoresins are one of the main conifer defenses against pathogens, avoiding the spread of infections. In this work we have detected DEGs encoding geranyl diphosphate synthase and geranylgeranyl diphosphate synthase (Contig 03270 and Ppnisotig 10634, respectively; cluster H), involved in the synthesis of mono and diterpenes, induced at H3. In the same way, in the same cluster, overexpressed at H3, appears Contig 08417, encoding an abietadienol/abietadienal oxidase–like protein, which catalyzes several oxidative steps in diterpenol biosynthesis . Contig 13499, encoding a (-)-camphene tricyclene synthase-like is also induced at H3, appearing in cluster C. This monoterpene synthase is involved in the synthesis of different monoterpenes, as camphene, tricyclene, limonene or myrcene .
Also induced at H2/H3 appear several genes presumably involved in ethylene synthesis, although this hormone is supposed to act in the first steps of the response . This is the case of ACS (1-aminocyclopropane-1-carboxylic acid synthase) or ACO (1-aminocyclopropane-1-carboxilate oxidase). In spruce and Douglas fir, for instance, multiple ACS genes and a single ACO gene were found to be induced during the immediate response to wounding . Conversely, we found multiple ACO genes overexpressed during the whole response (from H1 to H3), such as Contig 00524 (cluster C) and Contig 16100 (cluster I). This result suggests a different response in P. canariensis compared to Picea and Pseudotsuga, which could be related with the efficient healing response of Canary Island pine.
On the contrary, transcript levels of many other defensive genes overexpressed at H1 decrease to normal levels at H2/H3, or are even repressed at H3, when latewood is forming in other parts of the tree. During latewood formation defensive genes are expressed, as previously reported in P. canariensis  or C. japonica . This constitutive upregulation of these genes could account for the comparatively lower expression levels detected for traumatic wood formation. This is the case of two HSPRO genes (Contigs 02906 and 19857), induced at H1 and related to nematode resistance, or a noduline (Contig 12353), presumably involved in plant-microbe interactions, which show lower transcription levels at H2/H3 in the healing borders than in controls.
As exposed previously for the immediate response, other non-annotated DEGs were significantly overexpressed at H2 and/or H3. Thus, 26 non-annotated DEGs were included in cluster C, and 80 more in cluster H. For instance, Contig 21346, with FC values over 7 in H2 and H3, or Contig 12627, with a FC value of 4.42 in H2 and 8.62 in H3). Important numbers of non-annotated sequences were found in clusters G and L, where 58 and 70 DEGs, respectively, showed underexpression for H2 and/or H3 (f.i., Contigs 20076, 23934 or 24690 in cluster G were strongly repressed at H3, with FC values close to -10, or Contigs 14134 and 20478, repressed for H2 and H3 with FC values from -9.5 to -3.12 and from -5.43 to -2.76, respectively). As well, some poorly annotated contigs were remarkably repressed. For instance, Contig 02798 (cluster L) and Contig 12685 (cluster G), homologous to uncharacterized LOC101210414, and LOC101213469 from Cucumis sativus , showed FC values of -5.06 and -6.43 at H3, respectively.
Wounding induces a complete rearrangement of the transcriptional programme in the cambial zone close to the injuries. In particular, a considerable percentage of genes presumably involved in xylogenesis show an altered transcription pattern in response to wound and during healing.
At the first instance, radial growth is stopped in the vicinity of the wound, and a complete set of defensive genes, mostly related to biotic stress, are induced, as a barrier against opportunistic pathogens. Interestingly, some of these genes have also been reported to be preferentially transcribed in differentiating late wood. Later on, cambial activity is restored in the lateral borders of the wound, even at a higher rate than in other parts of the stem. This fast growth, which is dependent on the general health and reserves of the tree, eventually leads to the complete healing of the wound and restoration of the cambial ring. Anatomically, we have not detected a well-defined contrast among early and late wood in the traumatic wood formed during 18 months after wounding. During this period, most of the genes preferentially expressed during normal late wood development do not change their expression pattern, described in Chano et al. . However, a subset of genes shows their transcription levels significantly altered by wound and healing. Among them, it is noteworthy the presence of genes involved in cell wall formation. Thus, genes coding for CAZymes and cellulose synthases overexpressed in normal late wood formation are comparatively repressed in traumatic wood. Conversely, similar genes typical of early wood keep their high transcription levels in traumatic wood, even at the moment of late wood formation. On the contrary, an early wood CCoAOMT, involved in lignin biosynthesis, is also repressed in traumatic wood. These genes, together with many others non-annotated yet, but showing similarly modified transcription patterns in healing tissue, probably underlie the anomalous characteristics of traumatic wood. In the same way, we cannot discard that other genes not detected as DEG due to the restrictive criteria used in this work could still play a biologically significant role in the wound wood formation process.
Our results suggest that the tree, after the synthesis of defensive molecules against eventual pathogens, and once cambial activity is restored at the wound borders, produces a fast growing traumatic wood. This tissue, in which annual rings are not clearly distinguished, at least the first year, could be less efficient as preventive barrier than normal late wood regarding secondary wall lignification, but it presents a high proportion of resin ducts, and also provides a good way to heal the wound in the shortest possible time. Further investigations are needed to clarify this point.
Plant material and wounding
For this work we used 3 Pinus canariensis trees, 5 years old. Pines were grown in greenhouse, using 650 ml conical containers with 3:1 (v/v) peat:vermiculite. After the first year, trees were transferred to soil in experimental garden at UPM facilities, and grown under environmental conditions. At the moment of the beginning of this experiment, trees were approximately 2 m high and 7-10 cm diameter at the base. Using a sterile scalpel, we performed two wounds, removing bark, phloem, vascular cambium and first rows of xylem from a rectangular window 10 cm high and spanning half the circumference of the stems (Fig. 7). Wounds were performed in opposite sides of the stem and with an interval of approximately three wound heights.
Samples were collected according to the described seasonal growth and healing patterns described previously for the species [9, 36] Wounding was performed on April 9th, when cambial activity was ongoing. One week later we collected a frame of tissue from the wound margins in both wounds (H1); at this moment formation of a first traumatic tissue can be expected . On June 25th, when the callous tissue was emerging from the margins of wounds and the trees were at the end of the early-wood development period, we collected the tissues growing in the margins of one wound per tree (H2). Later on, on July 9th, concurring with the late-wood development period , we collected the thick callous tissues growing in the frame of remaining wounds (H3).
For each sample, controls were collected at the same sampling dates, from branches away from the wounds, in order to distinguish transcriptomic changes of the vegetative growth of those caused by wounding. The tissue collected for control samples included bark, phloem, vascular cambium and the most external layers of xylem. Collected samples were processed individually, immediately frozen in liquid nitrogen and stored a -80°C.
RNA isolation and sequencing
Total RNA was isolated from each sample, using the CTAB-LiCl precipitation method , and purified with the RNeasy Plant Mini Kit (Qiagen, CA, USA). Quantity of total RNA for each sample was measured with Nanodrop model ND-1000 (Thermo Scientific, MA, USA), and RNA quality was checked using Experion Bioanalyzer (Bio-Rad, CA, USA).
A set of 15266 contigs involved in meristematic activity of Pinus canariensis, selected from a previous work , was used for the design of a two-color 60K microarray (Agilent, USA). Furthermore, we added 2303 contigs from other cDNA libraries of P. pinea, as well ESTs and sequences of the loblolly pine from the Pine Gene Index Database (http://www.mgel.msstate.edu/dna_libs.htm). For each contig, one 60 bp long probe was designed and spotted at least 3 times on the slide. Probes designed for Populus, mouse and human ESTs available in public databases were included as negative controls.
For each sampling point (H1, H2, H3), the three biological replicates were hybridized (wound vs. control) following the two-color protocol provided by the manufacturer (Agilent Technologies, CA, USA), and images were captured with a GenePix 4000B (Axon, CA, USA), and spots were quantified using the GenePix software (Axon, CA, USA). Microarray data was uploaded to the NCBI’ Gene Expression Omnibus and are accessible through the GEO series accession number GSE102275 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE102275).
Background correction and normalization of expression data were performed using LIMMA (Linear Models for Microarray Data; ). For local background correction and normalization, the methods “normexp” and “loess” in LIMMA were used, respectively. To achieve similar distribution across arrays and consistency among arrays, log-ratio values were scaled using the median-absolute value as scale estimator.
The non-parametric algorithm “Rank Products”, available as a package for Bioconductor in R [84,85,86], was used for evaluation of Differentially Expressed Genes (DEGs). This method detects genes that are consistently high ranked in a number of replicated experiments independently of their numerical intensities. Results are provided in the form of P adjusted by False Discovery Rate (FDR), defined as the probability of a given gene is ranked in the observed position by chance.
Those probes with a FC above 2 and below -2, with a significance level FDR below 0.05, were selected as differentially expressed. Thus, technical replicates were merged into one value per contig, and a datamatrix formed by ratios between experimental and control measurements for selected Differentially Expressed Genes (DEGs), including time sampled and biological replicate, was created. Clustering was performed in R, and the heatmap was plotted using the heatmap.2 function of the gplots package . Enrichment analysis of DEGs was performed using Blast2GO v.2.7.2 as well.
The expression patterns of 12 DEGs covering the main profiles obtained from microarrays were confirmed by qRT-PCR using the same RNA employed for microarray hybridizations. First strand cDNA synthesis was performed using SuperScript™ III reverse transcriptase (Invitrogen, USA) following manufacturer’s instructions and using 4 μgr of total RNA and random hexamers. Gene specific primers were designed for twelve selected DEGs (Table 2) using the Primer3 software , with a melting temperature between 60 and 65° C, and producing amplicons between 80 and 120 bp. qRT-PCR was performed in a CFX96™ Real-Time PCR Detection System (Biorad, USA), using the SsoFast™ EVAgreen® Supermix (Biorad, USA), according to manufacturer’s protocol, and following the standard thermal profile: 95° C for 3 min, 40 cycles of 95° C for 10 s and 60° C for 10 s. In order to compare data from different qRT-PCR runs, the CT values were normalized using the Ri18S as housekeeping gene, whose specific primers were FW 5’-GCGAAAGCATTTGCCAAGG-3’ and REV 5’-ATTCCTGGTCGGCATCGTTTA-3’. This genes has been previously proved to be useful for this purpose in pine species (f.i., see Perdiguero et al. ). The expression ratios were then obtained using the delta-delta-CT method corrected for the PCR efficiency for each DEG .
Birnbaum KD, Alvarado AS. Slicing across Kingdoms: Regeneration in Plants and Animals. Cell. 2008:697–710.
Sugimoto K, Gordon SP, Meyerowitz EM. Regeneration in plants and animals: Dedifferentiation, transdifferentiation, or just differentiation? Trends in Cell Biology. 2011:212–8.
Sánchez Alvarado A, Yamanaka S. Rethinking differentiation: Stem cells, regeneration, and plasticity. Cell. 2014:110–9.
Sena G, Wang X, Liu H-Y, Hofhuis H, Birnbaum KD. Organ regeneration does not require a functional stem cell niche in plants. Nature. 2009;457:1150–3.
Sena G, Birnbaum KD. Built to rebuild: In search of organizing principles in plant regeneration. Current Opinion in Genetics and Development. 2010:460–5.
Stobbe H. Developmental Stages and Fine Structure of Surface Callus Formed after Debarking of Living Lime Trees (Tilia sp.). Ann Bot. 2002;89:773–82.
Pang Y, Zhang J, Cao J, Yin SY, He XQ, Cui KM. Phloem transdifferentiation from immature xylem cells during bark regeneration after girdling in Eucommia ulmoides Oliv. J Exp Bot. 2008;59:1341–51.
Zhang J, Gao G, Chen JJ, Taylor G, Cui KM, He XQ. Molecular features of secondary vascular tissue regeneration after bark girdling in Populus. New Phytol. 2011;192:869–84.
Chano V, López R, Pita P, Collada C, Soto Á. Proliferation of axial parenchymatic xylem cells is a key step in wound closure of girdled stems in Pinus canariensis. BMC Plant Biol. 2015;15:64.
Arbellay E, Stoffel M, Sutherland EK, Smith KT, Falk DA. Changes in tracheid and ray traits in fire scars of North American conifers and their ecophysiological implications. Ann Bot. 2014;114:223–32.
Zajaczkowska U. Regeneration of Scots pine stem after wounding. IAWA J. 2014;35:270–80.
Zajaczkowska U. Overgrowth of Doublas fir (Pseudotsuga menziensii Franco) stumps with regenerative tissue as an example of cell ordering and tissue reorganization. Planta. 2014;240:1203–11.
Fahn A, Werker E, Ben-Tzur P. Seasonal effects of wounding and growth substances on development of traumatic resin ducts in Cedrus libani. New Phytol. 1979;82:537–44.
Bollschweiler M, Stoffel M, Schnewly D, Bourqui K. Traumatic resin ducts in Larix decidua stems impacted by debris flows. Tree Physiol. 2008;28:255–63.
Stoffel M. Dating past geomorphic processes with tangential rows of traumatic resin ducts. Dendrochronologia. 2008;26:53–60.
Nagy NE, Franceschi VR, Solheim H, Krekling T, Christiansen E. Wound-induced traumatic resin duct development in stems of Norway spruce (Pinaceae): Anatomy and cytochemical traits. Am J Bot. 2000;87:302–13.
Luchi N, Ma R, Capretti P, Bonello P. Systemic induction of traumatic resin ducts and resin flow in Austrian pine by wounding and inoculation with Sphaeropsis sapinea and Diplodia scrobiculata. Planta. 2005;221:75–84.
Rodríguez-García A, López R, Martín JA, Pinillos F, Gil L. Resin yield in Pinus pinaster is related to tree dendrometry, stand density and tapping-induced systemic changes in xylem anatomy. For Ecol Manage. 2014;313:47–54.
Stoffel M, Klinkmüller M. 3D analysis of anatomical reactions in conifers after mechanical wounding: first qualitative insights from X-ray computed tomography. Trees. 2013;27:1805–11.
García-Iruela A, Esteban L, de Palacios P, García-Fernández F, de Miguel Torres Á, Vázquez-Iriarte E, Simón C. Resinous wood of Pinus pinaster Ait.: physico-mechanical properties. Bioresources. 2016;11:5230–41.
Gärtner H, Heinrich I. The formation of traumatic rows of resin ducts in Larix decidua and Picea abies (Pinaceae) as a result of wounding experiments in the dormant season. IAWA J. 2009;30:199–215.
Schneuwly DM, Stoffel M, Dorren LKA, Berger F. Three-dimensional analysis of the anatomical growth response of European conifers to mechanical disturbance. Tree Physiol. 2009;29:1247–57.
Schneuwly DM, Stoffel M, Bollschweiler M. Formation and spread of callus tissue and tangential rows of resin ducts in Larix decidua and Picea abies following rockfall impacts. Tree Physiol. 2009;29:281–9.
Ballesteros JA, Stoffel M, Bodoque JM, Bollschweiler M, Hitz O, Díez-Herrero A. Changes in Wood Anatomy in Tree Rings of Pinus pinaster Ait. Following Wounding by Flash Floods. Tree-Ring Res. 2010;66:93–103.
Mullick DB. A new tissue essential to necrophylactic periderm formation in the bark of fou conifers. Can J Bot. 1975;53:2443–57.
Oven P, Torelli N. Wound response of the bark in healthy and declining silver firs (Abies alba). IAWA J. 1994;15:407–15.
Oven P, Torelli N. Response of the cambial zone in conifers to wounding. Phyton-Ann. REI Bot. 1999;39:133–7.
Xu J, Hofhuis R, Sauer M, Friml J, Scheres B. A molecular framework for plant regeneration. Science. 2006;311:386–8.
Sugimoto K, Jiao Y, Meyerowitz EM. Arabidopsis regeneration from multiple tissues occurs via a root development pathway. Dev Cell. 2010;18:463–71.
Wan X, Landhäusser SM, Lieffers VJ, Zwiazek JJ. Signals controlling root suckering and adventitious shoot formation in aspen (Populus tremuloides). Tree Physiol. 2006;26:681–7.
Asahina M, Azuma K, Pitaksaringkarn W, Yamazaki T, Mitsuda N, Ohme-Takagi M, Yamaguchi S, Kamiya Y, Okada K, Nishimura T, Koshiba T, Yokota T, Kamada H, Satoh S. Spatially selective hormonal control of RAP2.6L and ANAC071 transcription factors involved in tissue reunion in Arabidopsis. Proc Natl Acad Sci U S A. 2011;108:16128–32.
Fäldt J, Martin D, Miller B, Rawat S, Bohlmann J. Traumatic resin defense in Norway spruce (Picea abies): Methyl jasmonate-induced terpene synthase gene expression, and cDNA cloning and functional characterization of (+)-3-carene synthase. Plant Mol Biol. 2003;51:119–33.
Krokene P, Nagy NE, Solheim H. Methyl jasmonate and oxalic acid treatment of Norway spruce: anatomically based defense responses and increased resistance against fungal infection. Tree Physiol. 2008;28:29–35.
Mckay SAB, Hunter WL, Godard K, Wang SX, Martin DM, Bohlmann J, Plant AL. Insect attack and wounding induce traumatic resin duct development and gene expression of (-)-pinene synthase in Sitka spruce. Plant Physiol. 2003;133:368–78.
Zulak KG, Bohlmann J. Terpenoid biosynthesis and specialized vascular cells of conifer defense. Journal of Integrative Plant Biology. 2010:86–97.
Chano V, López de Heredia U, Collada C, Soto Á: Transcriptomic analysis of juvenile wood formation during the growing season in Pinus canariensis. Holzforschung 2017, aop:1–19. doi: 10.1515/hf-2017-0014
Mishima K, Fujiwara T, Iki T, Kuroda K, Yamashita K, Tamura M, Fujisawa Y, Watanabe A. Transcriptome sequencing and profiling of expressed genes in cambial zone and differentiating xylem of Japanese cedar (Cryptomeria japonica). BMC Genomics. 2014;15:219.
Shigo AL. Compartmentalization: A Conceptual Framework for Understanding How Trees Grow and Defend Themselves. Annu Rev Phytopathol. 1984;22:189–214.
Zhong R, Richardson EA, Ye ZH. The MYB46 transcription factor is a direct target of SND1 and regulates secondary wall biosynthesis in Arabidopsis. Plant Cell. 2007;19:2776–92.
Huang D, Wang S, Zhang B, Shang-Guan K, Shi Y, Zhang D, Liu X, Wu K, Xu Z, Fu X, Zhou Y, Gibberellin-Mediated A. DELLA-NAC Signaling Cascade Regulates Cellulose Synthesis in Rice. Plant Cell. 2015;27:1681–96.
Xu Z, Gongbuzhaxi, Wang C, Xue F, Zhang H, Ji W: Wheat NAC transcription factor TaNAC29 is involved in response to salt stress. Plant Physiol Biochem 2015, 96:356–363.
Levine A, Tenhaken R, Dixon R, Lamb C. H2O2 from the oxidative burst orchestrates the plant hypersensitive disease resistance response. Cell. 1994;79:583–93.
Diehn SH, Burkhart W, Graham JS. Purification and partial amino acid sequence of a wound-inducible, developmentally regulated anionic peroxidase from soybean leaves. Biochem Biophys Res Commun. 1993;195:928–34.
Mohan R, Vijayan P, Kolattukudy PE. Developmental and tissue-specific expression of a tomato anionic peroxidase (TAP1) gene by a minimal promoter, with wound and pathogen induction by an additional 5’-flanking region. Plant Mol Biol. 1993;22:475–90.
Gunnar Fossdal C, Sharma P, Lönneborg A. Isolation of the first putative peroxidase cDNA from a conifer and the local and systemic accumulation of related proteins upon pathogen infection. Plant Mol Biol. 2001;47:423–35.
Wang JE, Liu KK, Li DW, Zhang YL, Zhao Q, He YM, Gong ZH. A novel peroxidase CanPOD gene of pepper is involved in defense responses to Phytophtora capsici infection as well as abiotic stress tolerance. Int J Mol Sci. 2013;14:3158–77.
Hiraga S, Sasaki K, Ito H, Ohashi Y, Matsui H. A Large Family of Class III Plant Peroxidases. Plant Cell Physiol. 2001;42:462–8.
Valério L, De Meyer M, Penel C, Dunand C. Expression analysis of the Arabidopsis peroxidase multigenic family. Phytochemistry. 2004;65:1331–42.
Passardi F, Cosio C, Penel C, Dunand C. Peroxidases have more functions than a Swiss army knife. Plant Cell Rep. 2005;24:255–65.
Thornalley PJ. Glutathione-dependent detoxification of α-oxoaldehydes by the glyoxalase system: involvement in disease mechanisms and antiproliferative activity of glyoxalase I inhibitors. Chem Biol Interact. 1998;111–112:137–51.
Seki M, Narusaka M, Abe H, Kasuga M, Yamaguchi-Shinozaki K, Carninci P, Hayashizaki Y, Shinozaki K. Monitoring the expression pattern of 1300 Arabidopsis genes under drought and cold stresses by using a full-length cDNA microarray. Plant Cell. 2001;13:61–72.
Meyer Y, Siala W, Bashandy T, Riondet C, Vignols F, Reichheld JP. Glutaredoxins and thioredoxins in plants. Biochim Biophys Acta. 2008;1783:589–600.
Riggelman R, Fristensky B, Hadwiger L. The disease resistance response in pea is associated with increased levels of specific mRNAs. Plant Mol Biol. 1985;4:81–6.
Culley DE, Horovitz D, Hadwiger LA. Molecular characterization of disease-resistance response gene DDR206-d from Pisum sativum (L.). Plant Physiol. 1995;107:301–2.
Stanford A, Bevan M, Northcote D. Differential expression within a family of novel wound-induced genes in potato. Mol Gen Genet. 1989;215:200–8.
Guevara-Morato MÁ, De Lacoba MG, García-Luque I, Serra MT. Characterization of a pathogenesis-related protein 4 (PR-4) induced in Capsicum chinense L3 plants with dual RNase and DNase activities. J Exp Bot. 2010;61:3259–71.
Seo PJ, Lee AK, Xiang F, Park CM. Molecular and functional profiling of Arabidopsis pathogenesis-related genes: Insights into their roles in salt response of seed germination. Plant Cell Physiol. 2008;49:334–44.
Mbéguié-A-Mbeguié D, Gomez RM, Fils-Lycaon B: Sequence of an allergen-, stress-, and pathogenesis-related protein from apricot fruit (Accession No. U93165). Gene expression during fruit ripening (PGR 97-180). Plant physiol 1997, 115:1730.
Castro MS, Fontes W. Plant defense and antimicrobial peptides. Protein Pept Lett. 2005;12:11–6.
Tam JP, Wang S, Wong KH, Tan WL. Antimicrobial peptides from plants. Pharmaceuticals. 2015:711–57.
Yuan S, Wang Y, Dean JFD. ACC oxidase genes expressed in the wood-forming tissues of loblolly pine (Pinus taeda L.) include a pair of nearly identical paralogs (NIPs). Gene. 2010;453:24–36.
Hudgins JW, Ralph SG, Franceschi VR, Bohlmann J. Ethylene in induced conifer defense: cDNA cloning, protein expression, and cellular and subcellular localization of 1-aminocyclopropane-1-carboxylate oxidase in resin duct and phenolic parenchyma cells. Planta. 2006;224:865–77.
Rojo E, León J, Sánchez-Serrano JJ. Cross-talk between wound signalling pathways determines local versus systemic gene expression in Arabidopsis thaliana. Plant J. 1999;20:135–42.
O’Donnell PJ, Calvert C, Atzorn R, Wasternack C, Leyser HMO, Bowles DJ. Ethylene as a Signal Mediating the Wound Response of Tomato Plants. Science. 1996;274:1914–7.
León J, Rojo E, Sánchez-Serrano JJ. Wound signalling in plants. J Exp Bot. 2001;52:1–9.
Kumar D, Klessig DF. High-affinity salicylic acid-binding protein 2 is required for plant innate immunity and has salicylic acid-stimulated lipase activity. Proc Natl Acad Sci. 2003;100:16101–6.
Ward E, Uknes S, Williams S, Dincher S, Wiederhold D, Alexander D, Ahl-Goy P, Metraux J, Ryals J. Coordinate Gene Activity in Response to Agents That Induce Systemic Acquired Resistance. Plant Cell. 1991;3:1085–94.
Vidhyasekaran P. Salicyclic acid signalling in plant innate immunity. In: Plant Hormone Signalling Systems in Plant Innate Immunity. Dordrecht: Springer; 2015. p. 27–122.
Pena-Cortés H, Albrecht T, Prat S, Weiler EW, Willmitzer L. Aspirin prevents wound-induced gene expression in tomato leaves by blocking jasmonic acid biosynthesis. Planta. 1993;191:123–8.
Spoel SH, Koornneef A, Claessens SMC, Korzelius JP, Van Pelt JA, Mueller MJ, Buchala AJ, Métraux J-P, Brown R, Kazan K, Van Loon LC, Dong X, Pieterse CMJ. NPR1 modulates cross-talk between salicylate- and jasmonate-dependent defense pathways through a novel function in the cytosol. Plant Cell. 2003;15:760–70.
Banks JA, Nishiyama T, Hasebe M, Bowman JL, Gribskov M, dePamphilis C, Albert VA, Aono N, Aoyama T, Ambrose BA, Ashton NW, Axtell MJ, Barker E, Barker MS, Bennetzen JL, Bonawitz ND, Chapple C, Cheng C, Correa LGG, Dacre M, DeBarry J, Dreyer I, Elias M, Engstrom EM, Estelle M, Feng L, Finet C, Floyd SK, Frommer WB, Fujita T, et al. The Selaginella genome identifies genetic changes associated with the evolution of vascular plants. Science. 2011;332:960–3.
Varshney RK, Song C, Saxena RK, Azam S, Yu S, Sharpe AG, Cannon S, Baek J, Rosen BD, Tar’an B, Millan T, Zhang X, Ramsay LD, Iwata A, Wang Y, Nelson W, Farmer AD, Gaur PM, Soderlund C, Penmetsa R V, Xu C, Bharti AK, He W, Winter P, Zhao S, Hane JK, Carrasquilla-Garcia N, Condie JA, Upadhyaya HD, Luo MC, et al.: Draft genome sequence of chickpea (Cicer arietinum) provides a resource for trait improvement. Nat Biotechnol 2013, 31:240–246.
Mayer K, Schoof H, Haecker A, Lenhard M, Jürgens G, Laux T. Role of WUSCHEL in regulating stem cell fate in the Arabidopsis shoot meristem. Cell. 1998;95:805–15.
Ji J, Strable J, Shimizu R, Koenig D, Sinha N, Scanlon MJ. WOX4 Promotes Procambial Development. Plant Physiol. 2010;152:1346–56.
Ji J, Shimizu R, Sinha N, Scanlon MJ. Analyses of WOX4 transgenics provide further evidence for the evolution of the WOX gene family during the regulation of diverse stem cell functions. Plant Signal Behav. 2010;5:916–20.
Haecker A, Gross-Hardt R, Geiges B, Sarkar A, Breuninger H, Herrmann M, Laux T. Expression dynamics of WOX genes mark cell fate decisions during early embryonic patterning in Arabidopsis thaliana. Development. 2004;131:657–68.
Miwa H, Kinoshita A, Fukuda H, Sawa S. Plant meristems: CLAVATA3/ESR-related signaling in the shoot apical meristem and the root apical meristem. J Plant Res. 2009;122:31–9.
Ro D-K, Arimura G-I, Lau SYW, Piers E, Bohlmann J. Loblolly pine abietadienol/abietadienal oxidase PtAO (CYP720B1) is a multifunctional, multisubstrate cytochrome P450 monooxygenase. Proc Natl Acad Sci U S A. 2005;102:8060–5.
Falara V, Akhtar TA, Nguyen TTH, Spyropoulou EA, Bleeker PM, Schauvinhold I, Matsuba Y, Bonini ME, Schilmiller AL, Last RL, Schuurink RC, Pichersky E. The Tomato Terpene Synthase Gene Family. Plant Physiol. 2011;157:770–89.
Ralph SG, Hudgins JW, Jancsik S, Franceschi VR, Bohlmann J. Aminocyclopropane carboxylic acid synthase is a regulated step in ethylene-dependent induced conifer defense. Full-length cDNA cloning of a multigene family, differential expression, and cellular and subcellular localization in spruce and Douglas fir. Plant Physiol. 2007;143:410–24.
Huang S, Li R, Zhang Z, Li L, Gu X, Fan W, Lucas WJ, Wang X, Xie B, Ni P, Ren Y, Zhu H, Li J, Lin K, Jin W, Fei Z, Li G, Staub J, Kilian A, van der Vossen E a G, Wu Y, Guo J, He J, Jia Z, Ren Y, Tian G, Lu Y, Ruan J, Qian W, Wang M, et al.: The genome of the cucumber, Cucumis sativus L. Nat Genet 2009, 41:1275–1281.
Chang S, Puryear J, Cairney J. A simple and efficient method for isolating RNA from pine trees. Plant Mol Biol Report. 1993;11:113–6.
Smyth GK: Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol 2004, 3:Article3.
Hong F, Breitling R, McEntee CW, Wittner BS, Nemhauser JL, Chory J. RankProd: A bioconductor package for detecting differentially expressed genes in meta-analysis. Bioinformatics. 2006;22:2825–7.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, Leisch F, Li C, Maechler M, Rossini AJ, Sawitzki G, Smith C, Smyth G, Tierney L, Yang JYH, Zhang J. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5:R80.
R Core Team: R: A Language and Environmental for Statistical Computing. R Foundation for Statistical Computing, Vienna (Austria) 2013.
Warnes GR, Bolker B, Bonebakker L, Gentleman R, Liaw WHA, Lumley T, Maechler M, Magnusson A, Moeller S, Schwartz M, Venables B: gplots: Various R Programming Tools for Plotting Data. R package version 3.0.1, 2016.
Untergasser A, Cutcutache I, Koressaar T, Ye J, Faircloth BC, Remm M, Rozen SG. Primer3-new capabilities and interfaces. Nucleic Acids Res. 2012;40:e115.
Perdiguero P, Barbero MDC, Cervera MT, Collada C, Soto A. Molecular response to water stress in two contrasting Mediterranean pines (Pinus pinaster and Pinus pinea). Plant Physiol Biochem. 2013;67:199–208.
Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29:e45.
Authors thank the editor and the two anonymous referees for their helpful comments and suggestions. This work has been funded through the projects AGL2009-10606 (Spanish Ministry of Science and Innovation) and SPIP2014-01093 (Spanish National Parks Agency, Ministry of Agriculture). VC had a pre-doctoral fellowship granted by the Spanish Ministry of Science and Innovation.
Availability of data and materials
The datasets generated during the current study are available in the NCBI’ Gene Expression Omnibus and are accessible through the GEO series accession number GSE102275 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE102275).
The authors declare that they have no competing interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.