All the plant materials used in this study were collected from the “Hongmanao” pomegranates planting base of Anhui Science and Technology University. Tips of 2–3 cm were retrieved from the tender stems of the pomegranates well grown and free of diseases and insect pests, and cultivated in water for 2–3 days. Then, the growth point of each stem tip was collected and rinsed with clean water in a beaker for 15 min, 75% alcohol for 30 s, and sterile water 2–3 times, to remove surface dust and other impurities. Rinsed materials were then disinfected with 0.1% HgCl₂ for 8 min and rinsed with sterile water 3–5 times. The growth point of the leaf primordium (approximately 0.3–0.5 mm) was collected from the disinfected materials and inoculated to a prepared culture medium. The explants were grown in the WPM medium (containing 0.8 mg/L IBA 3%, sucrose, and 0.55% agar under a pH of 6.5) at the culture chamber under controlled culture conditions: a 12/12 h light/dark photoperiod of 3,000 LX illumination, 25 ± 1°C temperature and 50% relative humidity.

After three cycles of subculture (35–40 days for each cycle), 0.8 mg/L IBA was added to the WPM culture media to induce the proliferation of stem tips. Then, the seedlings of the same growth condition were selected and assigned into four experimental groups, where the first three groups were treated with paclobutrazol, B9, and mannitol, respectively, and the remaining one group was treated with the same amount of distilled water as the control (CK). As shown in Supplementary Table 1, five concentration gradients were set for each drug, and each treatment/concentration included five plants as biological replicate. After 30, 60, and 90 days of treatment, we measured the height of the above-ground parts of each plant in each treatment. According to their growth conditions, plants treated at two desired concentrations were selected for each PGR (6 and 8 mg/L for paclobutrazol, 6 and 8 mg/L for B9, and 2.5 and 15 g/L for mannitol) for downstream experiments. Leaves were collected from each selected plant with the treatment of 90 days and activities of antioxidant enzymes, including superoxidase dismutase (SOD) and peroxidase (POD), and the reduced glutathione/oxidized glutathione (GSH/GSSG) ratio were measured following the methods described in Qian et al. (2020) to assess the physiological status of plants under treatment. Meanwhile, fresh leaves were also collected from two plants (biological replicates) of each treatment of 90 days (CK and PGRs), frozen in liquid nitrogen and stored at −80°C for RNA sequencing (RNA-seq).

Total RNA extraction, RNA-seq, and data preprocessing were carried out as described in Qian et al. (2020). High-quality reads were aligned to the P. granatum reference genome (Luo et al., 2020) using hisat2 program (Kim et al., 2019). The reads mapped to more than one genome location were filtered out and only those uniquely mapped were retained for downstream analyses. For each gene, the expression level was assessed by the abundance of mapped reads using the featureCounts program from the Subread package (Liao et al., 2014). DESeq2 (Love et al., 2014) was employed for differential expression analysis between each of the PGR treatments (paclobutrazol, B9, and mannitol of different concentrations) and CK. Genes whose |log₂ fold-change _treatment/CK| > 1 and adjusted p < 0.05 were considered as significantly differentially expressed (denoted as “DEGs”). Hereafter, DEGs are denoted as “upregulated” or “downregulated” if their expression after treatment was higher or lower than CK. Raw read counts of genes were normalized by the fragments per kilobase of exon per million fragments mapped (FPKM), to correct between-sample differences in sequencing depth.

Genome-wide functional annotations, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations, were adopted from our previous study (Qian et al., 2020). For each of the comparisons between PGR-treated plants and CK, GO and KEGG enrichment analyses were implemented for up- and downregulated DEGs, respectively. The significance levels were assessed by Fisher’s exact test, where GO terms and KEGG pathways were considered as significantly overrepresented if their p < 0.05.

To access the dose-effect of each PGR (paclobutrazol, B9, and mannitol), trend analysis was carried out across different treatment concentrations using Mfuzz package (Futschik and Carlisle, 2005; Kumar and Futschik, 2007). Briefly, for each treatment, highly variable genes were grouped in six clusters via mfuzz function, according to the transcriptional similarities across different treatments. Then, the significantly enriched clusters were identified by Fisher’s exact test as described in Qian et al. (2020).

