Comparative Transcriptome Analysis Reveal Candidate Genes Potentially Involved in Regulation of Primocane Apex Rooting in Raspberry (Rubus spp.)

Raspberries (Rubus spp.) exhibit a unique rooting process that is initiated from the stem apex of primocane, conferring an unusual asexual mode of reproduction to this plant. However, the full complement of genes involved in this process has not been identified. To this end, the present study analyzed the transcriptomes of the Rubus primocane and floricane stem apex at three developmental stages by Digital Gene Expression profiling to identify genes that regulate rooting. Sequencing and de novo assembly yielded 26.82 Gb of nucleotides and 59,173 unigenes; 498, 7,346, 4,110, 7,900, 9,397, and 4,776 differently expressed genes were identified in paired comparisons of SAF1 (floricane at developmental stage 1) vs. SAP1 (primocane at developmental stage 1), SAF2 vs. SAP2, SAF3 vs. SAP3, SAP1 vs. SAP2, SAP1 vs. SAP3, and SAP2 vs. SAP3, respectively. SAP1 maintains an extension growth pattern; SAP2 then exhibits growth arrest and vertical (downward) gravitropic deflection; and finally, short roots begin to form on the apex of SAP3. The Kyoto Encyclopedia of Genes and Genomes enrichment analysis of SAP1 vs. SAP2 revealed 12 pathways that were activated in response to shoot growth arrest and root differentiation, including circadian rhythm—plant (ko04712) and plant hormone signal transduction (ko04075). Our results indicate that genes related to circadian rhythm, ethylene and auxin signaling, shoot growth, and root development are potentially involved in the regulation of primocane apex rooting in Rubus. These findings provide a basis for elucidating the molecular mechanisms of primocane apex rooting in this economically valuable crop.


