Transcriptome analysis of substrate temperature effects on adventitious root formation in peach rootstocks

Plant material and sample preparation

This experiment was conducted under controlled conditions in a greenhouse at the Gansu Academy of Agricultural Sciences to investigate the hardwood cutting propagation of the peach GF677 rootstocks. The greenhouse temperature is 18/5 °C (day/night), the photoperiod is nine hours of light/15 h of darkness, and the humidity is 55–75%. We will take the rootstock branches in the middle of November 2024. The cuttings were selected from the middle and upper portions of current-year branches, ensuring they were free of pests and diseases and had diameters ranging from 0.5 to 1.0 cm, and after being stored in sand for one month, hardwood cuttings will be taken for propagation. The optimal length of the cuttings was 15 to 20 cm. A cutter was used to create a 40 to 45° bevel at the base near the basal bud. The top of the cutting was cut flat, and it was advisable to apply a plant wound healing agent to the flat surface. The substrate consisted of a 1:1:1 volume ratio of perlite, peat moss, and vermiculite. The cuttings were dipped in a solution of 200 mg L⁻¹ indole-3-butyric acid (IBA) for 20 s before planting, while the control group was dipped in water (El-Boray et al., 1995; Eliwa & Wahba, 2018; Zhang & Wang, 2018; Zhang et al., 2013). They were heated using electric heating elements (ambient temperature (CK), 19 °C, 22 °C, 25 °C, and 28 °C) and the rooting rates of the five treatments were counted after 40 days. Sampling was conducted on the 7th day (labeled as CT1, CT2, CT3, CT4, and CT5), the 14th day (labeled as, CT6, CT7, CT8, CT9 and CT10), and the 21st day (labeled as, CT11, CT12, CT13, CT13 and CT14) for RNA sequencing. CT0 indicates that the cuttings were treated with IBA and were planted after 0 days at ambient temperature. CT1, CT6 and CT11 represent cuttings treated with IBA and planted after 7, 14 and 21 days respectively at ambient temperature. For detailed information on the experimental design of this study, please refer to Table 1. In this study, three biological replicates were set for each treatment to ensure the reliability of the results. Each replicate contained 30 randomly selected healthy cuttings, which were taken from multiple plants in the greenhouse to reduce genetic variability.

RNA library construction and high-throughput sequencing

We sampled the phloem from two cm of the base of each cutting, mixed them evenly and extracted RNA. Following the manufacturer’s instructions, total RNA was extracted from the plants using a Pure Plant RNA Extraction Kit (Tiangen, China). The concentration and purity of the extracted RNA were measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The Hieff NGS Ultima Dual-mode mRNA Library Prep Kit (compatible with Illumina platform, Yeasen, China) was utilized to construct the sequencing library. The library fragments were purified using AMPure XP beads (Beckman Coulter, USA). Subsequently, the cDNA products were amplified by polymerase chain reaction (PCR), and finally sequenced on an Illumina HiSeq2500 genomic sequencer by Biomarker Biotechnology Co., Ltd (Song et al., 2013).

RNA-sequencing data analysis

Three biological replicates were set for each sample during RNA-sequencing. The statistical power of this experimental design, calculated in RNASeqPower is 0.84. By processing the raw data, we eliminated reads containing adapters, poly-N sequences, and low-quality reads to obtain high-quality clean reads. All subsequent downstream analyses were conducted based on these high-quality clean reads. Subsequently, we aligned these clean reads to the reference genome sequence (Prunus persica) and annotated them accordingly. To achieve rapid and accurate alignment of clean reads against the reference genome and to obtain their precise genomic coordinates, we utilized the HISAT2 v2.0.5 software (default parameters) (Mortazavi et al., 2008). The HISAT2 index was constructed from the reference genome Fasta file using default settings (hisat2-build). This configuration uses dynamic scoring for mismatch tolerance without explicit threshold specification. Furthermore, the SAMtools analysis tool is used to analyze the alignment status of specific regions (such as intergenic regions, exon regions, intron regions). Subsequently, transcript assembly was conducted using StringTie v1.3.6 (Pertea et al., 2015) in three stages: Per-sample assembly: stringtie <sample.bam>-p 4 -G <reference.gtf>-o <sample.gtf>; Merge assemblies across samples: stringtie –merge -G <reference.gtf>-l novel -o merged.gtf sample1.gtf sample2.gtf ...; Novel transcript identification: gffcompare -R -r <reference.gtf>-o gffcompare merged.gtf. Transcripts with no overlap to reference annotations (class code “u”) were retained as novel genes. The quantitative analysis was conducted using the featureCounts v1.5.0-p3 tool within the subread software (Liao, Smyth & Shi, 2014). Set the minimum mapping quality score to 10. Finally, gene expression levels were estimated based on fragments per kilobase of transcript per million mapped reads (FPKM) (Trapnell et al., 2010). The DESeq2 R package (v 1.20.0) was adopted to conduct a significant analysis of the gene expression differences of the rootstocks of peach under different treatment conditions (Anders & Huber, 2010). P-values were adjusted for the false discovery rate (FDR) using the Benjamini–Hochberg method. The genes that were identified by DESeq2 with an adjusted P-value ≤ 0.05 and |log₂FC| ≥ 1 were regarded as DEGs. The Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs sets were conducted using the clusterProfiler software (Yu et al., 2012).

