Differential gene expression analysis of the resprouting process in Pinus canariensis provides new insights into a rare trait in conifers

Resprouting is crucial in population dynamics in response to wildfires or herbivory. This trait, common in angiosperms, is considered rare in conifers, being Pinus canariensis one of the few species able to resprout. We analyzed gene expression during wound-induced resprouting in 5 years-old Canarian pines. RNA was extracted at three different stages, including control samples from remote branches, representing immediate response to wounding (R0), resprouting initiation (R1), and elongation (R2), and then hybridized to a microarray designed with 15,000 cDNAs from Pinus canariensis transcriptome of meristematic activity. We found 1164 Differentially Expressed Genes (DEGs), of which 903 were significant in R0, 278 in R1, and 261 in R2. Genes related to defense- and stress-response were mainly found to be overexpressed at R0, while DEGs coding for transcription factors such as GROWTH-REGULATING FACTOR (GRF), FLOWERING-PROMOTING FACTOR (FPF), and the HOMEODOMAIN LEUCINE ZIPPER Class IV (HD-ZIP IV), mainly related to outgrowth processes and lateral organogenesis in plants, showed overexpression at R1 and R2, when new shoots were emerging. This data was compared to differential gene expression during apical growth in P. canariensis, suggesting similarities between vegetative apical growth and regulation of lateral shoot development in response to mechanical wounding, but also yielding 32 DEGs exclusively overexpressed at some point during resprouting stages (R1 and/or R2). The implication of such transcription factor families, with important roles in leaf and inflorescence development also in flowering plants, suggests underlying homologies between common lateral organogenesis processes in angiosperms and resprouting in gymnosperms.


