Genome-wide alternative polyadenylation dynamics underlying plant growth retardant-induced dwarfing of pomegranate

Dwarfed stature is a desired agronomic trait for pomegranate (Punica granatum L.), with its advantages such as lower cost and increased yield. A comprehensive understanding of regulatory mechanisms underlying the growth repression would provide a genetic foundation to molecular-assisted dwarfing cultivation of pomegranate. Our previous study induced dwarfed pomegranate seedlings via exogenous application of plant growth retardants (PGRs) and highlighted the important roles of differential expression of plant growth-related genes in eliciting the dwarfed phenotype of pomegranate. Alternative polyadenylation (APA) is an important post-transcriptional mechanism and has been demonstrated to act as a key regulator in plant growth and development. However, no attention has been paid to the role of APA in PGR-induced dwarfing in pomegranate. In this study, we characterized and compared APA-mediated regulation events underlying PGR-induced treatments and normal growth condition. Genome-wide alterations in the usage of poly(A) sites were elicited by PGR treatments, and these changes were involved in modulating the growth and development of pomegranate seedlings. Importantly, ample specificities were observed in APA dynamics among the different PGR treatments, which mirrors their distinct nature. Despite the asynchrony between APA events and differential gene expression, APA was found to regulate transcriptome via influencing microRNA (miRNA)-mediated mRNA cleavage or translation inhibition. A global preference for lengthening of 3’ untranslated regions (3’ UTRs) was observed under PGR treatments, which was likely to host more miRNA target sites in 3’ UTRs and thus suppress the expression of the corresponding genes, especially those associated with developmental growth, lateral root branching, and maintenance of shoot apical meristem. Together, these results highlighted the key role of APA-mediated regulations in fine-tuning the PGR-induced dwarfed stature of pomegranate, which provides new insights into the genetic basis underlying the growth and development of pomegranate.

Introduction

Pomegranate (Punica granatum L.) is one type of the economic fruit trees that are widely cultivated across the globe. Because it is rich in vitamins and has antioxidant and anti-inflammatory properties in fruits, the health benefits of pomegranate are highly regarded, such as preventing or alleviating diseases and lowering high blood pressure or high cholesterol levels (National Center for Complementary and Integrative Health, NCCIH; Bourekoua et al., 2018; Shahamirian et al., 2019; Asrey et al., 2020; Turrini et al., 2020). With the fast-rising demand for pomegranate products, more and more attention has been paid to screen and breed pomegranate cultivars with the desired high fruit yield and quality. Dwarfing cultivation is one of the major focuses because of its advantages in plant photosynthetic efficiency, fruit production, and disease resistance compared to normal growing mode (Seleznyova et al., 2008; Foster et al., 2017; Wang et al., 2018; Zhou and Underhill, 2021). Qian et al. (2022) demonstrated that exogenous application of plant growth retardants (PGRs) can successfully elicit dwarfed pomegranate seedlings. Comparative transcriptome analysis further unraveled that PGR-mediated downregulation of plant growth hormone synthesis played a central role in inducing the dwarfed stature of pomegranate, providing new clues for molecular breeding of favorable dwarfed pomegranate varieties. Besides gene transcription, plant transcriptome is also under the regulation of post-transcriptional mechanisms, which have been demonstrated as a key contributor to the phenotypic plasticity of plants (Ye et al., 2019; Zhou et al., 2019; Singh and Roychoudhury, 2021). However, our current knowledge on the functional importance of post-transcriptional processes in pomegranate is still limited.

Polyadenylation [poly(A)] is an important post-transcriptional mechanism in eukaryotes that modulates mRNA maturation from the precursor mRNA (pre-mRNA). It includes two coupled steps: endonucleolytic cleavage at the 3’ end of pre-mRNA and the addition of a poly(A) tail at the cleavage sites (Colgan and Manley, 1997; Tian and Manley, 2017). More importantly, for many genes, the cleavage and poly(A) signal recognition occur at multiple positions, that is, giving rise to multiple mRNA isoforms with different lengths, which is referred to as alternative polyadenylation (APA). APA events have been demonstrated to be widespread across genomes; for example, over 70% of the Arabidopsis genes were found to possess more than one poly(A) site (Wu et al., 2011; Elkon et al., 2013). These APA events may alter the stability and translation of mRNA or the length of the resulting protein products; thus, APA serves as a key contributor to the complexity of eukaryotic transcriptome (Shen et al., 2008; Di Giammartino et al., 2011; Sun et al., 2012; Tian and Manley, 2017). Recent studies have highlighted the biological importance of APA in regulating plant growth, development, and resistance to environmental stresses (de Lorenzo et al., 2017; Zhou et al., 2019; Yu et al., 2022; Wang et al., 2023). For instance, Yu et al. (2022) performed a genome-wide investigation to APA dynamics underlying Arabidopsis leaf ontogeny and showed that the largest changes in poly(A) site usage occurred at the early stage of true leaf development, while the APA levels experienced a reduction along the developmental process. Furthermore, it was shown that these APA genes participated in modulating the biological processes associated with leaf development, for example, response to phytohormone. These findings highlighted the essential roles of APA-mediated post-transcriptional regulations in plant growth and development. However, the APA mechanisms underlying PGR-induced dwarfing have not been investigated in pomegranate.