Weighted gene co-expression network construction

The gene co-expression network was constructed using the WGCNA R package. The Pearson correlation coefficient is computed based on gene expression profiles across diverse samples, followed by transformation into an adjacency matrix using a weighting function (Wang et al., 2023). The scale-free topology criterion was used to determine the soft threshold power (β), thereby ensuring that the gene expression matrix met the requirements of a scale-free network. A minimal module size of 30, a merge cut height of 0.25, and a area assign threshold of 0. Subsequently, visualization of gene modules is performed using Cytoscape.

Root phenotype analysis

We conducted a 40-day phenotypic analysis of the root systems of peach rootstock cuttings propagated under different substrate temperatures of 19 °C, 22 °C, 25 °C and 28 °C. The results revealed that root development was optimal under the condition of 25 °C (Fig. 1A). We further conducted statistics and analysis on the rooting rate of the rootstocks at 40 days post-cutting (Fig. 1B). The results indicated that the highest rooting rate, 91%, was obtained at a substrate temperature of 25 °C. Consequently, we can conclude that the optimal substrate temperature for the propagation of peach rootstocks by cutting is 25 °C.

Overview of RNA-Seq data

We conducted an in-depth sequencing analysis on 45 libraries of peach rootstocks under three different cutting periods (7, 14, and 21 days) and five substrate temperature conditions (ambient temperature (CK), 19 °C, 22 °C, 25 °C, and 28 °C). This sequencing effort yielded over 2.1 billion raw reads, with an average of approximately 45 million reads per library (Table 2). After rigorous filtering process, RNA-Seq technology generated between 39 million and 48 million high-quality clean reads for each sample, with an average Q20 value exceeding 98.62% and a Q30 value exceeding 95.79% (Table 2). These high-quality data were utilized for all subsequent expression analyses. The filtered and cleaned reads were successfully aligned to the reference genome of Prunus persica, with an average genome alignment rate of 91.02% (Table 2). The multiple mapped category accounted for the largest proportion of 1.92% (CT0, CT1), while the unique mapped category had the lowest proportion of 82.1% (CT14). Furthermore, we conducted detailed statistical analysis on the alignment regions of each sample and found that exon regions accounted for up to 96.26% of the total number of reads (File S1).

We performed a comprehensive analysis of gene expression levels within our dataset using the featureCounts tool from the subread package. Specifically, File S2 presents a detailed read count expression matrix for 26,978 individual genes across various samples. Notably, the three biological replicates of the samples exhibited high Pearson correlation coefficients (R² ≥ 0.773, Fig. 2A), indicating strong reproducibility in our data. To gain further insights into the differences and relationships between the samples, we employed principal component analysis (PCA) separated samples into 16 distinct clusters corresponding to the 15 experimental treatments (5 substrate temperatures ×3 cutting periods) plus the untreated control (CT0). Among them, peach rootstock cuttings incubated for 7 days at different substrate temperatures (ambient temperature (CK), 19 °C, 22 °C, 25 °C, 28 °C) are classified as CT1, CT2, CT3, CT4, and CT5, respectively; those incubated for 14 days at different substrate temperatures are classified as CT6, CT7, CT8, CT9, and CT10; while another set of cuttings incubated for 21 days are classified as CT11, CT12, CT13, CT14, and CT15. The control group is classified as CT0 (Fig. 2B and File S3). Furthermore, we employed various visualization techniques, including violin plots and boxplots to conduct a comprehensive and meticulous examination of the gene expression distribution across the samples (Figs. 2C–2D).

Screening and analysis of DEGs

