Gene expression profiling during adventitious root formation in carnation stem cuttings

Sequencing and transcriptome assembly supports the discovery of novel genes expressed
in the stem cutting base

In a recent study 26] we characterized AR formation in a collection of 10 carnation cultivars. The 2003 R 8 and the 2101–02 MFR cultivars have been chosen for further studies due to their differences in rooting
performance and in their differential response to a mild auxin treatment during rooting
(Fig. 1b). The bad-rooting behaviour of the 2003 R 8 cultivar, which was mostly caused by a delay in AR initiation, was partially restored
by exogenous auxin application.

Several cDNA libraries prepared from stem cutting bases of mock-treated and auxin-treated
samples at particular time points during adventitious rooting (0, 6, 24 and 54 hAP)
were sequenced (see Methods; Fig. 1a). As a result, 3,683 million of raw reads were obtained. The amount of expression
data generated in our study had the potential to transform the boundaries and extent
of previous feature annotations in the carnation genome 27]. Genome-guided assemblies were performed with the purpose of serving as sequence
evidence for a genome re-annotation. Thus, our updated, evidence-based annotation
comprised 59,396 transcripts, corresponding to 57,641 genes, with an average length
of 2,856 bp (Table 1). We were able to merge exons from genes that had been previously predicted as separate
coding sequences 27]; these merges resulted in a more complete or contiguous annotation for 394 genes
and their corresponding transcripts (see Additional file 3). In order to quantify the improvements of the new annotation in relation to the
former annotation, we compared the proportion of reads mapping to each of them and
also generated a number of descriptive statistics (Table 1).

Table 1. Comparison between the previously published annotation and the updated genome annotations

Time-dependent comparison of the auxin treatment identifies 1,286 differential expressed
genes (DEGs) in response to the auxin stimulus

We then tested the effects in gene expression of factors like cultivar, treatment
and time point using different models and contrasts (Table 2). Three questions were addressed: i) for which genes does the cultivar factor have a significant effect? (Table 2, test 1), ii) which genes change their expression over each pair of time points? (while accounting
for cultivar-specific effects; Table 2, test 2–13) and iii) how does the auxin treatment affect gene expression distinctively at each time point
for each cultivar? (Table 2, test 14–21). To reduce the complexity of the model, in this last case, we made subsets
of samples belonging to each cultivar and estimated the parameters separately.

Table 2. Differential expression tests

In all cases, expression models were fitted to our time-course study by treating each
time point as a different “experimental group”, even though the inherent ordering
and spacing provided by time points is ignored then.

To investigate the dependencies between treatments and time, explicitly addressing
the question of when a gene is differentially expressed, we modeled the interaction
between time and treatment as a covariate. Of 57,641 genes, we filtered out genes
that were not expressed (0 counts), resulting in a total of 37,936 genes tested for
the subset of cultivar 2101–02 MFR and 37,849 for the cultivar 2003 R 8 subset. The factorial analysis identified a total of 1,286 distinct genes as differentially
expressed between auxin-induced and control cuttings over different time points (Table 2, test 14–21). Most auxin-related expression changes took place in the initial time
points (0 hAP vs. 6 hAP). Among them, DEGs of 2101–02 MFR were associated (Fisher exact test) to functions like photosynthesis (GO:0015979;
P 0.001) and chlorophyll binding (GO:0016168; P 0.001). As for the same comparison in the cultivar with poor rooting performance,
2003 R 8, the functions associated to DEG were translational initiation factor activity (GO:0003743,
P?=?0.0025), and negative regulation of signal transduction (GO:0009968, P?=?0.0024), among others.

Clustering of time-course expression profiles reveals co-expression of functionally
related genes