In this study, we reanalyzed the published RNA-seq datasets (Qian et al., 2020) and characterized the genome-wide APA dynamics in the pomegranate seedlings treated with three kinds of PGRs, paclobutrazol, B9, and mannitol, to decipher the biological significance of APA-mediated mechanisms underlying PGR-induced dwarfing in pomegranate. Furthermore, we also compared the APA regulation to the gene expression changes, with the aim of dissecting the different contributions of transcriptional and post-transcriptional mechanisms to growth repression in pomegranate. Our findings will broaden our understanding of the genetic basis behind the PGR-elicited dwarfed stature of pomegranate and provide a foundation for future molecular-assisted dwarfing cultivation of pomegranate.

Materials and methods

Plant materials and data preprocessing

In our previous study, gene expression was characterized for the seedlings untreated (control group, CK) and treated with each of the three PGRs at different concentrations (paclobutrazol: 6 and 8 mg/L; B9: 6 and 8 mg/L; mannitol: 2.5 and 15 g/L) (Qian et al., 2022). Here, we reanalyzed the 14 RNA-seq datasets (two biological replicates for each scenario), which were deposited in the Gene Expression Omnibus (GEO) database of the National Center for Biotechnology Information (NCBI) under the accession number GSE195722, to investigate genome-wide poly(A) usage dynamics under the PGR treatments over CK. Data preprocessing was performed following the pipelines described in Qian et al. (2022). Briefly, for each dataset, low-quality bases, whose quality score < 20, and adapter contamination were first trimmed from the end of reads using Trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Simultaneously, the reads with either error rate > 0.1 or ambiguous/N bases > 15 were discarded from the dataset. Finally, the sequences with length after trimming < 50 bp were also excluded from the downstream analysis. The clean reads were mapped to the pomegranate reference genome (the soft-seeded pomegranate cultivar “Tunisia”) via HISAT2 (Kim et al., 2019; Luo et al., 2020). The reads uniquely aligned to the genome were extracted and converted into bedgraph format using the sub-command genomecov of the BEDTools suite for downstream analysis (Quinlan and Hall, 2010).

Differentially expressed alternative polyadenylation analysis