We conducted a detailed comparative analysis on 12 groups of DEGs (P-value ≤ 0.05, —log₂FC—≥ 1). We compared CT1 with CT2, CT3, CT4, and CT5; CT6 with CT7, CT8, CT9, and CT10; CT11 with CT12, CT13, CT14, and CT15 respectively. Overall, the number of DEGs increased with the rise of substrate temperature at 7 days and 21 days of cutting. However, during the 14-day cutting treatment period, the increase in DEGs continued until 25 °C (Fig. 3A). Specifically, the highest count of DEGs in peach rootstock cuttings was observed on the 7th and 21st days after cutting at substrate temperatures of 25 °C and 28 °C. Conversely, the peak number on the 14th day post-cutting was noted at substrate temperatures ranging from 22 °C to 25 °C (Fig. 3A). From the significance analysis of upregulated and downregulated genes, the optimal substrate temperatures range for peach rootstock cuttings at 7- and 21-days post-insertion was determined to be 25 °C to 28 °C, whereas for cuttings at 14 days post-insertion, the optimal range shifted to 22 °C to 25 °C. Notably, the samples collected on the 21st day after cutting showed more prominent variations in differential gene expression levels under different substrate temperatures (Fig. 3A). Based on the statistical analysis of rooting rates for cuttings at various substrate temperatures (Fig. 1) and the statistical analysis of the number of DEGs, we conclude that the optimal substrate temperature for peach rootstock cuttings is 25 °C. Subsequently, we performed a heatmap analysis on the 14,380 DEGs identified in various comparisons (Fig. 3B).

More notably, there were 231 DEGs shared among the four comparisons, while 1,814, 926, 720, and 1,312 DEGs were uniquely expressed in the comparisons of CT1 vs CT0, CT4 vs CT1, CT9 vs CT6, and CT14 vs CT11, respectively. Among them, the meanings represented by each treatment are: CT0: control at 0 d ambient; CT1: control at 7 d ambient; CT4: 7 d at 25 °C; CT6: control at 14 d ambient; CT9: 14 d at 25 °C; CT11: control at 21 d ambient; CT14: 21 d at 25 °C. Furthermore, in the comparison between CT1 and CT0, there were 2,629 upregulated genes and 1,512 downregulated genes. In the comparison between CT4 and CT1, there were 1,678 upregulated genes and 1,709 downregulated genes. Similarly, in the comparison between CT9 and CT6, there were 1,635 and 1,109 single genes were upregulated and downregulated, respectively. In the comparison between CT14 and CT11, there were 571 single genes were upregulated, with 1,112 single genes downregulated (Fig. 3C).

GO enrichment and KEGG pathway analysis of DEGs

Based on the above conclusions, we will focus on conducting in-depth research on cuttings that have been propagated for 7, 14, and 21 days at a substrate temperature of 25 °C. These cuttings will be compared with the control group, and subjected to GO and KEGG pathway analysis. When comparing CT1 vs CT4, the most significantly enriched GO terms in biological processes were “microtubule-based movement”, “cellular process”, and “photosynthesis”. Among the cellular components, three GO terms “thylakoid part”, “photosynthetic membrane”, and “photosystem” were highly enriched. In the comparisons of CT6 vs CT9, the most significantly enriched GO terms in biological processes were “response to stress” and “defense response”. When comparing CT11 vs CT14, the most significantly enriched GO terms in biological processes were “movement of cell or subcellular component”, “microtubule-based process”, and “DNA-dependent DNA replication” (Fig. 4A). Among them, the meanings represented by each treatment are: CT1: control at 7 d ambient; CT4: 7 d at 25 °C; CT6: control at 14 d ambient; CT9: 14 d at 25 °C; CT11: control at 21 d ambient; CT14: 21 d at 25 °C. Across these comparisons, “copper ion binding” and “microtubule motor activity” were significantly enriched in terms of molecular function (Fig. 4A).

KEGG enrichment analysis was also conducted to identify the enriched metabolic pathways. The results indicated that significantly enriched KEGG pathways were mainly associated with photosynthesis processes, specifically including “Starch and sucrose metabolism” (Fig. 4B). Furthermore, several other enriched pathways related to hormones and metabolites were identified, including “Flavonoid biosynthesis”, “Plant hormone signal transduction”, and “Zeatin biosynthesis” (Fig. 4B).

Identification DEGs involved in the formation of AR-related pathways