Introduction
Resprouting is a common process in woody angiosperms, recognized as an important feature in plant ecology and population dynamics (Noble and Slatyer 1980;Loehle 2000;Hoffmann et al. 2009;Matula et al. 2020; see also review by Clarke et al. 2013). However, resprouting has been considered as a rare feature in gymnosperms (Pausas and Keely 2017) and not many conifers have been recognized as good resprouters, according to the current literature and despite evidence of vegetative outgrowth and resprouting in fossil records (see review by Burrows 2021). Resprouting is a tolerance trait especially important in response to traumatisms and mechanical stress caused by forest fires or herbivory, among other disturbance factors, and it is produced due to the presence of vegetative latent buds that are able to produce regrowth and de novo develop into a new epicormic shoot (Clarke et al. 2013). Depending on the position of resprouts, species can be classified as stem, apical, basal or root resprouters (Clarke et al. 2013). For the latter, belowground tissues are well protected and isolated by soil. However, intense fires such as crown-fires can damage aerial organs in above-ground tree resprouters, causing necrosis and stem hydraulic failures (Fernandes et al. 2008;Midgley et al. 2011). Many of these tree species present self-pruning trunks, preventing canopy fires and therefore avoiding damages to their shoot apical meristems (SAM), located in the top part of the crown. On its side, thick barks confer a better protection for lateral meristematic tissues and latent buds, allowing them to resprout (Schwilk and Ackerly 2001).
Stem resprouting is usually produced from preformed latent buds located in leaf axils that are embedded and protected by bark (Clarke et al. 2013), but for some species like Wollemi pine (Wollemia nobilis W.G.Jones, K.D.Hill & J.M.Allen) buds are not preformed and epicormic shoots arise from small meristems located in these axils (Burrows et al. 2003). Moreover, epicormic resprouting has been proposed to occur in response to the physiological disequilibrium caused by environmental factors, including mechanical damage. For instance, insect attack or intense crown fire may cause removal of the apical shoot or defoliation, causing therefore a reduced photosynthetic capacity and the necessity of expanding foliar coverage and recovering canopy (Meier et al. 2012). Latent buds, which were in dormancy by auxin produced in the apex and leaves, are then active, inducing outgrowth of epicormic shoots (Ishii et al. 2011). Additionally, growth by cambial activity may compete with resprouting for nutrients and energy, therefore the reduction or cessation of cambial growth has been suggested to induce resprouting (Bachelard 1969;Nicolini et al. 2001;Ishii et al. 2011). Subsequently, mechanical damage which does not imply a major threat for the plant but causing a reduction of cambial activity (Larson 1994) may also induce resprouting.
Molecular basis of lateral organ initiation has been addressed in model organisms. For instance, some studies in crop science have focused on the plant architecture of Brassica napus L., including genome-wide association studies (GWAS) using a single nucleotide polymorphism (SNP) array of candidate genes for the regulation of flower development (Li et al. 2016). Moreover, transcriptomics of the axillary bud development under sugar and hormone signaling in B. napus revealed important roles of key regulators and different transcription factor families Li et al. 2020). Similarly, gene expression of axillary bud outgrowth in response to exogenous application of strigolactone, auxin and auxin transport inhibitor naphtyl phthalamic acid was studied in the ornamental crop Chrysanthemum morifolium Ramat. (Dierck et al. 2018). As well, the role of auxin, and its interaction with signaling pathways and regulatory genes during phyllotaxis and architectural shoot development from SAM has been also widely studied in A. thaliana (Vernoux et al. 2011; see also review by Traas 2019). However, little is known about the regulatory mechanisms triggering resprouting in tree species in response to external stimuli, despite the great interest that this process and outgrowth of new branches may arouse in forestry, playing a key role in forest and woodland management through coppicing (Dihn et al. 2019), and showing a great potential in tree breeding programs to enhance production.
In gymnosperms, resprouting has been suggested to be scarce by many authors, as reviewed recently by Burrows (2021). Not without reason, resprouting has been considered an ancestral condition that may have played a key role in adaptation to recurrent fires in the Paleozoic Era (He et al. 2016), but probably lost in modern conifers (Bond and Midgley 2003). One of the few conifers which is considered as a good resprouter in both juvenile and adult stages is the Canary Island pine (Pinus canariensis C.Sm. ex DC; Pausas and Keely 2017). In this species, the resprouting capability has been proposed to be linked to its evolutionary history in a volcanic environment under selective pressures exerted by perturbation regimes and driving the dynamics of its populations (Climent et al. 2004). Moreover, small injuries promote the cessation of cambial activity in the wounded area (Chano et al. 2015), including the repression of genes involved in radial growth (Chano et al. 2017a, b). In this work, we aim at a better understanding of the gene expression dynamics of local resprouting response to mechanical damage in P. canariensis. For this purpose, we used a microarray designed using a normalized transcriptome of the meristematic activity (Chano, Lopez de Heredia et al. 2017). We studied gene expression profiles after wounding and during the very early stages of the resprouting process and compared these results with those obtained during apical growth, co-analyzing expression profiles of genes suggested to be involved in one or both processes as a first approximation to identify resprouting specific genes.