To get some insight into the specific pathways regulated at different time points
during AR formation in the two cultivars studied, we performed a GO-enrichment analysis
for the sets of DEGs shown in Fig. 2a. In the 2003 R 8 cultivar, the GO categories “protein amino acid phosphorylation” (GO:0006468; 105
genes; P??0.001) and “transmembrane transport” (GO:0055085; 50 genes; P??0.001) were specifically and significantly enriched at 6–24 and 24–54 hAP, respectively.
Interestingly, the “auxin-activated signalling pathway” (GO:0009734; 11 genes; P??0.001) category was found significantly enriched among DEGs shared between 0–24
hAP in this cultivar. Conversely, in the 2101–02 MFR cultivar, the GO category “hormone-mediated signalling pathway” (GO:0009755; 58 genes;
P??0.001) was specifically enriched at 0–6 hAP. Moreover, the “auxin-activated signalling
pathway” (GO:0009734; 34 genes; P??0.001) and “cell cycle” (GO:0007049; 79 genes; P??0.001) categories were found significantly enriched for DEGs shared between 0 and
24 hAP. The GO categories “glucose catabolic process” (GO: 0006007; 41 genes; P??0.001), and “cellulose biosynthetic process” (GO:0030244; 22 genes; P?=?0.002) were significantly enriched specifically at 6–24 hAP, whereas at 24–54 hAP
the most significant GO-enrichment was found for genes assigned to the “response to
stress” (GO:0006950; 46 genes; P?=?0.007) category.

Fig. 2. Differentially expressed genes (DEGs) over time during AR formation. a Venn diagram illustrating DEGs (P??0.05, Benjamini-Hochberg correction) in the stem cutting base of the 2003 R 8 and 2101–02 MFR cultivars over time. b Gene ontology (GO) classification of genes whose expression changes more than 1.25
log-fold over the time-course with respect to biological process (BP), cellular component
(CC) and molecular function (MF). Asterisks indicate significant enrichment of genes
(P??0.05)

To obtain a more general view of the functions involved in the early stages of AR
formation, we transformed the GO functional annotation into its cut-down version (Plant
GO slim) for the 14,554 DEGs in the 2101–02 MFR cultivar between 0 hAP and 6 hAP (Fig. 2b). The Biological Process (BP) classification of DEGs highlighted a significant enrichment
(P??0.001) for the following GO categories: “biosynthetic process” (GO:0009058) and
“transport” (GO:0006810). Among the Cellular Component (CC) categories, “plastid”
(GO:0009536) and “plasma membrane” (GO:0005886) were the most significantly enriched
ones. In terms of Molecular Function (MF), a significant enrichment was found for
genes at the categories “transferase activity” (GO:0016740) and “nucleotide binding”
(GO:0000166). These results indicated that large expression changes are taking place
at transcriptome level in the stem cutting base during the initial stages of AR formation
in this cultivar.

Validation of expression of some of the genes detected during AR formation

The reliability of our transcriptome profiling dataset was validated by examining
the expression of selected genes by using qRT-PCR and by comparing them to the normalized
data obtained in the RNA-Seq analysis (see Methods). We found highly significant and
positive correlations between qRT-PCR and RNA-Seq results for both cultivars in all
time points and treatments (Fig. 3a-b). Additional statistical analysis revealed that the variation observed between qRT-PCR
and RNA-Seq results depended largely on the expression levels of the studied genes
(Additional file 4: Figure S2). Thus, for genes with very low or very high numbers of RNA-Seq reads,
the qRT-PCR validation was less accurate. Representative examples of the results obtained
for genes with contrasting expression profiles are shown in Fig. 3c-f. While the expression levels of Dca5879 were only varying over time (Fig. 3c-d), those of Dca29160 were also depending on cultivar and the auxin treatment (Fig. 3e-f).

Fig. 3. qRT-PCR validation of RNA-Seq results. a–b Correlation analysis of qRT-PCR and RNA-Seq data from selected genes in the 2003 R 8 (a) and 2101–02 MFR (b) cultivars. Each dot represents the relative expression data for a given gene and
a given sample. c–f Bars represent the mean of the relative expression level of qRT-PCR (black) or RNA-Seq
(grey) data relative to mock-treated samples at 0 hAP. Error bars indicate the standard
deviation for the mean data shown

Gene set enrichment analysis in the 2003 R 8 cultivar

As the number of DEGs between auxin-treated and mock-treated samples in the 2003 R 8 cultivar was scarce (see above), we did not distinguish between treatments in this
cultivar when performing a GO-enrichment analysis using STEM (see Methods). The expression
of 7,341 genes was found to be specifically regulated during AR formation in the 2003 R 8 cultivar and 4,599 of these genes were clustered along different expression profiles
and further classified into four major groups based on their expression pattern between
0–6 hAP and 6–54 hAP: DownDown (DD), DownUp (DU), UpDown (UD), and UpUp (UU) (Fig. 4a).