We systematically screened for DEGs that may be involved in AR formation under various cutting substrate temperatures and cutting periods. Previous studies have shown that ARFs regulate the development of adventitious roots by specifically binding to auxin response elements (Tao et al., 2023). In this study, we successfully identified five key ARFs (gene id: 18775811, 18774376, 18772269, 18791569, and 18786733). During the process of adventitious root formation, these genes increased by 1.33 to 3.46 times compared with the control group. Notably, the LATERAL ORGAN BOUNDARIES DOMAIN (LBD) genes, as the primary downstream targets of ARFs, have been shown to be involved in the process of adventitious root formation (Kirolinko et al., 2024). In this screening, we found six differentially expressed LBD genes (gene id: 18774258, 18767664, 18777846, 18789146, 18786806, and 18774612) with varying expression patterns across different samples. Notably, the expression levels of LBDs in the cutting group were 1.36 to 6.85 times higher than those in the control group. The SAUR genes constitute the largest group of specific auxin-responsive factors participating in root development (Zhou et al., 2024). Among the screened 10 SAURs (gene id: 18770094, 18787544, 18793348, 18783824, 18789144, 18768105, 18783913, 18,784,823, 18792443, and 18785952), their expression levels exhibited significant differences compared to control samples, with the highest fold change being 42.09. On the other hand, members of the auxin-responsive GH3 family play a crucial role in regulating auxin homeostasis through the synthesis of auxin conjugates in higher plants (Caño et al., 2018). Notably, the GH3 gene (gene id: 18768891) showed prominent differential expression across various sample tests (Fig. 5). In summary, these genes closely associated with auxin signaling may collectively participate in the formation of AR.

In A. thaliana, genes related to root formation are influenced by auxin, such as LRP1, RGF, and AIR12 (Singh et al., 2020; Shinohara, 2021; Gibson & Todd, 2015). During the period of adventitious root formation, the expression level of the root-related gene AUR3 (gene ID: 18786669) was significantly increased compared to the control group, with an increase of up to 3.8 times. Concurrently, during the AR formation period, the maximum expression levels of LRP1 (gene ID: 18786210, and 18781220) increased by 18.2-fold and eight-fold, respectively. Additionally, four RGF1 genes (gene ID: 18769572, 18790800, 18777415, and 18775914) exhibit increased expression during AR formation. Furthermore, the rooting-related gene AIR9 (gene ID: 18780861) was also significantly upregulated during the formation of AR. Its expression reached the maximum value 7 days after cutting at 28 °C, increasing by 2.58 times compared to the control group (Fig. 5). These findings suggest that cellular division and proliferation activities are significantly enhanced during the development of ARs in peach cuttings.

Furthermore, the NAM/ATAF1/2/CUC2 (NAC) transcription factor and the APETALA2 (AP2) genes are also widely involved in multiple aspects of plant organ development, such as root stem cell development and cell differentiation (Li et al., 2019; Guyomarc’h, Boutté & Laplaze, 2021). Through in-depth analysis, we have identified eight differentially expressed NAC genes (gene ID: 18784579, 18787178, 18768537, 18770366, 18790018, 18772621, 18789474, and 18790535) and six differentially expressed members of the AP2 genes (gene ID: 18777768, 18773161, 9978767, 18773090, 18766891, and 18773658). During the process of adventitious root formation, the maximum fold changes of NAC and AP2 genes compared with the control group were 4.61-fold and 7.79-fold, respectively (Fig. 5). The changes in the expression of these genes indicate that their regulatory roles in the formation of adventitious roots.

Impact of substrate temperatures on adventitious root formation in peach hardwood cuttings via WGCNA

In order to precisely screen out the genes closely related to the response of cutting substrate temperature and period, we first evaluated the number of candidate genes. If the total number of genes did not exceed 45,000, we used the expression data of all genes as the basis for the subsequent weighted gene co-expression network analysis. After this screening procedure, we finally determined 26,978 genes, which became the core for constructing the weighted gene co-expression network. We selected the power value corresponding to R² of 0.8 as the soft threshold (File S4).

Based on the correlation of expression levels between genes, the construction of clustering trees, and the stability of modules, WGCNA analysis precisely divided all genes into 68 modules (Fig. 6A), with the number of genes in each module ranging from 34 to 4,072. Among these, 44 modules exhibited a high correlation with the samples mentioned in this study (R > 0.90, Fig. 6B). Notably, we found an extremely close correlation between the CK1 (CT0) sample and the MEmediumorchid module (R = 0.994572), and a tight correlation between the CK3 (CT0) sample and the MEdarkolivegreen module (R = 0.996928). Additionally, significant correlations were also observed between the T6-1 sample and the MEpaleturquoise module (R = 0.994412), as well as between the T9-2 sample and the MEplum module (R = 0.994312). Based on the clear heatmap of FPKM values, four modules exhibited distinct expression patterns at specific time points (Fig. 6C). Further analysis revealed that these four modules contained various transcription factors involved in the formation of AR-related genes, including four members of the WRKY family, eight members of the ERF family, six members of the NAC family, 11 members of the bHLH family, five members of the bZIP family, and 14 members of the MYB family (File S5). These transcription factors play crucial regulatory roles in the formation of AR, and their high expression levels during the corresponding periods fully support this point.