Plant material, wounding, and sample collection
Using a sterile scalpel, we produced mechanical injuries in three 5-year-old Canary Island pines growing in an experimental garden at UPM facilities under environmental conditions and regular watering. At the beginning of the experiment, trees were about 2.5 m high and the diameter at breast height was 5 to 7 cm. Wounding was performed by making two rectangular windows 10 cm high and covering half of the stem circumference in each tree, placed on opposite sides of the main trunk and with a vertical separation of three times the window height between them. We removed bark, phloem, vascular cambium and first layers of xylem (Online Resource Figure S1). Wounds were produced on April 9th, when cambial activity was ongoing, and sampling was carried out at three time points after wounding: (i) one week after wounding, a frame of tissue from the margins of both wounds was collected on April 16th (R0; this sample was analogous to H1 sample from Chano et al. (2017) ; (ii) the newly emerged resprouts close to wound margins, and surrounding tissue (phloem, cambium and first layers of xylem) were harvested on June 3th (R1), when trees were still forming early wood (Chano, Lopez de Heredia et al. 2017); and (iii) other new resprouts that were left from R1 and continued developing were collected on June 25th (R2). Control samples were also harvested from remote branches of the same trees, and all the collected samples were frozen in liquid nitrogen and stored at -80 ºC until further processing.

RNA isolation and microarray analysis
RNA extraction from each individual sample was performed as described in Chang et al. (1993). Total RNA was then purified with the RNeasy Plant Mini Kit (Qiagen, Dussldorf, Germany), and quantified with Nanodrop ND-1000 (Thermo Scientific, MA, USA). RNA from the three biological replicates of each sampling point (R0, R1 and R2) was hybridized with a two-color microarray (Agilent, USA) designed with a set of 15,266 contigs selected from a de novo assembled transcriptome representing cambial activity in P. canariensis, as described in Chano, Lopez de Heredia et al. (2017). In addition, 2303 contigs from libraries representing apical shoots in Pinus pinea were also included. The analysis of the output included background correction and normalization of expression values using LIMMA (Linear Models for Microarray Data; Smyth and Speed 2003), and selection of Differentially Expressed Genes (DEGs) by means of "Rank Products" algorithm available for R Bioconductor (Gentleman et al. 2004). Those transcripts showing a FC value (Fold Change, defined as the logarithmic ratio between experimental and control conditions for each biological replicate and sampling time) over 2 (overexpressed) or below − 2 (underexpressed), with a significance level p-value corrected by False Discovery Rate (FDR) below 0.05, were considered as DEGs. Expression values were clustered using the pheatmap package in R (Warnes et al. 2015). Microarray data was uploaded to the Gene Expression Omnibus database (GEO; https://www.ncbi. nlm.nih.gov/geo; Accession number GSE205852). By using Omicsbox following the Blast2GO methodology (Conesa and Götz 2008), we also performed an enrichment analysis based on Fisher's exact test of GO terms associated to these DEGs and using the whole P. canariensis' transcriptome as background.

Validation of expression profiles by qPCR
Validation of microarray expression data was performed using quantitative Real Time PCR (qRT-PCR) analysis of 11 DEGs (Table 1) covering the main profiles obtained from the analysis. The gene specific primers were designed with the Primer 3 webtool (Untergasser et al. 2012) with melting temperatures between 60 and 65 ºC, and producing amplicons ranging between 80 and 120 bp. We also included the Ri18S as housekeeping gene for CT normalization between different qRT-PCR runs ( Table 1). First strand cDNA synthesis was performed using the SuperScript™ III reverse transcriptase (Invitrogen™) following the protocol from the manufacturer, and the SsoFast™ EVAgreen® Supermix (Biorad, USA) was used as fluorescent DNA stain. The qPCR reactions were performed in a CFX96™ Real-Time PCR Detection System (Biorad, USA), following a thermal profile of 95 ºC for 3 min, 40 cycles of 95 ºC for 10 s, and 60 ºC for 10 s. Expression values were obtained using the Δ-Δ-CT method with PCR efficiency correction (Pfaffl 2001), and the statistical analysis of qPCR data included Shapiro-Wilk test for distribution, and Bartlett test for homoscedasticity. Statistical differences between groups were inferred following ANOVA. As well, correlation between microarray and qPCR data was assessed by means of Pearson's correlation test.

Differentially expressed genes during resprouting in response to wounding
Transcriptomic analysis of resprouting induced by mechanical wounding in P. canariensis was performed by means of a two-color cDNA microarray (Agilent, USA). Samples were taken during three stages over the responsive period, namely R0, R1 and R2, and transcripts with a Fold Change over 2 and below − 2 and significance level adjusted by FDR below 0.05, at least at one sampled date over the time series were selected for downstream analysis (Fig. 1). In total, 1164 DEGs were identified, of which 903 were found in R0, 278 in R1, and 261 in R2. Comparative analysis by groups resulted in 63, 67 and 76 DEGs that were exclusively shared between R0 and R1, R1 and R2, and R0 and R2, respectively, and 36 that were common for the three sampled dates (Fig. 2). On its side, the GO enrichment analysis by using Omicsbox did not yield any statistically enriched term.  between both qPCR and microarray hybridization, the general profiles obtained in this study were validated (Fig. 4).

Clustering of DEGs by Gene expression profiles
According to their transcriptional variations during the analyzed period, six co-expressed gene clusters were identified ( Fig. 3; clusters 1 to 6). Cluster 1 showed genes that were mainly repressed at R0, and roughly recover basal expression at R1 and R2. Genes included in cluster 5 presented a similar expression profile, although showing slight overexpression at R2. On the other hand, clusters 2 and 3 showed induced genes at R0, declining at R1 and R2. Cluster 4 included genes mainly repressed during the whole response. And finally, genes grouped in cluster 6 showed overexpression at R1 and R2 from a basal level at R0. The complete list of DEGs can be found in Online Resource Table S1, while a selection of 134 DEGs is shown in Table 2, classified by their functional role, and including the top-hit description from BLASTx, the FC values, and significance level after the statistical analysis.
Moreover, we selected 11 DEGs for validation of the expression profiles resulted from microarrays (Table 1), including 4 genes with key role in cell wall formation and lignin deposition (pectinesterase-like, CeSA-like, CCoAOMT-like and PAL-like proteins), coding for transcription factors involved in developmental processes (MYB46-like, ATHB13-like, NAC2-like, WRKY51-like, and YAB5-like), and encoding proteins related to stress response (EXORDIUM-like, and Major Allergen Pru AR1like proteins). According to the correlation coefficients

Resprouting stages
A total of 436 genes were differentially expressed R1 and R2, when the new epicormic shoots arose from stem, significantly up-or downregulated. From the latter, 244 DEGs showed a significant repression during both R1 and R2, or at least at one of these stages. On the other hand, cluster 6 included 21 DEGs overexpressed at R1 and, to a greater extent, R2, with high FC values ranging from 8 to over 100. Among these DEGs, we should mention those transcription factors belonging to the HD Zip Class IV family, such as PROTODERMAL FACTOR1-like (PDF1; Ppnisotig08430) and MERISTEM LAYER1-like proteins (ML1; Ppnis-otig05462), or the homologous of the HOMEODOMAIN GLABROUS 11-like protein (HDG11; Ppnisotig07853). Moreover, the transcript Ppnisotig04954, encoding for a P. canariensis orthologous member of the GROWTH-REGULATING FACTOR (GRF) family, showed increasing FC values in R1 and especially in R2. Additionally, a found in cluster 2 to be overexpressed at R0. Moreover, several TIFY-like transcription factors were also induced at R0, including one transcript coding for a TIFY 10a-like protein (Contig06361) found in cluster 3 highly overexpressed.
Conversely, some genes showed repression at R0, as those included in clusters 1 and 5. In this regard, we found genes coding for catalytic subunits of the cellulose synthase (f.i., Contig00654, -20,721, and − 41,981), members of the CAZymes superfamily (f.i., Contig05066,436,778 or Ppnisotig01389), a COBRA-like protein (Con-tig01405) or genes involved in phenylpropanoid biosynthesis (Contig02447 and − 06476). In addition, transcription factors of the MYB (Contig02588, -12,050, and − 12,416) and bHLH families (Contig12421 and − 15,411) were also found in cluster 1 to be repressed at R0, together with the homologous of Arabidopsis ATHB14 (Contig01739). Also in cluster 1 we found a homologous of auxin efflux carrier component 1c-like (PIN1C), encoded by Contig19183, that was repressed in R0.    overexpression at R1 and R2 was detected for the sequence Ppnisotig05961, which is coding for an auxin-responsive protein IAA33. Finally, in cluster 2 we also found a gene encoding for a ATHB13-like protein (Contig20304), a YUCCA4-like, and three members of the PLATZ-like transcription factor family (Contig13870, -24,637, and Ppnis-otig07889), all of them induced at R2. The Venn diagram in Fig. 5 illustrates a comparative analysis by groups between DEGs in the resprouting process FLOWERING-PROMOTING FACTOR1-like (FPF1) protein, encoded by Contig03401, had a similar expression profile than GFR during R1 and R2.
Among the transcription factors found in cluster 5, Contig13239, encoding for a YABBY5-like protein, was significantly overexpressed in phase R1, while a homologous of LEAFY (LFY; Contig03423) was induced at R2, as well as a PHD (plant homeodomain) finger protein (Con-tig03217). Also in cluster 5, repression at R0 followed by   during R1/R2 (including 48 overexpressed DEGs in at least one sampling date), and 91 were expressed just in R1/R2 (where 32 DEGs were induced at some point).

Discussion
Microarray hybridizations of RNA samples covering the resprouting response to wounding yielded a total of 1164 DEGs, of which 903 DEGs were found in R0, the largest set, while R1 and R2 included 278 and 261 DEGs, respectively. In addition, these 1164 DEGs were grouped in 6 clusters according to their expression profiles along the response. Considering these significant differences, two main periods were covered by the sampling design, (i) R0, with DEGs and DEGs found during the seasonal apical growth in the Canary Island pine. From the 1164 resprouting-DEGs, 571 were shared with apical growth (49.05%), of which 295 DEGs were underexpressed during the whole response and 276 overexpressed at some point of the process. From the latter, 108 DEGs were induced during the resprouting phases R1 and/or R2. Moreover, 170 DEGs were shared between R1, R2 and apical development, but not with R0. Correspondingly, 6599 DEGs were exclusive for apical growth, while another set of 593 genes were resproutingexclusive DEGs (included as DEGs in R0, R1 and R2, but not in apical growth). From them, 407 were uniquely expressed at R0 (34.97% of total), including 188 induced and 219 repressed genes. On the other hand, 95 resprouting exclusive DEGs (16.84% of total) were expressed also  1997). Moreover, the peroxidases and glutathione-S-transferase-like proteins found in clusters 2 and 3 are involved in ROS (Reactive Oxygen Species) detoxification during the response to pathogenesis (Gunnar Fossdal et al. 2001), and in the anti-oxidative plant defense (Meyer et al. 2008). Consistently, transcription factors with active roles in response to stress were also overexpressed during R0, like NACs and WRKYs (Zhang and Wang 2005;Hu et al. 2010). On the other hand, several genes involved in xylogenesis, like CesA or CAZymes, or in the phenylpropanoid biosynthesis, showed repression at R0, together with transcription factors related to cell wall biosynthesis like MYBs and bHLH, or the HD Zip Class II family member ATHB14/ PHABULOSA (PHB) (Ariel et al. 2007). The latter has been described as mediator during adaxial/abaxial polarization in ovule primordia (Sieber et al. 2004), the establishment of the SAM and apical bilateral symmetry (Prigge et al. 2005), and the perception of radial positional information in leaf primordia (McConnell et al. 2001). Interestingly, several of these genes were related to earlywood formation during xylogenesis in the Canary Island pine (Chano, Lopez de Heredia et al. 2017), and its repression seems to be a major trend during the rearrangement of the xylogenesis transcriptional program in the immediate response to wounding, as described also in Chano et al. (2017).

Resprouting stages
Some of the induced DEGs at R1 and/or R2, like those included in clusters 2, 5 and 6, were presumably involved in developmental processes and de novo organogenesis. For instance, YUCCA4-like transcription factor (cluster 2, mainly related to the immediate response to wounding, including 728 DEGs (62.54% of total) that were exclusive for this time, and (ii) R1 and R2, with genes presumably involved in the emerging and de novo development of epicormic shoots. Additionally, even though GO enrichment analysis yielded no statistically enriched terms, we were able to identify terms related to resprouting or outgrowth processes that were associated to R1/R2 DEGs (Online Resource Fig. S2). Thus, we found remarkable GO terms related with developmental activity, such as "growth", "anatomical structure development" or, more specifically, terms such as "shoot system development" and "flower development", indicative of de novo lateral organogenesis.

Immediate response
In accordance with our findings during the immediate response described in the course of the wound-healing process (Chano, Lopez de Heredia et al. 2017), R0 included genes related to pathogenic and biotic stress response, such as the disease resistance response protein PI206, or the putative PR-4-like protein, both found during the infection responses in Pisum sativum L. (Culley et al. 1995) and in Solanum tuberosum Poepp. Ex Walp. (Stanford et al. 1989), respectively. We also found the overexpression of chitinases and endochitinases-like proteins at this stage, with important roles in plant defense against fungi and insects (Seo et al. 2008), PAL-like proteins, involved in the phenylpropanoid biosynthesis pathway and in plant defense (Ohl et al. 1990;Kervinen et al. 1998;Kim and Hwang 2014), antimicrobial peptide 1-like with antifungal and antibacterial properties (Castro and Fontes 2005), defensin-like protein involved in defense-related processes and plant development (Tam et al. 2015) and a major allergen pru ar1 homologous which is at R2) has been reported to be implicated in the repression of trichome outgrowth in Arabidopsis (Nakamura et al. 2006). Interestingly, these three co-expressed HD Zip class IV members (PDF1-, ML1-, and HDG11-like) were downregulated along the apical seasonal growth in the Canary Island pine (Chano et al. 2021). In addition, FPF1 and GRF1 (cluster 6, both overexpressed at R1 and R2) have also been described to participate in flowering (Kania et al. 1997;Hoenicka et al. 2012), as well in lateral organogenesis and multiple developmental processes (Omidbakhshfard et al. 2015). GRFs and part of the microRNA regulator machinery of GRFs expression were also reported to be involved in floral organ development, determining sepal/petal identity (Pajoro et al. 2014).
In the stem, epicormic resprouting is developed from precursor meristematic cells or directly by the SAM in already preformed buds. The SAM is the meristematic tissue located at the tip of lateral branches and in the apex of the main trunk (Murray et al. 2012). Therefore, many subprocesses involving cell differentiation, proliferation and growth from SAM are common between epicormic resprouting and seasonal apical development. The identification of those genetic regulators triggering epicormic resprouting in response to stress, and therefore not involved in seasonal apical growth in the Canary Island pine (Chano et al. 2021), yielded 32 DEGs showing overexpression at one or both resprouting stages (R1 and R2). Among these 32 DEGs, we found the overexpression of ATHB13-like protein in R2, which has been suggested to be involved in regulation of cotyledon and leaf development in Arabidopsis (Mattsson et al. 2003;Henriksson et al. 2005). This gene was also induced during the immediate response to wounding in the Canary Island pine (Chano et al. 2017). In A. thaliana this gene was reported to have a negative role in stem elongation, while it is usually expressed in tapetum development during pollen formation (Ribone et al. 2015), and upregulated in response to abiotic stress conditions and pathogenesis .

Concluding remarks
Many details should be considered during the activation of dormant buds, including the initiation of epicormic resprouting, and the elongation of these newly formed sprouts. The big picture of lateral organogenesis and the central role played by the SAM along the whole process is complex, with many transcription factor families implicated by activating or deactivating expression of other gene families, with specific roles in cell wall synthesis, cell wall loosening, or microtubule organization, or the maintenance of meristematic activity of SAM's initial cells, among other processes (Traas 2019). overexpressed in R2) has been described to play an important role in auxin increase during floral organogenesis in Arabidopsis thaliana L. (Cheng et al. 2006;Mungía-Rodríguez et al. 2020). In P. canariensis, this gene was reported to participate during primary growth in the early stages of apical shoot development (Chano et al. 2021). Additionally, YABBY5-like protein (cluster 5, significantly induced at R1) is implicated in the regulation of lateral organ growth in seed plants (Yamada et al. 2011), as well as LEAFY-like transcription factor (cluster 5, induced at R2), a plant-specific transcription factor controlling flower development from axillary meristems (Schultz and Haughn 1991;Weigel et al. 1992). These two transcription factors are relevant candidates to be involved in the initiation and elongation, respectively, of resprouting in conifers, as they differ in the moment of their activity. YAB5 was reported to participate in the fate and maintenance of meristem activity during lateral organogenesis in rice (Tanaka et al. 2012). As well, YAB5 acts redundantly in Arabidopsis with FIL/YAB3 and YAB2 in the polar establishment of abaxial surface, activating laminar programs, repressing SAM and forming the marginal domain in leaves (Sarojam et añl. 2010). Moreover, the homologous YAB5 in Cabomba coraliniana was expressed in the abaxial tissues of the distal leaf primordium, playing an important role during the proliferating cell process, as well in the floral bud development. In vegetative shoots, YAB5 displays a similar expression pattern than FIL/YAB3, expressed in marginal and abaxial tissues in lobe primordial and procambial strands (Yamada et al. 2011). Regarding LEAFY, it has been found recently in the lycopsid genus Isoetes to be involved in the development of both reproductive and vegetative tissues during lateral organogenesis (Yang et al. 2017). These two genes were found to be expressed during the initiation of apical growth in the Canary Island pine but repressed at later stages (Chano et al. 2021).
Furthermore, PDF1 (cluster 6, overexpressed at R1 and strongly induced at R2) is a proline-rich cell wall protein expressed exclusively in the L1 layer of shoot meristems in Arabidopsis, controlling cell differentiation in the epidermis of new buds (Abe et al. 2003). In Gossypium spp. it was found that the GbPDF1 protein has an important role during fiber initiation (Deng et al. 2012). Moreover, the functional role of PDF1 was found under regulation of two HD Zip Class IV members in Arabidopsis (Abe et al. 1999(Abe et al. , 2003. One of them is PROTODERMAL FACTOR2 (PDF2), which has not been identified as DEG in this work, while the other is the Arabidopsis thaliana ATML1 gene, found to be co-expressed together with PDF1 (also in cluster 6, overexpressed at R1 and R2). This gene was first suggested to be expressed only in the first layer (L1) of the SAM (Lu et al. 1996). In addition, HDG11 (cluster 6, overexpressed close.

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
The molecular basis of the epicormic shoot production process is still poorly studied in gymnosperms. This study sheds some light on this process in the Canary Island pine, wherein resprouting can be induced after wounding. Indeed, resprouting has been suggested to occur following a reduction of cambial activity (Bachelard 1969;Ishii 2011), and mechanical wounding was also suggested to cause cambial growth cessation in P. canariensis (Chano et al. 2015(Chano et al. , 2017a. The examination of transcriptional changes during this process allowed us to determine genes required for cell signaling, cell identity and cell proliferation sub-processes, including indumentum development, laminar outgrowth, and bud emerging. Our data provide evidence in favor of the resprouting activity of YABBY and LEAFY transcription factors, and members of the class IV of the Homeodomain Leucine Zipper such as PDF1, ML1 and HDG11. As well, other transcription factors involved were GRF1 and FPF1. The implication of such transcription factor families, with important roles in leaf and inflorescence development also in flowering plants, suggests underlying homologies between cell proliferation and tissue differentiation patterns in many lateral organogenesis processes.
These are the first insights into the transcriptomic basis of resprouting in conifers. Many of the genes reported here showed significant overexpression during resprouting as well as intense transcriptomic activity at the beginning of apical growth, suggesting non-exclusive mechanisms for resprouting in P. canariensis, at very least during development of epicormic shoots. Nevertheless, we cannot exclude an exclusive genetic control of resprouting initiation: genes such as ATHB13 or PHD-finger, together with some nonannotated sequences, were not found during seasonal shoot development. Thus, further research is needed to disentangle the functionalities of these candidate genes for comprehensive analysis of the resprouting capability in conifer species.
Authors contribution AS and CC designed the experiment. AS and VC carried out the experiment and sample collection. VC performed lab work, data analysis, and wrote the first version of the manuscript. AS and OG edited and reviewed the final version. All authors have read and approved the final manuscript.
Funding This work was funded via projects AGL2009-10606 (Spanish Ministry of Science and Innovation), SPIP2014-01093 (Spanish National Parks Agency, Ministry of Agriculture) and FP7-KBBE-2011-5 289841 (European Commission, Seventh Framework Programme). Open Access funding enabled and organized by Projekt DEAL.

Competing interests
The authors have no competing interests to dis-