Fig. 4. Model profiles identified during AR formation in carnation cultivars. a A combined auxin- and mock-treated-dataset is shown for the 2003 R 8 cultivar. b The mock-treated dataset for the 2101–02 MFR cultivar. Profiles were classified into four groups: DownDown (DD), DownUp (DU),
UpDown (UD), and UpUp (UU) and further ordered based on their P-value (bottom left-hand
corner). The number of genes assigned to each profile is shown in the top left-hand
corner

1,286 genes were gradually repressed (DD group) during AR formation in these conditions.
Some of the most significantly-enriched GO categories within this group were “response
to auxin” (GO:0009733; 15 genes; P??0.001) and “ion transport” (GO:0006811; 46 genes; P??0.001). We assigned 1,381 additional genes to the DU group. About two-thirds of
these genes showed an early repression and later became upregulated above their expression
at 0 hAP (labelled as “a” in Fig. 4a). The remaining genes in this group, which were quickly downregulated and whose levels
were more-or-less restored to initial levels at later time points (labelled as “b”
in Fig. 4a), showed enrichment in “photosynthetic membrane” (GO:0034357; 30 genes; P??0.001) encoding genes. Some of the most significantly-enriched GO categories within
the DU group as a whole were “cell wall organization or biogenesis” (GO:0071554; 51
genes; P??0.001), “cytoskeleton” (GO:0005856; 46 genes; P??0.001), and “cellular carbohydrate metabolic process” (GO:0044262; 63 genes; P??0.001). Another 790 genes showing a biphasic response were classified into the
UD group. Finally, 1,142 genes within the UU group were ranked for GO enrichment:
“cell division” (GO:0051301; 32 genes; P??0.001), “microtubule” (GO:0005874; 52 genes; P??0.001), and “cell wall organization or biogenesis” (28 genes; P??0.001) among others. Interestingly, the expression of most genes included within
the “cell division” and “cell wall organization or biogenesis” categories peaked after
6 hAP in agreement with the timing of cell cycle re-activation in the cambium observed
for this cultivar, as it is shown later.

Expression profiling in the 2101–02 MFR cultivar and in the response to auxin

12,525 genes were identified in the 2101–02 MFR cultivar as being specifically expressed during AR formation without exogenous auxin
treatment (Fig. 4b). 3,586 genes were assigned to the DD group where one of the most significantly-enriched
GO categories was “protein serine/threonine kinase activity” (GO:0004674; 175 genes
[13.8 % of the protein serine/threonine kinase encoding genes with dynamic expression
profiles (EGs)]; P??0.001). 1,737 genes showing a biphasic response were classified into the DU group.
Those upregulated at later time points (labelled as “a” in Fig. 4b) encoded proteins enriched in “cytoskeleton” (GO:0005856; 41 genes; P??0.001) and “cell division” (17 genes; P??0.001). Similarly to that found previously for the 2003 R 8 cultivar, the “photosynthetic membrane” (35 genes; P??0.001) category was found enriched among genes whose expression levels were restored
to basal levels (labelled as “b” in Fig. 4b). Among the biphasic genes that were assigned to the UD group (1,385), one of the
significantly-enriched GO categories found was “carbohydrate derivative metabolic
process” (GO:1901135; 52 genes; P?=?0.002). Finally, 2,103 genes were included within the UU group, where the most
significantly-enriched GO categories were “cellular carbohydrate metabolic process”
(69 genes; P??0.001) and “cell wall organization or biogenesis” (54 genes; P??0.001). In addition, we found specific GO-enriched categories in profile P29 (Fig. 4b). On the one hand, enriched genes upregulated after 6 hAP (P29) encoded proteins
related to “microtubule” (38 genes; P??0.001), “cell division” (24 genes; P??0.001), and “regulation of cell cycle” (GO:0051726; 23 genes; P??0.001). On the other hand, genes encoding putative chromatin-related functions
such as “histone H3 lysine 9 methylation” (GO:0051567; 14 genes; P??0.001), or “DNA packaging” (GO:0006323; 14 genes; P??0.001) were also found significantly enriched.