INTRODUCTION
Raspberries (Rubus sp.) are an important economic fruit crop that grows in all temperate regions of the world. The raspberry fruit contains abundant polyphenol antioxidants, including anthocyanin pigments that are important ingredients of health products and can potentially prevent various human diseases (Skrovankova et al., 2015). Raspberries have a unique rooting characteristic, and its biennial shoots can grow several meters and root at the stem apex in autumn (Heslop-Harrison, 1959). In order to prevent excessive vegetative spread via rooting at stem apices, newly planted red raspberries often require staking to hold them upright. The apex rooting process of the stem apex can be divided into three successive stages: in stage 1, the stem exhibits an elongation growth pattern; in stage 2, elongation growth ceases and is followed by gravitropic curvature of the stem; and in stage 3, the root primordium differentiates from the stem apex and the root extends. After the formation and growth of new root at the stem apex, negatively geotropic shoots form at the rooting boss. New plants are generated in the next growth season. Thus, not only root primordium differentiation but also growth arrest and geotropism deflection are involved in formation of adventitious roots (ARs) at the stem apex, in contrast to adventitious roots derived from ordinary shoot cutting.
In Arabidopsis thaliana, the lateral root (LR) initiates from the primed LR founder cells in the xylem pole pericycle. After receiving particular signals, founder cells are activated and cell division is induced. The primordial LR begins to form and develops into the LR (Malamy and Benfey, 1997). Auxin plays an import role in this process (Benková et al., 2003;Dubrovsky et al., 2011;Sukumar et al., 2013;Villacorta-Martín et al., 2015). The polar transport of auxin from root tip to the aboveground part affects LR initiation, while auxin transport in the opposite direction influences LR germination (Reed et al., 1998;Casimiro et al., 2001). Many external factors affect LR initiation and growth by altering auxin distribution and polar transport, including gravity, bending, and mechanical stimulation De-Smet et al., 2007;Lucas et al., 2008). LR number decreased by ∼90% in transport inhibitor response (tir) 1 and tir1/auxin-signaling fbox (afb)2/afb3 triple mutants, indicating important roles for auxin receptor and auxin signal transduction in LR development (Kepinski and Leyser, 2005;Péreztorres et al., 2008). De novo organogenesis of root primordium can be divided into two steps (Hu and Xu, 2016): the initial transition from regenerationcompetent cells to root founder cells is accelerated by auxininduced upregulation of WUSCHEL-RELATED HOMEOBOX (WOX)11 and WOX12, after which WOX5/7 and its positive regulator WOX11/12 regulate root primordium initiation in de novo root organogenesis. These results suggest that early stages of root initiation are tightly regulated at the physiological as well as the genetic level. Most previous studies have focused on AR derived from the lower end of the stem, and it remains unclear how the stem apex of raspberries change from negative to positive geotropism and initiates rooting.
To address this issue, we analyzed gene expression in the stem apex during the early stages of adventitious rooting in primocane (i.e., the stem of the current season's growth) and floricane by Illumina HiSeq4000 analysis. We identified differentially expressed genes (DEGs) that may be involved in the regulation of AR formation in Rubus. Quantitative real-time (qRT)-PCR analysis was performed to validate the expression of some DEGs. Comparing gene expression patterns at three developmental stages of stem apex provided important insight into the molecular mechanisms of AR formation, for which we propose a model.

Plant Material
Red raspberry (Sijihong), which exhibits the typical spontaneous rooting from the stem apex of primocane, was used in the present study. The experimental orchard was located in Siping City, Jilin Province, China. At the start of spring growth, there are two classes of stem-primocanes and floricanes-that will develop during that season (Bailey and Howard, 1941). Primocanes are the annual main stems; these are sterile, and apex removal will induce the formation of three or more new sterile primocanes ( Figure 1A). Floricanes bear flowers and fruits, and apex removal will induce the formation of more floricanes that bear inflorescences either on the floricane itself or on their lateral branches ( Figure 1B). By the end of summer, the height of primocanes may be more than 1 m long, and the apex begins to exhibit geotropism before subsequently differentiating. The sampling date of primocane is divided into the following three stages according to morphological observations (Figure 2): SAP1 (stage 1, during which foliage on the tapering apex remains active); SAP2 (stage 2, when apex foliage withers, the apex becomes blunt, and its diameter increases dramatically); and SAP3 (stage 3, when the apex becomes chlorotic and short roots begin to form). Primocane samples (0.5 cm at stem apex) were collected between 9:00 and 10:00 a.m. on July 10, August 10, and September 10. At the same times, three floricane samples (SAF1, SAF2, and SAF3, 0.5 cm at stem apex) were collected as control materials ( Figure 2B). There were three biological replicates for each primocane and floricane sample. The samples were immediately frozen in liquid nitrogen and stored at −80 • C until use.

RNA Extraction, Library Construction, and Sequencing
In order to investigate global gene expression changes in the transcriptome during AR of the apex, we constructed 18 Digital Gene Expression (DGE) profiling libraries using the abovementioned samples. Total RNA was extracted from the dissected tissue using the RNA Easyspin Isolation System (Aidlab Biotech, Beijing, China). RNA was dissolved in diethylpyrocarbonatetreated water and stored at −70 • C until used for next-generation sequencing and quantitative real-time (qRT)-PCR validation. A total of 8 µg of total RNA was treated with DNase I and subjected to oligo (dT) magnetic bead adsorption to purify mRNA, which was fragmented by mixing with fragmentation buffer; the fragments were used for cDNA synthesis. Short fragments were purified and dissolved in elution buffer for end reparation and adenine addition, and then connected with adapters; suitable fragments were selected for PCR amplification. The Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and ABI Step One Plus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) were used for quantitative analysis and validation of sample libraries, and all 18 libraries were sequenced using Illumina HiSeq 4000 (San Diego, CA, USA) (Cheng et al., 2015). All raw transcriptome data were deposited in the Sequence Read Archive (https://www.ncbi.nlm. nih.gov/SRA, accessions Number: SRP105309).
FIGURE 1 | Diagrams of typical late season stems of "Sijihong" red raspberry (modified from Heslop-Harrison, 1959). The state of the buds is as indicated in the key. (A) Characteristics of primocane after terminal bud was removed (stem of the current season's growth). Shoot 2, 3, and 4 are newly induced lateral branches after the terminal bud 1 was cut off; (B) A 2-year-old floricane, terminal bud was cut during its first growth season.

Data Analysis and Mapping of DGE Tags
Sequencing reads of low quality that were adaptor-polluted or had a high content of unknown base (N) reads were filtered. De novo assembly with clean reads was carried out to obtain the final unigenes, which were annotated by querying against National Center for Biotechnology Information (NCBI) nucleotide and protein (NR) databases and Clusters of Orthologous Groups, Kyoto Encyclopedia of Genes and Genomes (KEGG), and SwissProt databases (Altschul et al., 1990). Blast2GO with NRwas used for Gene Ontology (GO) annotation (Conesa et al., 2005), and InterProScan5 was used for InterPro annotation (Quevillon et al., 2005). Unigenes that could not be matched to any database mentioned above were identified by ESTScan with BLASTpredicted coding DNA sequences as a model (Yu et al., 2016). In total, six pairs of DGE profiles of different sample libraries (SAP1 vs. SAF1, SAP2 vs. SAF2, SAP3 vs. SAF3, SAP2 vs. SAP1, SAP3 vs. SAP2, SAP3 vs. SAP1, where the former was used as control and the latter as the experimental group) were compared to assess gene expression changes during AR formation in Rubus. A strict algorithm for identifying DEGs between two samplesthe false discovery rate (FDR)-was used to determine the threshold P-value in multiple tests and analyses. FDR ≤ 0.001 and |log 2 Ratio| ≥ 1 were used as thresholds to identify DEGs. The possible functions of DEGs were determined by searching the GO database (http://www.geneontology.org/). Web Gene Ontology Annotation Plot (WEGO) was also used for GO classification of genes identified in each DGE library (Ye et al., 2006). To further characterize gene function, pathway enrichment analysis of the DGE results was performed by BLAST searches of the KEGG database (http://www.kegg.jp/kegg/). Cluster analysis of gene expression patterns was performed using Multi Experiment Viewer (http://mev.tm4.org/#/welcome). A Q ≤ 0.05 was selected as the threshold for significant enrichment of gene sets (Cheng et al., 2015).

qRT-PCR Analysis
DGE results were verified by qRT-PCR analysis. RNA samples used for qRT-PCR were identical to those used for DGE experiments. cDNA was synthesized using the PrimeScript RT Reagent kit (Takara Bio, Otsu, Japan) according to the manufacturer's instructions, and qRT-PCR was carried out using the SYBR Premix Ex Taq kit (Takara Bio) with gene-specific primers designed using Premier 6.0 software (Premier Biosoft International, Palo Alto, CA, USA) and synthesized commercially (Invitrogen, Carlsbad, CA, USA); the sequences are shown in Supplemental Table 1. The reaction mixture (final volume: 20 µL) contained 1.0 µL cDNA template, 2.0 µL 10 × reaction buffer, 0.5 µL dNTPs (10 mM), 1.0 µL each primer (10 µM), 10.0 µL 2× SYBR Green Supermix (Takara Bio), and 0.2 µL rTaq DNA polymerase (5.0 U/µL; Takara Bio), and the reaction was carried out on a Rotor-Gene 2000 thermocycler (Corbett Research, Sydney, Australia). Samples were prepared in triplicate in each experiment, and each biological sample consisted of two technical repeats. PCR conditions were as follows: 10 min at 95 • C, followed by 40 cycles of 15 s at 95 • C and 1 min at 60 • C. Amplification data were analyzed using real-time analysis software (Corbett Research). Relative expression levels of genes were calculated with the 2 − Ct method (Cheng et al., 2015).

Illumina Sequencing and Sequence Assembly
To detect global gene expression changes during AR formation of the apex, 18 DGE libraries were sequenced using the Illumina HiSeq 4000 platform, generating ∼50 million raw reads for each library. After filtering the reads containing adapter sequences and unknown nucleotides as well as those of low quality, 268.19 million clean reads with 26.82 Gb nucleotides were obtained ( Table 1). De novo assembly yielded 59,173 unigenes with a mean length of 1,017 bp for all 18 libraries. Clean reads were deposited in the NCBI Sequence Read Archive (SRR3741688). Sequence annotation based on seven different public nucleotide/protein databases yielded 46,673 (78.88%) annotated unigenes in at least one of the databases, with ∼21% of unigenes unmapped in any existing database ( Table 2).

Global Gene Expression in the Stem Apex of Primocanes and Floricanes
In total, 56,010 and 55,714 unigenes were expressed at the stem apex of primocanes and floricanes, respectively ( Figure 3A). Of those expressed in primocanes, 41,354 were constitutively expressed at all three developmental stages, whereas 6,983 and  Asterisk represents the number of unigenes which has matched annotation in at least one functional database. 7,673 were specifically expressed in one and two developmental stages, respectively ( Figure 3B). Of the unigenes expressed in floricanes, 42,578 were constitutively expressed at all three developmental stages, whereas 5,309 and 7,827 were specifically expressed in one and two developmental stages, respectively ( Figure 3C). These changes in gene expression indicated that primocane apex rooting of Rubus is an involute biological process regulated by many genes. Compared to SAP2, SAP3, SAF1, SAF2, and SAF3, the greatest abundance of stagespecific transcripts (3,972) was detected in SAP1 of primocanes ( Figures 3B,C), suggesting that a large number of specific genes is needed to coordinate the complex process of stem extension. Compared to SAP1, SAP3, SAF1, SAF2, and SAF3, 930 stage-specific transcripts were detected in SAP2 of primocanes ( Figures 3B,C), indicating that transcripts regulating growth arrest and gravitropic deflection of primocane are expressed at this stage and initiate the growth transition from negative to positive gravitropism.

DEGs in the Stem Apex of Primocanes and Floricanes
FDR ≤ 0.001 and |log 2 Ratio| ≥1 were used as thresholds to identify DEGs (Figure 4; Table 3). Comparisons of primocane and floricanes at the same developmental time points revealed 7,346 (SAP2 vs. SAF2) and 4,110 (SAP3 vs. SAF3) DEGs; this number was much higher than the 498 DEGs in SAP1 vs. SAF1, in which case the number of up-and downregulated genes were similar. In contrast, in SAP2 vs. SAF2, there were more downregulated than upregulated DEGs, suggesting that foliar bud-related genes were silenced as the plant prepared for root differentiation. In SAP3 vs. SAF3, there were more upregulated than downregulated DEGs, possibly because root differentiation-related genes needed to be activated for the regulation of root development.
Similarly, the analysis of floricanes at different time points revealed that most DEGs were downregulated in SAP2 vs. SAP1, while the majority of those in SAP3 vs. SAP2 were upregulated.

Analysis of DEG Clustering and Pathway Enrichment
Hierarchical clustering analysis of differentially expressed transcripts in three SAF vs. SAP paired comparisons revealed that only 131 transcripts were common to all three developmental stages (Figure 5), indicating that stage-specific transcripts were likely responsible for the regulation of AR differentiation. The analysis of DEGs in SAP1 vs. SAP2, SAP2 vs. SAP3, and SAP1 vs. SAP3 paired comparisons revealed 1,365 transcripts common to all three developmental stages in primocane, which grouped into eight categories (Figure 6). Apart from categories D and F, the other six categories (comprising 1,162 DEGs) showed similar trends of up-or down-regulation in SAP1 vs. SAP2 and SAP2 vs. SAP3 paired comparisons, indicating that AR induction and differentiation are closely and sequentially regulated at the level of gene expression.
To obtain functional information about the 1,365 common DEGs, we carried out a literature search and annotated biochemical and biological functions according to the WEGO database. Results of GO functional enrichment are shown in Figure 7. For cellular position, the categories with considerable enrichment and largest number of DEGs were cell (33.3%) and cell part (33.3%). For biological process, the categories with the greatest enrichment and largest number of DEGs were metabolic process (66.7%), cellular process (50.0%), and response to stimulus (14.6%). For molecular function, binding (60.4%) and catalytic activity (54.2%) were the most highly represented categories.
FIGURE 5 | Hierarchical clustering (HCE) analysis of differentially expressed transcripts between stem apex of primocane and floricane at three developmental stages (SAP1-VS-SAF1, SAP2-VS-SAF2, and SAP3-VS-SAF3, "a" was the control and "b" was experimental group in "a-VS-b"). The clusters from A to H indicate the eight major clusters resulting from HCE analysis. Each line refers to data from one gene. The color bar represents the log 10 RPKM and ranges from blue (low expression) to red (high expression).
FIGURE 6 | Hierarchical cluster analysis of DEGs between different developmental stages for stem apex of primocane and floricane (SAP1 vs. SAP2, SAP2 vs. SAP3, SAP1 vs. SAP3, "a" was the control and "b" was experimental group in "a vs. b"). Each line refers to data from one gene. The color bar represents the log 10 RPKM and ranges from green (low expression) to red (high expression). The clusters from A to H indicate the eight major clusters resulting from HCE analysis.
Stage 2 is the transition period from bud development to root formation; as such, DEGs at this stage are presumed to play an important role in determining the developmental course of Rubus primocane apex. We performed a KEGG pathway classification and functional enrichment analysis for DEGs identified in SAP1 vs. SAP2. A total of 12 pathways were significantly enriched, with a Q < 0.05 (Table 4). These DEGs were associated with flavonoid biosynthesis; circadian rhythm-plant; biosynthesis of secondary metabolites; stilbenoid, diarylheptanoid, and gingerol biosynthesis; cyanoamino acid metabolism; sesquiterpenoid and triterpenoid biosynthesis; porphyrin and chlorophyll metabolism; photosynthesis-antenna proteins; limonene and pinene degradation; plant hormone signal transduction; phenylpropanoid biosynthesis; and ascorbate and aldarate metabolism. These results suggest that changes in the expression of genes related to metabolism, environmental adaptation, and signal transduction are involved in the transition from bud development to root formation at the primocane apex of Rubus.

DEGs Associated with Primocane Apex Rooting in Rubus
Based on prior knowledge of their putative involvement in shoot and root development and enriched KEGG pathways of genes that were differentially expressed between primocane and floricane stem apices at three development stages, several subgroups were determined to be related to root induction and  Unigene13079_All, 4.00) in the cysteine and methionine metabolism and plant hormone signal transduction pathways. This latter suggested that ethylene biosynthesis and ethylene and auxin signaling play important roles in regulating shoot elongation arrest and initiation of root differentiation. We also identified a set of shoot and root development-related genes, including SHOOT MERISTEMLESS (STM) (CL6648.Contig1_All, 1.16; CL6648.Contig3_All, 1.84) related to shoot growth, and SHOOT GRAVITROPISM (SGR)5 (Unigene12164_All, −1.00), which is responsible for changes in shoot gravitropism during the transition from bud to root; seven DEGs were involved in regulating root primordium, root phototropism, lateral root primordium (LRP), and root cap, including CL6906.Contig3_All (−1.25), Unigene15061_All (−2.63), Unigene16263_All (2.99), Unigene17114_All (5.13), Unigene20360_All (−1.47), Unigene331_All (−1.49), and Unigene7065_All (4.53). Taken together, their up/downregulation in the SAP1 vs. SAP2 paired comparison was presumed to be beneficial for stem elongation growth arrest, stem gravitropic deflection, and root primordium differentiation.
Stage 3 is the major period in which the root primordium differentiates from stem apex and the root extends. In the SAP2 vs. SAP3 paired comparison, we identified 13 DEGs specific to this period that may facilitate rooting. In the plant hormone signal transduction pathway, two DEGs encoding ARFs (CL1090.Contig5_All and Unigene15505_All) were highly upregulated (5.91 and 4.53 fold, respectively); at the same time, three highly DEGs encoding CTR1 (Unigene24307_All, −3.10; CL5715.Contig2_All, 4.03) and one encoding ETHYLENE-INSENSITIVE (EIN)3 (Unigene1657_All,

DEGs Validated by qRT-PCR
To validate the DEGs identified in the SAP1 vs. SAF1, SAP2 vs. SAF2, and SAP3 vs. SAF3 paired comparisons, six genes were selected for validation by qRT-PCR analysis, including ACS, ARG7 (INDOLE-3-ACETIC ACID-INDUCED PROTEIN ARG7), ERF1, SGR5, RPT3, and RPD1 (ROOT PRIMORDIUM DEFECTIVE 1), all of which have been implicated in the regulation of root development. The expression of all six genes as determined by qRT-PCR was consistent with the DGE patterns (Figures 8A,B).

Ethylene and Auxin Play Essential Roles in the Regulation of Rubus Primocane Apex Rooting
AR formation is a quantitative genetic trait regulated by both environmental (e.g., temperature, light, and relative humidity) and endogenous (e.g., levels of hormone, sugar, mineral salt, and other molecules) factors (Tang et al., 2016;Zerche et al., 2016). Hormones play an important role in this process as plants respond to the changing environment, and influence cell fate specification via regulation of downstream gene expression. This was confirmed by the results of the DEG pathway enrichment analysis, which showed that plant hormone signal transduction (ko04075) and cysteine and methionine metabolism (ko00270) were significantly enriched ( Table 4).
Auxin is the principal phytohormone that initiates rooting and is critical for the first phase of AR formation (Ribeiro et al., 2016). ARF act as an important positive regulator of auxin signaling; two DEGs encoding ARFs were highly upregulated while six encoding ARFs were downregulated in AP2-VS-SAP3 and SAP1-VS-SAP2 paired comparisons, respectively, implying that auxin signaling tended to be activated in stage 3 and inactivated in stage 2. In Arabidopsis, auxin is a positive regulator of ethylene-mediated response in root growth, and ARF17, ARF6, and ARF8 act as negative or positive regulators of adventitious rooting (Rahman et al., 2001;Gutierrez et al., 2009). The upregulation of ARF-encoding DEGs (CL1090.Contig5_All, Unigene15505_All) in the SAP2 vs. SAP3 paired comparison is presumed to be beneficial for AR formation in stage 3, while their downregulation in the SAP1 vs. SAP2 paired comparison is expected to induce stem elongation arrest in stage 2. Aux/IAA acts as a negative regulator of auxin signaling. In accordance with the activation and inactivation of auxin signaling, highly downregulated DEGs (CL130.Contig7_All, −7.70) and upregulated DEGs encoding Aux/IAA (Unigene18450_All, 2.3) were also identified in the SAP1 vs. SAP2 and SAF2 vs. SAP3 paired comparisons, respectively. Thus, auxin may play an important role in regulating stem elongation in stage 2 and root differentiation in stage 3 by inhibiting and inducing the expression of ARFs and AUX/IAA.
Ethylene is thought to interact with auxin in the control of adventitious rooting in stems or stem cuttings (Da et al., 2013). ACS is a key enzyme in ethylene biosynthesis; a DEG (Unigene20575_All) encoding ACS was upregulated by 8.07 and 7.24 fold in the SAP1 vs. SAP2 and SAP1 vs. SAP3 paired comparisons, respectively, implying that ethylene biosynthesis was activated at stages 2 and 3. Wounding and other abiotic stress factors activate ethylene biosynthesis, and AR formation can be stimulated by ERFs (Druege et al., 2014. CTR1 is a negative regulator of the ethylene response pathway in Arabidopsis (Solano et al., 1998). Ethylene signaling was active at stages 2 and 3, as evidenced by the upregulation of positive regulators ERF1 (Unigene13079_All) and EIN3 (Unigene1657_All) and downregulation of CTR1 (Unigene24307_All). Taken together, these results suggest that activation of ethylene biosynthesis via ACS, EIN3, and ERF1 signaling contributes to the termination of bud development at stage 2, and subsequently induces AR formation at stage 3 in Rubus.

The Circadian Clock is an Important System Controlling AR Differentiation in Rubus
Stem elongation in Rubus is enhanced by long-day conditions (Sønsteby and Heide, 2009). The rooting process was previously thought to be a response to the shortening day length in late summer in European Rubi of the subgenus Eubatus (Heslop-Harrison, 1959). The three floricane samples were collected on July 10, August 10, and September 10, and our data suggest that rooting at primocane stem apex is regulated by photoperiod. The significant enrichment of the KEGG circadian rhythm-plant pathway implies that AR formation at the Rubus stem apex involves circadian rhythm-related genes. DEGs in the circadian rhythm-plant pathway included PHYB, APRR5, ELF3, and CDF1, all of which were downregulated in SAF2 vs. SAP2 and SAP1 vs. SAP2 paired comparisons, while HEC2 expression showed the opposite trend. APRR5 is associated with the inhibition of leaf expansion in the Arabidopsis shade-avoidance response (Takase et al., 2013). The depressed expression of HEC showed pin-shaped inflorescences in Arabidopsis thaliana, suggesting that HEC genes are involved in auxin-mediated control of gynoecium patterning (Gremski et al., 2007). In rice, PHYB mutants exhibited reduced total leaf area per plant (Liu et al., 2012). ELF3 is a nuclear protein that is expressed rhythmically and interacts with PHYB to control plant development and flowering, and is an important component of the core circadian clock independent of light conditions. We found that PHYB, APRR5, ELF3, CDF1, and HEC2 were highly expressed at stage 2, corresponding to the shortening of daylight time. PHYB along with ELF3 may regulate circadian rhythm, FIGURE 8 | Validation of DEGs by qRT-PCR analysis. Note: the relative expression levels of six differential expressed genes in three development stages of primocanes and floricanes were obtained by RNA-seq using RPKM method (A) and by qRT-PCR using the 2 − CT method (B). Bars represent mean ± SE (n = 6).
while PHYB and APRR5 may inhibit leaf growth and expansion at the stem apex and induce its final withering. The arrest in bud differentiation may create beneficial physiological conditions for root primordium differentiation. Thus, HEC2 may influence root induction and differentiation via auxin.

Direct Regulator of Primocane Apex Rooting in Raspberry
Roots play important roles in plant growth and development, including water and nutrient uptake and anchoring the plant in soil. In Rubus, AR formation implies the formation of new daughter plants. In Arabidopsis thaliana, RPD1 is associated with the maintenance of active cell proliferation, and its mutation was found to impair axis formation and cellular patterning at early stages of root primordium development (Konishi and Sugiyama, 2006). Tropism induces root growth and determines the eventual vertical growth trajectory of plants. Roots exhibit a negative phototropic growth pattern in response to blue and white light, which is mediated by phytochrome, which may regulate positive phototropism in roots (Kiss et al., 2003). The loss-of-function mutant of RPT2 revealed a role in phot1-induced phototropic response and stomatal opening of roots (Inada et al., 2004). RHD3 may be regulated by auxin and ethylene, and itself regulates cell enlargement and root development (Wang et al., 1997). AIR12 regulates the development of epidermal cells surrounding the emerging lateral root (Costa et al., 2016). RC acting as a gravity-sensing site can sense environmental cues to control the growth direction of the root tip as well as regulate root growth and protect internal cells (Wang et al., 2005;Suzuki et al., 2016). ARFs in auxin signaling are indispensable for root cap development (Wang et al., 2005). In the present study, DEGs encoding RPD1, RC, LRP, RPT2, RPT3, AIR12, and RHD3 were identified in AP1-VS-SAP2 and SAP2-VS-SAP3 paired comparisons, and we concluded that their up-or downregulation facilitate axis formation and cellular patterning, cell enlargement, negative phototropic growth, stomatal opening, and gravitropic sensing and control growth direction of roots, thereby contributing to root differentiation and development in stages 2 and 3. Considering the potential regulatory roles of phytochrome, auxin, and ethylene in root phototropism, we speculate that there is crosstalk between the circadian rhythmplant and plant hormone signal transduction pathways and root development.

Model for Regulation of Rooting Genes at the Stem Apex of Rubus
Accumulating evidence suggests reciprocal interactions between the circadian clock, metabolism, and stress signaling in the control of Arabidopsis growth (Müller et al., 2014). In sunflower (Helianthus annuus L.), gibberellins, auxin, and ethylene directly regulate light-mediated changes in shoot growth. The first two function as growth promoters, whereas ethylene is a growth inhibitor that probably interacts with gibberellins in controlling stem and leaf growth along the sunflower shoot (Kurepin and Pharis, 2014). Normal levels of abscisic acid (ABA) are important for shoot development and suppress ethylene synthesis in tomato (Lenoble and Sharp, 2004). Thus, ABA deficiency during shoot growth may be at least partly attributable to increased ethylene production.
During AR formation at the primocane apex, a series of physiological changes occurred according to the developmental time point (Figure 2), including shoot growth arrest, geotropic growth of the shoot, root primordium differentiation at the apex (stage 2), root development, soil entry of roots, and daughter plant formation (stage 3). Our transcriptome comparisons revealed that a set of DEGs related to circadian rhythmplant, plant hormone signal transduction, and shoot and root development acted coordinately to regulate this unique rooting process. Based on this evidence, we propose a simple model of AR formation at the primocane apex in which photoperiod acts as the environmental inducer (Figure 9). At stage 1, the primocane exhibits a normal extending growth pattern under a long photoperiod. The shorter photoperiod during stage 2 (approach of autumn) induces changes in the expression of genes related to circadian rhythm, including PHYB, ELF3, APRR5, CDF1, and HEC2. These DEGs inactivate auxin signaling via downregulation of ARFs and activate ethylene signaling via upregulation of ACS, EIN2, and ERF1/2. Inactive auxin and active ethylene cause growth arrest and changes in geotropism that may be necessary for subsequent root differentiation, which is induced by the expression of root development-related genes such as RC, LRP, RPD1, and RHD3. At stage 2 and 3, auxin signaling is active while ethylene signaling is inactive, and root differentiation continues and root formation begins under the control of RHD3, RC, RPD1, RPT2, and RPT3. In conclusion, as a plant with a unique rooting habit, not only typical rooting-related DEGs but also components of auxin and ethylene signaling pathways in addition to shoot growth-and circadian rhythm-related genes modulating shoot growth arrest and growth direction were found to be associated with the rooting process.

AUTHOR CONTRIBUTIONS
JL and YC contributed to study conception and design, collection and/or assembly of data, data analysis and interpretation, and manuscript writing. YM, YZ, JX, and YS prepared samples, collected and/or assembled data, and analyzed and interpreted data.