The gene model file for the reference genome in 12-column bed format (bed12) was converted from the GTF-format genome annotations file using the UCSC tools, gtfToGenePred and genePredToBed (https://genome.ucsc.edu/). The alignment result in bedgraph format and the gene model file were used as the inputs for the APA dynamics analysis using the APAtrap toolkit (Ye et al., 2018). Specifically, the annotated 3’ untranslated regions (3’ UTRs) were first refined and novel 3’ UTRs or 3’ UTR extensions were detected based on the mapping results of all the samples by the identifyDistal3UTR program. All the putative APA sites, as well as the usage level of APA sites, were predicted using predictAPA with default parameter settings. Differential usage analysis was performed for APA sites between CK and each of the treatment scenarios using the R package deAPA. The genes with an adjusted p-value < 0.05 and percentage difference (PD) ≥ 0.1 were considered to be significantly different in APA site usage between two groups, which were denoted as “differentially expressed APA genes (DAGs)”. The functional importance of the DAGs was assessed by Gene Ontology (GO) enrichment analysis using Fisher’s exact test, where the GO terms with p-value < 0.05 were considered to be significantly overrepresented compared to the genome background.

Prediction of putative microRNA target sites

The majority of APA events occur in 3’ UTRs, that is, producing mRNA isoforms with 3’ UTRs of different lengths. Changes in the length of 3’ UTRs may cause the presence or loss of cis-regulatory elements, and thus pose influences on the stability, nuclear export, and translation efficiency of mRNA (Shen et al., 2008; Di Giammartino et al., 2011; Sun et al., 2012; Tian and Manley, 2017). Here, for each comparison, the DAGs were first grouped into two categories based on the Pearson product moment correlation coefficient r: (1) DAGs with r < 0 were supposed to contain more abundant proximal poly(A) site/shortened 3’ UTR under the treatment than CK, while (2) DAGs with r > 0 were indicated to use more distal poly(A) site/lengthened 3’ UTR in the treatment scenario. For each DAG of each category, the DNA sequence of each APA isoform was extracted from the reference genome using the sub-command fastaFromBed of BEDTools, and the putative microRNA (miRNA) target sites were identified in the 3’ UTR by screening against the collected miRNA sequences in miRBase (Release 21) using the psRNATarget web server (http://plantgrn.noble.org/psRNATarget/). The maximum cutoff of complementary matching score was set to 4.0. The isoforms undergoing 3’ UTR lengthening were supposed to be under the extra regulation of the miRNAs whose target sites were located in the lengthened 3’ end, compared to those with shorter 3’ UTRs.

Comparison between differentially expressed alternative polyadenylation and differentially expressed genes

To further investigate the different regulatory roles of gene transcription and APA in PGR-induced dwarfing, we compared the DAGs to the differentially expressed genes (DEGs) detected in our previous study (Qian et al., 2022). The overlapping between DAGs and DEGs was visualized by a Venn diagram using the draw.pairwise.venn function of the R package VennDiagram. GO enrichment analysis was implemented with a cutoff p-value of 0.05 for the genes from each of the three categories: (1) the genes under the regulation of both differential expression and APA; (2) the genes specifically regulated by differentially expressed APA (DA-specific genes); and (3) the genes specifically regulated by differential expression (DE-specific genes).

Results

PGR-induced alternative polyadenylation changes play a substantial role in regulating pomegranate growth

Compared to CK, exogenous applications of PGRs elicited 289–2,553 DAGs with significant differentiations in APA usage (Figure 1A). Functional enrichment analysis showed that these PGR-responsive APA events were associated with the biological processes of plant growth and development (Figures 1B–D). For instance, the DAGs induced by 8 mg/L B9 were enriched in auxin transport, root development, and maintenance of shoot apical meristem identity (Figure 1B), and mannitol-responsive DAGs were predominantly involved in the GO terms of leaf development and senescence, photomorphogenesis, stomatal movement, and cellular response to osmotic stress (Figure 1C). The application of 6 mg/L paclobutrazol was found to affect growth regulation and cell wall biosynthesis, and DAGs under the 8 mg/L treatment was overrepresented in leaf development (Figure 1D).

FIGURE 1

Figure 1 Numbers of DAGs induced by three PGRs and the corresponding enriched GO terms. (A) Numbers of DAGS identified under the treatment of B9 (blue), mannitol (yellow), and paclobutrazol (red). (B–D) Featur GO terms significantly enriched for the DAGs responsive to B9 (B), mannitol (C), and paclobutrazol (D).

For all the PGRs, the treatment at high concentration could provoke more alterations in APA profiles than that at low level (Figure 1A), which was consistent with their larger effects on suppressing the growth of pomegranate seedlings (Qian et al., 2022). With regard to paclobutrazol, 341 DAGs were shared between the treatments at the two concentrations, which were overrepresented in the regulations of growth rate and leaf senescence (Figure 2A). Comparatively, 241 and 2,212 genes displayed 3’ UTR alterations specifically under 6 and 8 mg/L treatment, respectively (Figure 2A). In particular, the DAGs specifically induced by 6 mg/L paclobutrazol were enriched in leaf development and thylakoid, whereas those responsive to 8 mg/L treatment were involved in cell tip growth and stomatal movement regulation (Figure 2A). Similar concentration-level specificities were also observed in the treatments of B9 and mannitol (Supplementary Figure 1). For instance, the 6 mg/L B9 treatment altered the poly(A) site usage of the genes related to seedling development, shoot apical meristem development, and cell wall thickening, and the DAGs identified under the 8 mg/L treatment were overrepresented in auxin transport, developmental process, and maintenance of shoot apical meristem identity (Supplementary Figure 1B, left panel). Together, these functional specificities of APA events unraveled the dose–response relationships of PGR treatments, which may assist to determine the optimal concentration for PGR application.

FIGURE 2

Figure 2 Overlap of the DAGs across treatments. (A) DAGs overlapped between 6 and 8 mg/L paclobutrazol treatments, and the representative GO terms enriched for each category. Circle size and color represent the significance level (p-value) of enrichment. (B) DAGs overlapped among the treatments of 8 mg/L paclobutrazol, 8 mg/L B9, and 15 g/L mannitol, and the representative GO terms enriched for the DAGs of different categories (highlighted by different colors).

Different alternative polyadenylation regulations were elicited by different PGRs

The APA changes also showed substantial specificities among different PGR treatments (Figure 2B). In total, 1,308, 693, and 109 APA events occurred exclusively when exogenously applied with 8 mg/L paclobutrazol, 8 mg/L B9, and 15 g/L mannitol, respectively, while only 147 events were observed in all these three treatments. Functional enrichment analysis showed that, the commonly changed events were supposed to mainly affect leaf development and senescence, stomatal closure, hydrogen peroxide transmembrane transport, and cellular response to redox state (Figure 2B). Comparatively, the DAGs specifically elicited by 8 mg/L paclobutrazol were enriched in the GO terms of mitochondrial respiratory chain complex I, thylakoid, cellulose biosynthesis, and karrikin response, while those exclusively occurring under the 8 mg/L B9 treatment were overrepresented in auxin transport, xanthophyll biosynthesis, and phototropism (Figure 2B). The biological processes involved in leaf development, tricarboxylic acid transmembrane transport, signal transduction, and hydrogen peroxide biosynthesis were enriched for the specially induced DAGs by 15 g/L mannitol (Figure 2B).

Variations in the expression level of core polyadenylation factors, including polyadenylation machinery components, RNA-binding proteins, and transcription-related process, have been found to regulate APA (Hunt et al., 2012; Tian and Manley, 2017). Here, we first identified the genes encoding the subunits of four types of plant polyadenylation factors, cleavage stimulatory factor (CstF), cleavage and polyadenylation specificity factor (CPSF), poly(A) binding proteins (PABPs), and factor interacting with poly(A) polymerase (FIP1), in the pomegranate genome and compared their expression profiles under each treatment scenario to CK. The results showed that there were several polyadenylation factors significantly differentially expressed in response to the application of PGRs, which may play an important role in APA regulations (Figure 3). It is noteworthy that the expression profiles of polyadenylation apparatuses also displayed ample specificities among different treatments. Only one factor, PABP 7B, was commonly differentially expressed across all the three PGRs. In contrast, the expression of PABP 5 and 6 was specifically altered under the treatment of B9, and the gene that encodes PABP 9 and CPSF subunit 6 had an exclusively differential expression when treated with 8 mg/L paclobutrazol (Figure 3).

FIGURE 3

Figure 3 Heatmap illustrating the log2fold-change in the expression level of polyadenylation factor subunits under different treatments. * indicates the subunits of significantly differential expression under the treatment condition compared to CK.

Alternative polyadenylation and gene transcription variations play a relatively independent role in growth regulation

To explore the different roles of gene transcription and APA dynamics in the regulation of PGR-induced dwarfing in pomegranate, we compared DAGs identified in each treatment to DEGs of the corresponding scenario. When treated with PGRs, most genes were specifically under the regulation of either gene expression or APA (Figure 4 and Supplementary Figures 2, 3). For instance, 24 and 572 out of the 582 and 2,553 DAGs, which accounted for 4.1% and 22.4%, were also differentially expressed under the treatment of 6 and 8 mg/L paclobutrazol, respectively (Figures 4A, D). The DE-specific genes induced by 6 mg/L paclobutrazol were highly represented in cell proliferation, auxin biosynthesis, and brassinosteroid (BR) response (Figure 4B), while the DA-specific genes were predominantly involved in growth regulation, pectinesterase activity, and responses to salicylic acid and strigolactone (Figure 4C). When exposed to 8 mg/L paclobutrazol, genes related to superoxide dismutase activity and root hair elongation were likely to be exclusively differentially expressed compared to CK (Figure 4E), whereas those participating in cellulose biosynthesis, oxidative stress regulation, and responses to salicylic acid and karrikin showed different APA usages (Figure 4F).

FIGURE 4

Figure 4 Comparison between DAGs and DEGs under the treatments of paclobutrazol. (A) Venn diagram illustrating the overlap between DAGs and DEGs when treated with 6 mg/L paclobutrazol. (B, C) GO terms enriched for the DA-specific (B) and DE-specific genes (C) under treatment of 6 mg/L paclobutrazol. (D) Venn diagram illustrating the overlap between DAGs and DEGs when treated with 8 mg/L paclobutrazol. (E, F) GO terms enriched for the DA-specific (E) and DE-specific genes (F) under the treatment of 8 mg/L paclobutrazol.

A similar pattern was also observed under the treatment of B9 and mannitol: only 4.2%–14.2% of DAGs were overlapped with DEGs (Supplementary Figures 2A, E and 3A, E). When treated with 8 mg/L B9, for one example, the DA-specific genes were significantly enriched in the processes of root development, maintenance of shoot apical meristem identity, and auxin transport (Supplementary Figure 2F), while the DE-specific genes were overrepresented in cell wall catabolism and oxidative stress responses (Supplementary Figure 2G). With regard to the application of 2.5 g/L mannitol, the DA-specific genes were highly represented in cellular response to osmotic stress, callose deposition in cell wall, and leaf senescence (Supplementary Figure 3B), and the GO terms related to cell proliferation, growth, and development were found to be enriched for the genes specifically regulated by differential expression (Supplementary Figure 3C). Comparatively, when treated with 15 g/L mannitol, the light-mediated leaf development, leaf senescence, and photoperiodism were mainly regulated by APA events (Supplementary Figure 3F), whereas the genes involved in auxin metabolism, cell development-related programmed cell death, cell wall thickening, secondary shoot formation, superoxide radical removal, and L-ascorbic acid transmembrane transport were largely under the control of different expression (Supplementary Figure 3G). Taken together, these results suggested that, in many scenarios, APA and gene transcription regulate different aspects of the growth and development of pomegranate seedlings and together contribute to the PGR-induced dwarfed stature.

Changes in 3’ UTR length affect microRNA target sites

APA events were also found to substantially modulate gene expression at post-transcriptional and translational levels. Compared to CK, DAGs in the PGR-treated seedlings displayed a global preference for using distal poly(A) sites (Figures 5A, B and Supplementary Figures 4A, 5A). For example, under the treatment of 8 mg/L paclobutrazol, 2,360 DAGs exhibited a higher abundance of the isoforms with longer 3’ UTRs, while only 305 genes used more proximal poly(A) sites (Figure 5B). These lengthened 3’ UTRs were supposed to host more miRNA target sites, which can further modulate mRNA abundance by influencing their stability. Consistently, 65.99%–74.07% of the isoforms using longer 3’ UTRs were inferred to consist of extra miRNA target sites, compared to those using shorter ones, under all the treatment scenarios (Figures 5A, B and Supplementary Figures 4B, 5B). For up to 61 isoforms, more than 10 putative miRNA targets were under the impact of the changes in 3’ UTR length (Figure 5C and Supplementary Figures 4C, 5C). Most of these miRNAs, both constitutive (existing in isoforms with both short and long 3’ UTRs) and lengthened 3’ UTR-specific miRNAs, were identified to function in cleavage of the corresponding mRNA, while 9.00%–12.73% of the miRNAs specific to the extended 3’ UTRs were supposed to inhibit mRNA translation (Figure 5D and Supplementary Figures 4D, 5D). More important, up to 259 DAGs with lengthened 3’ UTRs targeted by miRNAs were significantly downregulated under the treatment of PGRs (Supplementary Figures 4E, 5E, 6). When treated with 8 mg/L paclobutrazol, the miRNA-mediated downregulated genes were overrepresented in the GO terms of developmental growth, lateral root branching, maintenance of shoot apical meristem identity, and cellular response to strigolactone, indicative of their important roles in regulating the growth and development of pomegranate seedlings (Figure 5E).

FIGURE 5

Figure 5 Overview of DAGs with lengthened/shortened 3’ UTR and the putative miRNAs targeting to the lengthened 3’ UTRs under the treatment of paclobutrazol. (A, B) Left panels: 3D volcano plot illustrating the DAGs displaying 3’ UTR lengthening (red) and shortening (blue) when treated with 6 (A) and 8 mg/L paclobutrazol (B). The bar above each volcano plot shows the number of DAGs using more longer (red) and shorter 3’ UTRs (blue). Right panels: Pie charts showing the proportion of isoforms with (blue) or without (gray) putative miRNAs targeting the lengthened area of 3’ UTRs under 6 mg/L (A) and 8 mg/L paclobutrazol treatment (B). (C) Frequency distribution illustrating the number of miRNA target sites identified specifically in lengthened 3’ UTRs across genes, when treated with 6 mg/L (green) and 8 mg/L paclobutrazol (orange), respectively. (D) Numbers of constitutive (existing in isoforms with both short and long 3’ UTRs) or lengthened 3’ UTR-specific miRNAs that were predicted with putative mRNA cleavage (green) and translation inhibition potentials (light green) under two paclobutrazol treatment scenarios, respectively. (E) Representative GO terms enriched for the DAGs that were with miRNAs targeting the lengthened 3’ UTR area and significantly downregulated in response to 8 mg/L paclobutrazol treatment.

Discussion

In the current study, we explored the PGR-induced APA dynamics using the RNA-seq data from our previous study and showed that all the PGR treatments, even at low concentrations, provoked genome-wide alterations in the usages of poly(A) sites (Figure 1A). These changes were found to substantially influence how pomegranate seedlings grew and developed. For example, the poly(A) site usages of the genes involved in auxin transport and growth regulation were significantly altered after the treatments (Figures 1B–D). In particular, one DAG encodes protein PIN-LIKES, which functions as an efflux carrier that mediates the unidirectional auxin flow (Polar auxin transport, PAT) among plant tissues (Křeček et al., 2009). The gene encoding mitogen-activated protein kinase kinase 2 (MKK2) also displayed significantly different APA profiles under the treatment. MKK2, together with mitogen-activated protein kinase 10 (MPK10), forms a module of mitogen-activated protein kinase (MAPK) signaling pathways that serves as a key regulator for PAT in plants (Jagodzik et al., 2018). These alterations in auxin transport may make contributions to the repressed growth and development in pomegranate. Correspondingly, the APA usages of the genes involved in shoot apical meristem identity maintenance and leaf development were also changed in response to PGR treatments (Figures 1, 2B). One of such genes is the calpain-type cysteine protease encoding gene DEK1. Studies in Physcomitrella patens highlighted the important function of DEK1 in controlling the cell fate transition from 2D to 3D growth, where DEK1 knockout in P. patens led to aberrant cell divisions and developmental arrest in buds (Demko et al., 2014; Johansen et al., 2016). Together, these widespread alterations in APA profiles under PGR treatments indicate the substantial significance of post-transcriptional mechanisms in modulating the dwarfed stature of pomegranate.

The APA dynamics display ample specificities among the treatments of different types/concentrations of PGRs (Figures 1, 2), which correspond to their distinct nature. The DAGs induced by 8 mg/L paclobutrazol were particularly associated with the response for karrikins, a type of plant growth regulator that controls plant development (Wang et al., 2020), while the genes with significant APA changes under 8 mg/L B9 treatment were overrepresented in phototropism and xanthophyll biosynthesis (Figure 2B). It is consistent with our previous observations from the transcriptome data that genes responsive to strigolactones, a type of plant signaling compound with similar biochemical properties and physiological activities to karrikins, were specifically downregulated when exposed to 8 mg/L paclobutrazol, whereas those involved in photosynthesis and photosystem II assembly/repair were suppressed by the application of B9 (Qian et al., 2022). Compared to B9 and paclobutrazol, mannitol treatments at both concentrations were found to elicit cellular response to osmotic stress (Figure 1C), corresponding to its specific mechanism that mannitol represses plant growth and development by increasing ambient osmatic pressure and causing drought stress to plants (Bhat and Chandel, 1993); thus, antioxidant reactions were activated to alleviate the oxidative damage. These results indicated that different regulatory mechanisms underlying the pomegranate dwarfing elicited by different PGRs were employed.

In the current study, we found that, in most scenarios, APA and transcriptional regulations are not synchronized and modulate PGR-induced growth repression via different routes, as manifested by both the little overlap between DAGs and DEGs and the differences in the pathways modulated by DAGs and DEGs (Figures 4A, D and Supplementary Figures 2A, E and 3A, E). The transcriptome data revealed that paclobutrazol obviously downregulated the genes related to the tryptophan-independent auxin biosynthetic process (Qian et al., 2022). Comparatively, APA events predominantly affected the polar transport of auxin among tissues. Similarly, when treated with 2.5 g/L mannitol, the cell wall modification process was modulated by the changes in both gene transcription and APA, although in distinct ways (Supplementary Figure 3). In particular, the gene encoding endoglucanase 8 (CEL1), which is a type of cellulose-hydrolyzing enzyme that regulates the cell wall relaxation associated with cell growth and expansion, was significantly downregulated (Tsabary et al., 2003; Shani et al., 2006). Consequently, the suppression of CEL1 would disrupt the differentiations of the plant vascular system and lead to shorter roots and shoots (Tsabary et al., 2003). In contrast, UTP-glucose-1-phosphate uridylyltransferase (also referred to as UDP-glucose pyrophosphorylase, UGPase), which supplies UDP-glucose substrate for the formation of secondary cell wall in plants, displayed substantial poly(A) usage variations. Moreover, Payyavula et al. (2014) showed that the maintenance of UGPase’s function was important for the normal growth of Populus deltoides. These results suggested that both APA and gene transcription make key contributions to the intricate regulatory network underlying the dwarfing stature of pomegranate.

Despite the independent function of APA in regulation, APA is able to modulate transcriptome via influencing the presence or absence of regulatory elements located in 3’ UTRs. Here, the DAGs induced by PGR treatments, especially at high concentrations, displayed a global preference for 3’ UTR lengthening (Figures 5A, B left panel and Supplementary Figures 4A, 5A), which anchor more miRNA target sites than the corresponding shorter 3’ UTRs (Figures 5A, B right panel and Supplementary Figures 4B, 5B). The majority of these extra “burdens” were supposed to suppress gene expression by causing the cleavage or destabilization of mRNA (Figure 5D and Supplementary Figures 4D, 5D). In particular, 259 DAGs with miRNAs specifically targeted in the lengthened 3’ UTRs were significantly downregulated in response to the 8 mg/L paclobutrazol treatment (Supplementary Figure 6). These genes were found to participate in developmental growth, lateral root branching, and maintenance of shoot apical meristem (Figure 5E). Of them, the gene encoding E3 ubiquitin-protein ligase (KEG) is known by its negative regulatory activity of abscisic acid (ABA) signaling, and the suppression of its expression has been shown to retard the growth of A. thaliana (Stone et al., 2006). As another example, the expression of the gene that encodes ammonium transporter 1 member 1 (AMT1;1) was also supposed to be suppressed. AMT1;1 plays an important role in ammonium uptake from soil solution by roots and the subsequent root-to-shoot transport of ammonium; thus, the inhibition of AMT1;1 would lead to nitrogen deficiency and growth defect in pomegranate seedlings (Mayer and Ludewig, 2006). Together, these results highlight the role of APA events in fine-tuning gene expression in response to PGR treatments, which makes a key contribution to the retarded growth and development in the dwarfed pomegranate seedlings.

Conclusion

In this study, we, for the first time, identified and characterized the APA dynamics underlying PGR-elicited dwarfing in pomegranate. Our findings highlight the biological importance of post-transcriptional mechanisms in modulating pomegranate growth and development, which adds a new dimension to the genetic basis of the agronomic trait of pomegranate. However, since our study is mainly based on the prediction from RNA-seq, we might be lacking in power to capture all of the signals and miss some of the true APA events. Thus, in the future, a more comprehensive investigation on the poly(A) usage alterations in pomegranate is essential using efficient technology to measure 3’ UTR dynamics, such as Poly(A) tag sequencing (PAT-seq).

Data availability statement

The datasets we used in this study can be found in online repositories National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) with accession number of GSE195722.

Author contributions

YY and JQ designed the study. XX, MF, and YY collected the data and performed the bioinformatic analyses. XX, YL, XC, JW, JQ, and YY wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (No. 32201420 to YY), the Key Research and Development Projects of Anhui province (No. 202204c06020062 to JQ), and the Production–Education–Research project (AKZY2022110 to JQ).

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2023.1189456/full#supplementary-material

References

Asrey, R., Kumar, K., Sharma, R., Meena, N. K. (2020). Fruit bagging and bag color affects physico-chemical, nutraceutical quality and consumer acceptability of pomegranate (Punica granatum l.) arils. J. Food Sci. Technol. 57, 1469–1476. doi: 10.1007/s13197-019-04182-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhat, S., Chandel, K. (1993). In vitro Conservation of musa germplasm: effects of mannitol and temperature on growth and storage. J. Hortic. Sci. 68 (6), 841–846. doi: 10.1080/00221589.1993.11516422

CrossRef Full Text | Google Scholar

Bourekoua, H., Różyło, R., Gawlik-Dziki, U., Benatallah, L., Zidoune, M. N., Dziki, D. (2018). Pomegranate seed powder as a functional component of gluten-free bread (Physical, sensorial and antioxidant evaluation). Int. J. Food Sci. Technol. 53 (8), 1906–1913. doi: 10.1111/ijfs.13777

CrossRef Full Text | Google Scholar

Colgan, D. F., Manley, J. L. (1997). Mechanism and regulation of mRNA polyadenylation. Genes Dev. 11 (21), 2755–2766. doi: 10.1101/gad.11.21.2755

PubMed Abstract | CrossRef Full Text | Google Scholar

de Lorenzo, L., Sorenson, R., Bailey-Serres, J., Hunt, A. G. (2017). Noncanonical alternative polyadenylation contributes to gene regulation in response to hypoxia. Plant Cell 29 (6), 1262–1277. doi: 10.1105/tpc.16.00746

PubMed Abstract | CrossRef Full Text | Google Scholar

Demko, V., Perroud, P.-F., Johansen, W., Delwiche, C. F., Cooper, E. D., Remme, P., et al. (2014). Genetic analysis of DEFECTIVE KERNEL1 loop function in three-dimensional body patterning in physcomitrella patens. Plant Physiol. 166 (2), 903–919. doi: 10.1104/pp.114.243758

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Giammartino, D. C., Nishida, K., Manley, J. L. (2011). Mechanisms and consequences of alternative polyadenylation. Mol. Cell 43 (6), 853–866. doi: 10.1016/j.molcel.2011.08.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Elkon, R., Ugalde, A. P., Agami, R. (2013). Alternative cleavage and polyadenylation: extent, regulation and function. Nat. Rev. Genet. 14 (7), 496–506. doi: 10.1038/nrg3482

PubMed Abstract | CrossRef Full Text | Google Scholar

Foster, T. M., McAtee, P. A., Waite, C. N., Boldingh, H. L., McGhie, T. K. (2017). Apple dwarfing rootstocks exhibit an imbalance in carbohydrate allocation and reduced cell growth and metabolism. Hortic. Res. 4, 17009. doi: 10.1038/hortres.2017.9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunt, A. G., Xing, D., Li, Q. Q. (2012). Plant polyadenylation factors: conservation and variety in the polyadenylation complex in plants. BMC Genomics 13, 641. doi: 10.1186/1471-2164-13-641

PubMed Abstract | CrossRef Full Text | Google Scholar

Jagodzik, P., Tajdel-Zielinska, M., Ciesla, A., Marczak, M., Ludwikow, A. (2018). Mitogen-activated protein kinase cascades in plant hormone signaling. Front. Plant Sci. 9, 1387. doi: 10.3389/fpls.2018.01387

PubMed Abstract | CrossRef Full Text | Google Scholar

Johansen, W., Ako, A. E., Demko, V., Perroud, P.-F., Rensing, S. A., Mekhlif, A. K., et al. (2016). The DEK1 calpain linker functions in three-dimensional body patterning in physcomitrella patens. Plant Physiol. 172 (2), 1089–1104. doi: 10.1104/pp.16.00925

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Paggi, J. M., Park, C., Bennett, C., Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37 (8), 907–915. doi: 10.1038/s41587-019-0201-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Křeček, P., Skůpa, P., Libus, J., Naramoto, S., Tejos, R., Friml, J., et al. (2009). The PIN-FORMED (PIN) protein family of auxin transporters. Genome Biol. 10 (12), 1–11. doi: 10.1186/gb-2009-10-12-249

CrossRef Full Text | Google Scholar

Luo, X., Li, H., Wu, Z., Yao, W., Zhao, P., Cao, D., et al. (2020). The pomegranate (Punica granatum l.) draft genome dissects genetic divergence between soft-and hard-seeded cultivars. Plant Biotechnol. J. 18 (4), 955–968. doi: 10.1111/pbi.13260

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayer, M., Ludewig, U. (2006). Role of AMT1;1 in NH4+ acquisition in arabidopsis thaliana. Plant Biol. 8 (4), 522–528. doi: 10.1055/s-2006-923877

CrossRef Full Text | Google Scholar

Payyavula, R. S., Tschaplinski, T. J., Jawdy, S. S., Sykes, R. W., Tuskan, G. A., Kalluri, U. C. (2014). Metabolic profiling reveals altered sugar and secondary metabolism in response to UGPase overexpression in populus. BMC Plant Biol. 14, 1–14. doi: 10.1186/s12870-014-0265-8

CrossRef Full Text | Google Scholar

Qian, J., Wang, N., Ren, W., Zhang, R., Hong, X., Chen, L., et al. (2022). Molecular dissection unveiling dwarfing effects of plant growth retardants on pomegranate. Front. Plant Sci. 13, 866193. doi: 10.3389/fpls.2022.866193

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, J., Zhang, X., Yan, Y., Wang, N., Ge, W., Zhou, Q., et al. (2020). Unravelling the molecular mechanisms of abscisic acid-mediated drought-stress alleviation in pomegranate (Punica granatum l.). Plant Physiol. Biochem. 157, 211–218. doi: 10.1016/j.plaphy.2020.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinlan, A. R., Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26 (6), 841–842. doi: 10.1093/bioinformatics/btq033

PubMed Abstract | CrossRef Full Text | Google Scholar

Seleznyova, A. N., Tustin, D. S., Thorp, T. G. (2008). Apple dwarfing rootstocks and interstocks affect the type of growth units produced during the annual growth cycle: precocious transition to flowering affects the composition and vigour of annual shoots. Ann. Bot. 101 (5), 679–687. doi: 10.1093/aob/mcn007

PubMed Abstract | CrossRef Full Text | Google Scholar

Shahamirian, M., Eskandari, M. H., Niakousari, M., Esteghlal, S., Hashemi Gahruie, H., Mousavi Khaneghah, A. (2019). Incorporation of pomegranate rind powder extract and pomegranate juice into frozen burgers: oxidative stability, sensorial and microbiological characteristics. J. Food Sci. Technol. 56, 1174–1183. doi: 10.1007/s13197-019-03580-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Shani, Z., Dekel, M., Roiz, L., Horowitz, M., Kolosovski, N., Lapidot, S., et al. (2006). Expression of endo-1, 4-β-glucanase (cel 1) in arabidopsis thaliana is associated with plant growth, xylem development and cell wall thickening. Plant Cell Rep. 25, 1067–1074. doi: 10.1007/s00299-006-0167-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, Y., Ji, G., Haas, B. J., Wu, X., Zheng, J., Reese, G. J., et al. (2008). Genome level analysis of rice mRNA 3'-end processing signals and alternative polyadenylation. Nucleic Acids Res. 36 (9), 3150–3161. doi: 10.1093/nar/gkn158

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, A., Roychoudhury, A. (2021). Gene regulation at transcriptional and post-transcriptional levels to combat salt stress in plants. Physiol. Plant 173 (4), 1556–1572. doi: 10.1111/ppl.13502

PubMed Abstract | CrossRef Full Text | Google Scholar

Stone, S. L., Williams, L. A., Farmer, L. M., Vierstra, R. D., Callis, J. (2006). KEEP ON GOING, a RING E3 ligase essential for arabidopsis growth and development, is involved in abscisic acid signaling. Plant Cell 18 (12), 3415–3428. doi: 10.1105/tpc.106.046532

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Fu, Y., Li, Y., Xu, A. (2012). Genome-wide alternative polyadenylation in animals: insights from high-throughput technologies. J. Mol. Cell Biol. 4 (6), 352–361. doi: 10.1093/jmcb/mjs041

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, B., Manley, J. L. (2017). Alternative polyadenylation of mRNA precursors. Nat. Rev. Mol. Cell Biol. 18 (1), 18–30. doi: 10.1038/nrm.2016.116

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsabary, G., Shani, Z., Roiz, L., Levy, I., Riov, J., Shoseyov, O. (2003). Abnormalwrinkled'cell walls and retarded development of transgenic arabidopsis thaliana plants expressing endo-1, 4-[beta]-glucanase (cell) antisense. Plant Mol. Biol. 51 (2), 213. doi: 10.1023/A:1021162321527

PubMed Abstract | CrossRef Full Text | Google Scholar

Turrini, F., Boggia, R., Donno, D., Parodi, B., Beccaro, G., Baldassari, S., et al. (2020). From pomegranate marcs to a potential bioactive ingredient: a recycling proposal for pomegranate-squeezed marcs. Eur. Food Res. Technol. 246, 273–285. doi: 10.1007/s00217-019-03339-4

CrossRef Full Text | Google Scholar

Wang, T., Liu, L., Wang, X., Liang, L., Yue, J., Li, L. (2018). Comparative analyses of anatomical structure, phytohormone levels, and gene expression profiles reveal potential dwarfing mechanisms in shengyin bamboo (Phyllostachys edulis f. tubaeformis). Int. J. Mol. Sci. 19 (6), 1697. doi: 10.3390/ijms19061697

CrossRef Full Text | Google Scholar

Wang, L., Xu, Q., Yu, H., Ma, H., Li, X., Yang, J., et al. (2020). Strigolactone and karrikin signaling pathways elicit ubiquitination and proteolysis of SMXL2 to regulate hypocotyl elongation in arabidopsis. Plant Cell 32 (7), 2251–2270. doi: 10.1105/tpc.20.00140

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T., Ye, W., Zhang, J., Li, H., Zeng, W., Zhu, S., et al. (2023). Alternative 3′-untranslated regions regulate high-salt tolerance of spartina alterniflora. Plant Physiol. 191 (4), 2570–2587. doi: 10.1093/plphys/kiad030

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X., Liu, M., Downie, B., Liang, C., Ji, G., Li, Q. Q., et al. (2011). Genome-wide landscape of polyadenylation in arabidopsis provides evidence for extensive alternative polyadenylation. Proc. Natl. Acad. Sci. U.S.A. 108 (30), 12533–12538. doi: 10.1073/pnas.1019732108

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, C., Long, Y., Ji, G., Li, Q. Q., Wu, X. (2018). APAtrap: identification and quantification of alternative polyadenylation sites from RNA-seq data. Bioinformatics 34 (11), 1841–1849. doi: 10.1093/bioinformatics/bty029

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, C., Zhou, Q., Wu, X., Ji, G., Li, Q. Q. (2019). Genome-wide alternative polyadenylation dynamics in response to biotic and abiotic stresses in rice. Ecotoxicol. Environ. Saf. 183, 109485. doi: 10.1016/j.ecoenv.2019.109485

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Z., Hong, L., Li, Q. Q. (2022). Signatures of mRNA alternative polyadenylation in arabidopsis leaf development. Front. Genet. 13, 863253. doi: 10.3389/fgene.2022.863253

CrossRef Full Text | Google Scholar

Zhou, Q., Fu, H., Yang, D., Ye, C., Zhu, S., Lin, J., et al. (2019). Differential alternative polyadenylation contributes to the developmental divergence between two rice subspecies, japonica and indica. Plant J. 98 (2), 260–276. doi: 10.1111/tpj.14209

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Underhill, S. J. R. (2021). Differential transcription pathways associated with rootstock-induced dwarfing in breadfruit (Artocarpus altilis) scions. BMC Plant Biol. 21 (1), 261. doi: 10.1186/s12870-021-03013-6

PubMed Abstract | CrossRef Full Text | Google Scholar