Additionally, the expression of 9,645 genes was found specifically altered during
AR formation in the 2101–02 MFR cultivar after exogenous auxin treatment and 5,568 of these genes were significantly
clustered to different expression profiles and grouped as described above (data not
shown). We found a substantial overlap between EGs of auxin-treated and mock-treated
samples (61.2 % for the auxin-treated EGs and 79.5 % for the mock-treated EGs). Consistently,
no significant differences in the overall trends of EGs were found between auxin-
and mock-treated samples for the 2101–02 MFR cultivar (Additional file 5: Figure S3). However, for a small number of genes assigned to specific profiles in
mock-treated samples, we found some changes in their expression profiles after the
auxin treatment. About half of the EGs-encoding proteins belonging to “microtubule”,
“cell division” and “histone H3 lysine 9 methylation” were reciprocally assigned to
profiles P29 (UU) or P21 (DU) in mock-treated and in auxin-treated samples, respectively.
We also found that the expression levels of most genes assigned to the “cellular carbohydrate
metabolic process” category were complementary at earlier time points in auxin-treated
vs. mock treated samples.

Comparative transcriptome profiling of AR formation between carnation cultivars for
selected GO categories

To identify genes whose expression correlates with the different stages of AR formation
and that could be used as markers, we selected EGs belonging to the GO categories
“cell division” and “response to auxin” from the different samples studied (see above).
We then built heat map representations from log expression data for all these genes
(Fig. 5). On the one hand, the expression of some genes encoding proteins related to cell
division were clustered together along the time point series with higher expression
at earlier time points, independently of cultivar and treatment (Fig. 5a). Interestingly, genes encoding mitotic cyclins (A-type and B-type) showed a peak
of expression between 24 hAP and 54 hAP (Fig. 5a), which is in light with the cellular changes observed in the stem cutting base during
AR rooting (see next section). Noteworthy, the majority of these genes moderately
respond to the auxin treatment by increasing their expression levels at earlier time
points (Fig. 5a). On the other hand, a small number of genes displayed contrasting expression profiles
between cultivars. Examples for the latter are the Dca37619 and Dca642 genes, which are respectively upregulated and downregulated in 2003 R 8 compared to the 2101–02 MFR cultivar (Fig. 5a).

Fig. 5. Analysis of expression of transcripts related to cell division (a) and response to auxin (b) during AR formation. Heat map drawing and clustering was done as described in Methods

Considering the expression profiles of genes assigned to the “response to auxin” category,
the effect of the auxin treatment was quite small irrespective of the cultivar, and
was mainly restricted to earlier time points (Fig. 5b). However, we found striking differences in the expression of a few of these genes
between cultivars, such as Dca1208 and Dca39239, which makes them candidates for further studies to analyze their role in the differential
response in auxin-mediated AR initiation between these two cultivars.

Cellular changes in the stem cutting base during AR formation reflects the effect
of the auxin signal

In a previous study we found that cell divisions within the cambial region of the
stem cutting base took place between 12 hAP and 24 hAP in a good-rooting cultivar
used as a reference 43]. We next characterized the cellular changes occurring within the cambium region in
the stem cutting base of the 2003 R 8 and 2101–02 MFR cultivars both in mock- and auxin-treated samples to understand the differential
responses observed in these two cultivars during AR formation. Although the cambial
ring of the 2003 R 8 cultivar displayed a very organized cellular pattern at 0 hAP, we found that some
regions within the cambium displayed subtle tissue disorganization (Fig. 6a and Additional file 6: Figure S4). Interestingly, we observed an increase in the number of disorganized
regions within the cambial ring at later time points, which could reflect local activation
of cell divisions (Fig. 6b). In addition, stem cutting bases of the 2003 R 8 cultivar treated with auxin contained an increased number of these disorganized regions
already at 0 hAP (Fig. 6c and 6d). In the 2101–02 MFR cultivar, we observed a higher frequency of small cell clusters within the cambial
ring at 0 hAP, which at later time points developed as large clusters of meristematic
cells with a disorganized internal structure (Fig. 6f). In most cases, these cell clusters appeared juxtaposed but physically isolated
by collapsed neighbouring cells (arrowheads in Fig. 6f). In auxin-treated samples however, cell clusters became apparent already at 0 hAP
(Fig. 6g), which is indicative of an early activation of cell division in the 2101–02 MFR cultivar. At 54 hAP, cell clusters were clearly evident and their numbers were higher
than in mock-treated samples (Fig. 6h).

Fig. 6. Morphological changes in the ultrastructure of the basal stem of carnation cuttings
during adventitious rooting. Light micrographs were taken from cross-sections of stem
cutting basal regions in the 2003 R 8 (a–d) and 2101–02 MFR (e–h) cultivars, either treated with mock (a–b, e–f) or with auxin (c–d, g–h) at different time intervals (0 and 54 hAP). ca, cambium; cc, cell clusters; pl,
phloem; xl, xylem. Scale bars: 50 ?m