Compared to untreated conditions (CK), all the PGRs were found to significantly suppress the growth of pomegranate seedlings (Figures 1A,B and Supplementary Figure 1). After 30 days of treatments, plants treated with paclobutrazol showed the most reduction in height (reduced by 74.0–79.4% compared to CK; 1.2–1.6 cm in average for five concentrations vs 6.0 cm for CK; p < 0.01 for all the five concentrations), followed by mannitol and B9 (Supplementary Figure 1A). When extended to 60 days, the height of the paclobutrazol-treated plants was only 13.2–18.0% of that of CK (p < 0.01), and the plant height under the treatment of mannitol was slightly higher than that under B9 treatment (Supplementary Figure 1B). After 90 days, the height ranged from 2.1 to 4.1 cm with the paclobutrazol treatment from 5 to 9 mg/L, which were 64.2–81.6% shorter than that of CK (11.3 cm in average) (p < 0.01; Figure 1A). Comparatively, the average height of B9- and mannitol-treated seedlings was 4.7–7.9 and 6.4–8.2 cm, respectively, which were also significantly shorter than CK (p < 0.01; Figure 1A).

Under the treatments of both paclobutrazol and B9, the activities of SOD and POD, as well as the GSH/GSSG ratio, in seedlings were found to be substantially enhanced for each of the PGR treatments than CK (p < 0.05 or <0.01 for all comparisons between treated plants and CK; Figures 1C,D and Supplementary Figure 2). For both paclobutrazol and B9, the enzymes were more activated as the increasing of treatment concentrations that SOD and POD activity increased by 169.9–301.0 and 54.1–164.8% over CK, respectively. In contrast, the activities were lower with the mannitol treatment of high concentration (15 g/L) than the low concentration (2.5 g/L). It might be due to the poor growth condition when treated with over-dose mannitol. These results suggested that, besides dwarfing, these exogenous PGR supplies may also lead to adverse impacts on plant normal physiological status.

In a total, 23.87–30.76 million raw reads were obtained from RNA-seq, where 23.81–30.61 million reads passed the quality control. After alignment, 747 and 3,917 DEGs were detected under 6 and 8 mg/L treatments of paclobutrazol, respectively (Supplementary Figures 3A,B), and 3,132 and 2,793 genes were differentially expressed when treated with 6 and 8 mg/L B9 (Supplementary Figures 3C,D). We further identified 1,132 and 815 DEGs under the treatments of 2.5 and 15 g/L mannitol (Supplementary Figures 3E,F). The results of qPCR validation showed a high correlation between the gene expression levels from RNA-seq and qPCR (Pearson correlation coefficient = 0.83; p = 6.29e-7; Figure 2), indicating a reliable inference of gene expression by our RNA-seq data.

For paclobutrazol, after the treatment of 6 mg/L, there were 394 up- and 353 downregulated DEGs, respectively (Supplementary Figure 3A). Comparatively, 1,844 and 2,073 genes were over- and underexpressed when treated with 8 mg/L (Supplementary Figure 3B). For downregulated DEGs, 176 were shared between the two treatment concentrations. GO enrichment analysis showed that these commonly suppressed genes were mainly involved in shoot apical meristem development and auxin biosynthesis and metabolism (Figure 3A). Comparatively, the genes specifically inhibited by 6 mg/L treatment were enriched for the cell proliferation process, while those specific to 8 mg/L treatment were relevant to developmental growth, actin filament bundle assembly, and shoot system development (Figure 3A). We then investigated the expression levels of three key genes along the tryptophan-dependent IAA biosynthetic pathway, including those encoding tryptophan aminotransferase-related protein (TAR), L-tryptophan–pyruvate aminotransferase 1 (TAA1), and indole-3-pyruvate monooxygenase (YUCCA) (Mashiguchi et al., 2011; Mano and Nemoto, 2012). The results showed that, in general, the expression of these genes was suppressed by paclobutrazol treatments, at both concentrations (Figure 3B). In particular, most of the TAA1 and YUCCA genes were downregulated under the treatment of 8 mg/L paclobutrazol. Comparatively, there was no consistent pattern across different members of indole-3-acetic acid-amido synthetase (GH3) genes, which respond to an excess amount of auxin in plants (Mashiguchi et al., 2011).

