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.

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 above-mentioned samples. Total RNA was extracted from the dissected tissue using the RNA Easyspin Isolation System (Aidlab Biotech, Beijing, China). RNA was dissolved in diethylpyrocarbonate-treated 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).

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 BLAST-predicted 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 samples—the false discovery rate (FDR)—was used to determine the threshold P-value in multiple tests and analyses. FDR ≤ 0.001 and |log₂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).

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).

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).

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 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 stage-specific 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.

FDR ≤ 0.001 and |log₂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.

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.

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.

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 differentiation (Table 5) in SAP1-VS-SAP2 and SAP2-VS-SAP3 paired comparisons of Rubus. In the former, we identified 24 unigenes involved in the regulation of root induction and differentiation at the primocanes apex, whose expression changed significantly; these included PHYTOCHROME (PHY)B (CL1524.Contig3_All, −1.21), EARLY FLOWERING (ELF)3 (CL2434.Contig1_All, −1.71), ARABIDOPSIS PSEUDO-RESPONSE REGULATOR (APRR)5 (CL5214.Contig6_All, −10.36), CYCLING DOF FACTOR (CDF)1 (Unigene16183_All, −3.45), and HECATE (HEC)2 (Unigene21192_All, 3.22) in the circadian rhythm–plant pathway; and 1-AMINOCYCLOPROPANE-1-CARBOXYLATE SYNTHASE (ACS) (Unigene20575_All, 8.07), AUXIN RESPONSE FACTOR (ARF) (CL1090.Contig5_All, −9.79; CL2444.Contig1_All, −7.96; CL1556.Contig5_All, −7.51; Unigene20581_All, −3.30; CL5079.Contig1_All, −3.16; CL1556.Contig2_All, −3.17; Unigene15505_All, 4.58.), THREONINE-PROTEIN KINASE CTR1 (CTR1) (Unigene24307_All, 7.91), and ETHYLENE-RESPONSIVE TRANSCRIPTION FACTOR (ERF)1 (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, 7.92) were identified. Auxin and ethylene signaling pathways were more active in stage 3 than in stage 2. A set of root development-related DEGs were identified, including three DEGs encoding ROOT CAP (RC) (Unigene17114_All, 1.53; Unigene7065_All, 3.10; Unigene811_All, −2.07), three encoding ROOT PHOTOTROPISM (RPT)2 (CL1554.Contig3_All, −1.65; CL1554.Contig1_All, −1.60; CL1554.Contig2_All, −1.26), and one each encoding ROOT HAIR DEFECTIVE (RHD)3 (CL4872.Contig3_All, 3.27) and AUXIN-INDUCED IN ROOT CULTURES (AIR)12 (Unigene10761_All, 1.22); their up/downregulation in SAP2 vs. SAP3 paired comparisons was presumed to facilitate to root cap and negative root phototropism formation.

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).

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, 2016). 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.

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, 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.

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 rhythm–plant and plant hormone signal transduction pathways and root development.

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 rhythm—plant, 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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.