To confirm our observations, we estimated some cellular parameters in the two contrasting
regions identified within the cambium (Additional file 6: Figure S4; see Methods). In the 2003 R 8 cultivar we observed that the cell division rates significantly differed between
organized and disorganized regions at the different time points studied, which seemed
not to be affected by the auxin treatment (Additional file 6: Figure S4B). These results suggested that auxin act as a trigger for cell division
within a certain population of responsive cambial cells. On the other hand, we observed
a significant increase in the number of cells within the cambium at later time points
for the 2101–02 MFR cultivar, both in organized and disorganized regions (Additional file 6: Figure S4C), which is indicative of a broad activation of cell division within the
cambium, as has been previously described for a good-rooting reference cultivar 43]. Interestingly, the division rate at a given time point was found unchanged irrespectively
of the auxin treatment (Additional file 6: Figure S4B). Taken together, these results suggested that auxin acts by promoting
divisions of quiescent cambial cells rather than by increasing the number of divisions
of already dividing cells, the former producing a net increase in the number of cell
clusters within the cambial ring.

Morphogenetic hormone levels in the stem cutting base during AR formation are correlated
with rooting performance

Several plant hormones play a crucial role in controlling AR formation, with auxin
and cytokinin playing opposite roles 1], 44]. In addition, wounding stimulates ethylene biosynthesis which is known to positively
influence AR formation in some species 45], 46]. We found high levels of endogenous indole-3-acetic acid (IAA) only in the 2101–02 MFR cultivar at 0 hAP, which were quickly downregulated to basal levels (Fig. 7), as previously described 43]. In addition, we found very low levels of trans-zeatin (tZ) in the stem cutting base of the 2101–02 MFR cultivar. In contrast, in the 2003 R 8 cultivar, tZ levels steadily increased during the time-course experiment (Fig. 7). Hence, the endogenous auxin/cytokinin ratio estimated as the proportion between
IAA and tZ levels was much higher in the 2101–02 MFR cultivar than in the 2003 R 8 cultivar for all time points, with the highest ratio found at 0 hAP. In addition,
the levels of the 1-aminocyclopropane-1-carboxylic acid (ACC) ethylene precursor were
higher in the 2101–02 MFR cultivar than in the 2003 R 8 cultivar (Fig. 7), which might reflect higher endogenous ethylene production in the 2003 R 8 cultivar.

Fig. 7. Changes in the concentration of some phytohormones in the basal stem of carnation
cuttings during adventitious rooting. Data for indole-3-acetic acid (IAA), trans-zeatin (tZ) and the ethylene precursor ACC are shown in ng/g fresh weight. Time refers
to hours after planting (0 hAP). Asterisks indicate significant differences (P??0.05) between treatments for a given sampling time

TF analysis

By modulating gene transcription at specific times and during specific processes,
TFs and their regulatory networks have important roles in development and stress response.
A number of transcripts that show significant expression changes during the time-course
experiment performed with STEM, encoded putative transcription factors that belonged
to 19 transcription factor families (http://planttfdb.cbi.pku.edu.cn/; Fig. 8). To identify the transcription factors that might regulate the differential rooting
responses observed between 2101–02 MFR and 2003 R 8 cultivars, we performed an enrichment analysis using Fisher’s T-test (see Methods). In the 2003 R 8 cultivar, a number of genes encoding C2C2-GATA transcription factors, such as Dca7186 and Dca56796, were significantly enriched and showed a clear downregulation of their expression
over time (Fig. 8). The WRKY transcription factor family was among the most highly downregulated transcription
factor genes irrespective of cultivar and treatment (Fig. 8). Of the 14 carnation genes encoding putative WRKY proteins that were shared between
cultivars and treatments, 10 were downregulated over the time-course experiment. Interestingly,
the expression of Dca28099, the putative carnation ortholog of the Arabidopsis PLETHORA5 gene 47] was found upregulated only in the 2101–02 MFR cultivar after auxin treatment, which emphasizes its value as an early marker for
adventitious root formation in this species.

Fig. 8. Expression of TF genes during AR formation. The number of upregulated (white bars)
and downregulated (grey bars) TFs is indicated for each family in 2003 R 8 (a) and 2101–02 MFR mock-treated (b) or auxin-treated (c) samples. Asterisks indicate significant enrichment of genes (P??0.05) within each TF family