Trend analysis showed different impacts across the treatments of different paclobutrazol doses. One of the two sufficiently overrepresented trends was the genes were of similar level of downregulation under the treatment of either 6 or 8 mg/L, which were mainly involved in cell cycle, DNA replication proteins and shoot apical meristem development (Supplementary Figure 4A). For the second scenario, the genes were of more reduced expression as the increase of paclobutrazol concentration, and these genes were found to participate in growth hormone synthesis and primary shoot apical meristem specification (Supplementary Figure 4B), indicating the plausible mechanism underlying the severer growth suppression induced by the increase of paclobutrazol dosage (from 6 to 8 mg/L). Meanwhile, for the treatments of both concentrations, genes responding to oxidative stress were observed to be upregulated (Supplementary Figures 4C,D), which was corresponding to the enhanced activities of antioxidant enzymes (Figures 1B,C). It may reflect a stress response in paclobutrazol seedlings under the treatment of paclobutrazol.

Concerning B9, 1,726 and 1,410 genes were significantly downregulated after 6 and 8 mg/L treatment (Supplementary Figures 3C,D, 5A). For both two treatments, genes involved in BR metabolism, chromocenter, and actin filament bundle assembly were suppressed compared to CK (Supplementary Figure 5B). Each treatment also has specific impacts on pomegranates: auxin biosynthetic regulation and cell wall modification during cell growth were downregulated under 6 mg/L B9 treatment (Supplementary Figure 5C), while cell division, developmental growth, and ABA degradation processes were exclusively depressed when treated with 8 mg/L B9 (Supplementary Figure 5D). Comparatively, genes related to oxidative stress resistance were significantly induced by the treatments (Supplementary Figures 5E–H). Trend analysis showed that plant growth and primary shoot apical meristem specification were more suppressed with the increase of concentration (from 6 to 8 mg/L), while the antioxidative processes, such as hydrogen peroxide response and peroxisome and glutathione metabolism, became stronger accordingly (Supplementary Figure 6).

We also assessed the different transcriptional responses between paclobutrazol- and B9-treated plants. In total, 62 genes were downregulated under all four treatments, and they mainly participated in regulating IAA synthetase and cell cycle, including the processes of DNA replication and chromocenter (Figure 4). Moreover, 395 genes were commonly suppressed when treated with 8 mg/L paclobutrazol and 6 and 8 mg/L B9, which were enriched for shoot apical meristem specification, cell cycle regulation, and plant hormone signal transduction (Figure 4). It is noteworthy that, besides these commonly depressed pathways, the processes involved in strigolactone response, shoot system development, and auxin metabolism were inhibited exclusively by 8 mg/L paclobutrazol, while carbohydrate-derivative metabolic process was suppressed specifically by 8 mg/L B9 treatment (Figure 4), suggesting different regulatory mechanisms underlying the dwarfing by different PGRs.

When treated with mannitol, 451 up- and 681 downregulated genes were induced at 2.5 g/L (Supplementary Figure 3E), while the expression level of 381 and 434 genes was increased or reduced at 15 g/L (Supplementary Figure 3F). For both treatments, the auxin biosynthetic process was significantly suppressed (Figures 5A,B). Corresponding to our expectation, the expression of genes related to cell division and plant developmental growth was even reduced when treated with 15 g/L compared to 2.5 g/L. In contrast, cell proliferation and DNA replication processes were more inhibited under the treatment of 2.5 g/L than 15 g/L (Figure 5C). It is noteworthy that, the pathways involved in the response to oxidative stress, such as the removal of reactive oxygen species (ROS) and the regulation of osmotic equilibrium, were largely enhanced (Supplementary Figure 7), suggesting the undesired side effect of mannitol on plants. Furthermore, this hazardous impact got even worse with the increase of treatment concentrations. Compared to 2.5 g/L treatment, genes of the negative regulation of cell death showed reduced expressions when treated with 15 g/L mannitol (Supplementary Figure 7).