Evolutionary and Integrative Analysis of the Gibberellin 20-oxidase, 3-oxidase, and 2-oxidase Gene Family in Paeonia ostii: Insight into Their Roles in Flower Senescence

Abstract

The brief longevity of tree peony blossoms constrains its ornamental value and economic worth. Gibberellins (GAs) are crucial in the modulation of flower senescence, and GA 20-oxidase (GA20ox), GA 3-oxidase (GA3ox), and GA 2-oxidase (GA2ox) catalyze the synthesis and deactivation of bioactive GAs. In Paeonia ostii, a total of three PoGA20ox, ten PoGA3ox, and twelve PoGA2ox proteins were identified and comprehensively analyzed. The analysis of the gene structures, conserved domains, and motifs revealed structural similarities and variances among the GA20ox, GA3ox, GA2ox-A, and GA2ox-B subfamilies. The synteny analysis indicated a scarcity of collinear blocks within the P. ostii genome, with no tandem or whole-genome duplication/segmental duplications found in PoGAoxs. The investigation into the binding of transcription factors to PoGAox promoters and the assessments of the expression levels suggest that PoGA2ox1 and PoGA2ox8.1 are promising candidate genes implicated in the regulation of floral senescence. Further, Pos.gene61099 (BPC6) and Pos.gene61094 (CIL2) appear to modulate PoGA2ox1 transcription in a positive and negative manner, respectively, while Pos.gene38359 (DDF1) and Pos.gene17639 (DREB1C) likely enhance PoGA2ox8.1’s expression. This study lays a foundation for an in-depth understanding of PoGAox functions and the development of strategies to delay flower senescence in tree peony.

1. Introduction

The tree peony (Paeonia section Moutan DC.), a species indigenous to China, is esteemed for its ornamental beauty, medicinal properties, and the edible oil derived from its seeds (Wang et al., 2019) [1]. There are over two thousand cultivars, distinguished by their sizable blooms, diverse color palette, and robust aroma (Luo et al., 2021; Yang et al., 2020) [2,3]. Its significant ornamental value and meaning of wealth and prosperity have led to its prevalent use in elaborate floral designs and bridal bouquets, signifying its expansive market potential. Nevertheless, a limitation in its horticultural application is the relatively short lifespan of its flowers, lasting only approximately seven days (Liu et al., 2017; Liu et al., 2024) [4,5]. Consequently, developing strategies to delay senescence and extend the blooming period is of paramount importance.

Flower senescence encompasses intricate processes orchestrated by structural and regulatory networks, primarily triggered by plant hormones (Wang et al., 2018; Wu et al., 2017) [6,7]. Gibberellins (GAs), a group of diterpenoid phytohormones, are instrumental across plant development stages and in mitigating stress, significantly influencing plant yield and quality (Toyomasu et al., 2015) [8]. Agricultural practices routinely manipulate GA concentrations for various purposes, including the following: enhancing seedless grapefruit growth; promoting sugarcane stem growth; augmenting fruit set in pome fruits; retarding fruit senescence in citrus; controlling vegetative expansion in canola, apple, and cotton; and terminating bud dormancy in various species (Fuentes et al., 2012; Gao et al., 2023; Sabir et al., 2022) [9,10,11]. GAs are also pivotal in the regulation of key physiological functions, such as seed germination, shoot elongation, leaf development, and fruit and flower senescence (Liao et al., 2019; Wang et al., 2016) [12,13]. Prior research has demonstrated that GA application can extend the lifespan of flowers in species like rose, iris, tobacco, allamanda, freesia, and peony (Ji et al., 2022; Ji et al., 2023a) [14,15]. GAs exhibit antagonism towards ethylene and abscisic acid (ABA) in senescing flowers, with a decline in bioactive GAs accelerating ethylene- and ABA-induced petal senescence in herbaceous peony and rose (Ji et al., 2023b; Lu et al., 2014) [16,17]. A recent study indicated that the PhOBF1 (ocs-element-binding protein 1) gene in petunia acts as a suppressor of ethylene-induced flower senescence by modulating GA biosynthesis (Ji et al., 2023a) [15]. Consequently, administering gibberellin application delays, whereas ABA or ethylene hastens, the senescence of peony petals.

To date, researchers have identified 136 GA molecules, though many serve as biosynthetic intermediates or breakdown products of bioactive GAs, which primarily include GA1, GA3, GA4, and GA7 (Hedden and Thomas, 2012) [18]. GA biosynthesis and deactivation in higher plants involve a three-stage enzymatic process (Wang et al., 2016) [13]. Initially, geranylgeranyl diphosphate is converted into the metabolite ent-kaur-16-ene by the ent-copalyl diphosphate synthase (CPS) and ent-kaurene synthase (KS) enzymes (Cho et al., 2016; Helliwell et al., 2001) [19,20]. In subsequent steps, catalysis by cytochrome P450-dependent monooxygenases, ent-kaurene oxidase (KO), and ent-kaurenoic acid oxidase (KAO), yield GA12 and GA53 from ent-kaur-16-ene (Cho et al., 2016; Helliwell et al., 2001; Wang et al., 2016) [13,19,20]. The final biosynthetic phase generates a variety of GA intermediates and bioactive GAs (GA1 and GA4) from GA12 and GA53 through the early-13-hydroxylation and the non-13-hydroxylation pathways, involving GA 20-oxidases (GA20oxs) and GA 3-oxidases (GA3oxs) in a sequence of oxidation reactions (Hedden and Thomas, 2012; Yamaguchi, 2008) [18,21]. Additionally, bioactive GAs (GA3 and GA7) may derive from GA9 and GA20 via oxidation by GA3ox enzymes in some species (Cho et al., 2016; Hedden and Thomas, 2012) [18,19]. GA2oxs uniquely regulate the degradation pathway, deactivating bioactive GAs (GA1 and GA4) and their immediate precursors (GA20 and GA9) (Olszewski et al., 2002; Rieu et al., 2008; Wang et al., 2016) [13,22,23].

Investigations into GA20ox, GA3ox, and GA2ox sequences have shown their affiliation with the 2-oxoglutarate-dependent dioxygenase (2-ODD) superfamily, sharing significant homology within their functional domains (Kaul et al., 2000; Sabir et al., 2022) [11,24]. The three GAoxs are evidenced by the conserved sequences NXYPXCXXP, three histidine residues, and LPWKET, which may be related to the binding of the common co-substrate 2-oxoglutarate, the Fe²⁺ and GA substrates, respectively (Wang et al., 2020) [25]. GA20ox, in particular, appears to modulate biosynthetic flow as a multicatalytic enzyme. The “green revolution” semidwarf1 rice variety is mutated in a GA20ox that is expressed in shoots but not reproductive tissues, leading to increased grain yields (Sasaki et al., 2002) [26]. For GA3ox, ga4 (GA4 locus) and le (Mendel’s stem length gene) mutants, which cause loss of function to mutate in the GA3ox-encoding gene, were first studied in Arabidopsis and pea, similarly showing that mutants can have a reduced height but not seed production (Chiang et al., 1995; Lester et al., 1997) [27,28]. Different GA2oxs can preferentially act on active GAs or their inactive precursors, with different developmental outcomes. Before the transition to flowering, OsGA2ox1 is expressed just below the shoot meristem, selectively excluding GA1 and GA4 (Sakamoto et al., 2001) [29]. It has been shown in Lolium that GA5 (which is resistant to deactivation by GA2oxs) can move into the meristem despite GA2ox activity (King et al., 2001) [30]. In Arabidopsis thaliana, the initial three cloned GA2ox genes have been annotated as AtGA2ox1, AtGA2ox2, and AtGA2ox3, followed by AtGA2ox4 and AtGA2ox6 (Hedden and Thomas, 2012; Thomas et al., 1999) [18,31]. Moreover, the distinct genes, AtGA2ox7, AtGA2ox8, AtGA2ox9, and AtGA2ox10, which differ from AtGA2ox1~6, hydroxylate C20-GA precursors (Lange et al., 2020; Schomburg et al., 2003) [32,33]. The study of GA3ox and GA2ox genes has shed light on their roles in the feedback regulation of GA biosynthesis and plant stature (Desgagne-Penix and Sponsel, 2008) [34].

With the recent publication of a high-quality genome (Yuan et al., 2022) [35], it has become possible to identify the GA20ox, GA3ox, and GA2ox gene families in tree peony. In this study, they are analyzed for the first time and their expression profiles in gradually opening flowers are shown, revealing the GA oxidase (GAox) genes related to flower senescence. These results provide fundamental information for further study of the function and regulatory mechanisms of these GAox genes and delaying flower senescence in the ornamental period in tree peony.

2. Materials and Methods

2.1. Plant Materials and Sequence Retrieval

The four-year-old cultivar of Paeonia suffruticosa, “Xueyingtaohua”, which was utilized in the study, was cultivated in a greenhouse at Qingdao Agricultural University (Qingdao, China) in a controlled environment at 20 °C with a natural photoperiod of 13 h light/11 h dark. Flowers with 15 cm peduncles were severed when the outer petals began to unfurl and immediately placed in water. Flower stems were recut 1 cm underwater and then transferred within one hour to flasks filled with deionized water. Flowers incubated at 20 °C for durations of 0, 24, 48, 60, and 72 h were in the states of initial opening (S1), half opening (S2), full opening (S3), initial senescence (S4), and senescence (S5), respectively. Petal samples from the S1~S5 flowers were harvested, promptly frozen in liquid nitrogen, and then stored at −80 °C for subsequent analysis.

Genomic and annotation gff3 files for Paeonia ostii were acquired from the China National Gene Bank Database (https://ftp.cngb.org/pub/CNSA/data5/CNP0003098/CNS0560369/CNA0050666/ (accessed on 9 December 2022)) (Yuan et al., 2022) [35]. Arabidopsis thaliana and Vitis vinifera genomic files were retrieved from the TAIR database (http://www.arabidopsis.org (accessed on 9 December 2022)) and Ensembl Plants (http://plants.ensembl.org (accessed on 9 December 2022)), respectively.

2.2. Identification and Phylogenetic Analysis

A BLASTP search (e-value < 1 × 10⁻⁵) was conducted using all GA20ox, GA3ox, and GA2ox proteins from A. thaliana (18) and Camellia sinensis (14) (Pan et al., 2017) [36] as queries to identify potential PoGA20oxs, PoGA3oxs, and PoGA2oxs in P. ostii (Table S1). Following BLASTP, HMM (hidden Markov model) searches with the domain DIOX_N (PF14226) and 2OG-FeII_Oxy (PF03171) were performed using the module Simple HMM search in TBtools (Chen et al., 2020) [37], complemented by inspections of nonredundant protein sequences in the NCBI CDD (https://www.ncbi.nlm.nih.gov/cdd (accessed on 8 January 2024)) and InterPro (http://www.ebi.ac.uk/interpro/ (accessed on 8 January 2024)) using the default settings. The construction of the phylogenetic tree in MEGA X (Kumar et al., 2018) [38] involved aligning protein sequences using MUSCLE and building a maximum likelihood tree with the Poisson model and 1000 bootstrap replicates. Further, the tree was visualized on the ITOL web server (Letunic and Bork, 2019) [39].

2.3. Properties, Gene Structure, and Conserved Motif Analysis

The physical locations and strain analyses of all PoGAox genes were completed using TBtools (Chen et al., 2020) [37] with the genome annotation files. ProtParam (Gasteiger et al., 2005) [40] was used to analyze protein characteristics like number of amino acids (NAA), molecular weight (MW), theoretical isoelectric point (pI), instability index (II), aliphatic index (AI), and grand average of hydropathy (GRAVY). Subcellular localizations were predicted using the NCCDipepCTDCCTDTQSO module via the Plant-mSubP website (Sahu et al., 2020) [41].

The gene structure and motif analyses were facilitated by TBtools (Chen et al., 2020) [37] and MEME v5.1.0 (http://meme-suite.org/tools/meme (accessed on 12 January 2024)), respectively. In MEME discover, the specified parameters were as follow: site distribution was set to “zero or one occurrence per sequence”, and the maximum number of motifs was 10. The results were compiled and visualized using TBtools (Chen et al., 2020) [37].

2.4. Duplication and Synteny Analysis of PoGAox Genes

Duplication and synteny analyses were carried out by comparing whole genome protein sequences of P. ostii, A. thaliana, and V. vinifera using BLASTP, with paralogous genes examined via MCScanX in TBtools. Gene duplication events were identified based on two established criteria: the sequence similarity among the gene pairs must exceed 75%, and the span of the homologous region must constitute more than 75% of the longer gene sequence (Gu et al., 2002) [42]. The visualizations of synteny, duplication, and Ka/Ks calculations also utilized TBtools (Chen et al., 2020) [37].

2.5. Cis-Element Prediction and Expression Analysis

Cis-element prediction in PoGAox gene promoter regions (i.e., 1.5 kb upstream regions starting from translation starting sites) were performed with PlantCARE (Lescot et al., 2002) [43]. The expression levels of PoGAoxs in gradually opening flowers (S1~S5) were obtained from the RNA-seq of “Xueyingtaohua” flowers (unpublished) and displayed using a heatmap. Furthermore, the promoters of key genes were submitted to the website PlantTFDB (http://planttfdb.gao-lab.org/ (accessed on 15 January 2024)), and the potential TFs (transcription factors) binding to their promoters (p-value < 10 × 10⁻⁶) were predicted in A. thaliana.

Furthermore, quantitative real-time PCR (qRT-PCR) was employed to validate the expressions of key PoGA2oxs in S1~S5 flowers. Total RNA extractions were performed utilizing an OminiPlant RNA Kit (DNase I)(Jiangsu Cowin Biotech Co.,Ltd., Taizhou, China). These extractions served as the foundation for the subsequent template cDNA synthesis, executed through the use of a HiFiScript cDNA Synthesis Kit (Jiangsu Cowin Biotech Co.,Ltd., Taizhou, China). The qRT-PCR analyses were conducted on the Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, USA), ensuring a minimum of three replicates for each reaction. The expression levels of the target genes were calculated employing the ΔΔCt method. The specific qRT-PCR primers utilized for PoGA2oxs and Ubiquitin, which served as the internal control, are detailed in Table S2.

3. Results

3.1. Identification and Characterization of PoGAox Genes

In the quest to pinpoint PoGAoxs within P. ostii, the protein sequences of all GA20oxs, GA3oxs, and GA2oxs in A. thaliana and Camellia sinensis (Pan et al., 2017) [36], served as the queries in the BLASTP search (Table S1). Next, an HMM search using the DIOX_N (PF14226) and 2OG-FeII_Oxy (PF03171) domains was performed. The subsequent analysis of the yielded nonredundant protein sequences involved validation using CDD and InterPro. Proteins devoid of the domains and not conforming to the annotation were excluded, culminating in the recognition of 3 PoGA20ox, 10 PoGA3ox, and 12 PoGA2ox proteins in P. ostii (Table S3). The study proceeded to dissect the gene subfamilies and phylogenetic relationships by constructing a comprehensive phylogenetic tree encompassing all PoGAox, AtGAox, and OsGAox protein sequences. This phylogenetic tree segregated PoGAoxs into four subfamilies—GA20ox, GA3ox, GA2ox-A, and GA2ox-B—and each gene was renamed in alignment with its homologous AtGAox genes (Figure 1; Table S3). The result show that GA20ox and GA3ox genes were clustered together, while GA2ox genes were divided into the two following categories: GA2ox-A subfamily, containing GA2ox1~6, and GA2ox-B subfamily, containing GA2ox7~10.

The predictive analysis extended to various characteristics of the PoGAox proteins, including NAA, MW, pI, II, AI, GRAVY, and subcellular localizations, all of which are tabulated in Table S4. The MW, which is commonly and directly proportional to the NAA, for PoGA20ox1, PoGA20ox2, and PoGA20ox5 proteins were 37.26 kDa, 39.87 kDa, and 43.40 kDa, respectively. The majority of PoGA3oxs fell within 28.02–40.89 kDa, except for PoGA3ox1.6 and PoGA3ox3 at 18.76 kDa and 21.07 kDa, respectively. The majority of PoGA2ox proteins spanned 30.01 kDa to 43.13 kDa, with specific proteins like PoGA2ox9.2 and PoGA2ox9.4 registering at 18.07 kDa and 26.41 kDa, respectively. The pI values spanned broad ranges, as follows: 6.03–6.86 for GA20oxs, 6.36–9.08 for GA3oxs, 5.33–8.37 for GA2ox-A proteins, and 4.65–8.46 for GA2ox-B proteins. Additionally, most PoGAoxs exhibited GRAVY values below zero, indicative of their predominantly hydrophilic nature, with the exception of PoGA3ox1.8 and PoGA2ox9.2, which are hydrophobic. On the basis of the instability index (II) values, it was observed that all three PoGA20oxs were likely to be stable, with II values ranging from 30.73 to 35.82. In the GA3ox subfamily, PoGA3ox1.1, PoGA3ox1.3, PoGA3ox1.5, and PoGA3ox1.8 were stable (II < 40), and the remaining six GA3oxs were predicted to be unstable (II > 40). The majority of PoGA2oxs are likely to be unstable, except PoGA2ox2, PoGA2ox9.1, and PoGA2ox9.2.

Regarding subcellular localization, the prediction indicated that the majority of PoGAoxs (72.0%) were localized in the plastid, suggesting their involvement in plastid-related functions. In the remaining proteins, PoGA20ox2, PoGA3ox1.7, and PoGA2ox8.1 were located in endoplasm, PoGA3ox1.8 and PoGA2ox9.2 were found in cell membrane, and PoGA3ox1.9 and PoGA2ox9.1 were predicted in the nucleus.

3.2. Gene Structure, Conserved Motifs, and Domains of PoGAoxs

The gene structure analysis of the PoGAox genes revealed some interesting patterns (Figure 2). The majority of PoGAox genes had 2–3 exons, indicating that they were related genes. Furthermore, the PoGAoxs in the same subfamily tended to contain the same number of exons. In detail, all PoGA20oxs and PoGA2ox1~6 (GA2ox-A subgroup) had three exons, whereas the majority of PoGA3oxs had two exons, except PoGA3ox1.7 (four exons) and PoGA3ox3 (one exons). In the GA2ox-B subfamily, there were variations in the gene structure, with PoGA2ox8.1, PoGA2ox8.2, PoGA2ox9.1, PoGA2ox9.2, PoGA2ox9.3, and PoGA2ox9.4 having five, two, three, one, three, and one exons, respectively.

In the identification using HMMER search, the proteins with the domain DIOX_N (PF14226, the highly conserved N-terminal region of 2-oxoglutarate/Fe(II)-dependent dioxygenase proteins) and 2OG-FeII_Oxy (PF03171, the 2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase to exert oxidoreductase activity) were regarded as candidates GAoxs, revealing that all GAoxs had DIOX_N and 2OG-FeII_Oxy domains. Moreover, the conserved domain analysis showed that the PLN02276, PLN02254, PLN02156, and PLN02984 superfamilies, which all contain the DIOX_N and 2OG-FeII_Oxy domains, presented in all proteins of the GA20ox, GA3ox, GA2ox-A, and GA2ox-B subfamilies, respectively (Figure 2). Further analysis using the MEME web server identified 10 conserved motifs containing 29 to 50 amino acids in PoGAoxs (Figure 2; Table S5). These motifs were found in varying numbers in different PoGAoxs, and some motifs were specific to certain subfamilies, such as motifs 4 and 5 in the GA3ox subfamily, motifs 8 and 9 in the GA2ox-A group, and motif 10 in the GA2ox-B group. These specific motifs may contribute to the functional diversity and specialization within different subfamilies of PoGAoxs.

3.3. Chromosomal Location and Synteny of PoGAox Genes

The synteny analysis provides insight into the evolutionary relationships and function between PoGAox gene family and their collinear genes in A. thaliana and V. vinifera. First, the physical location of all PoGAoxs in chromosomes was analyzed and displayed (Figure 3A). As shown, all chromosomes had PoGAox genes, and three PoGA20ox genes were in chromosomes Chr01, Chr04, and Chr05, respectively. Moreover, there were seven PoGA3ox (PoGA3ox1.1~1.7) densely clustered in a subregion of the Chr05. Then, all of the representative proteins of P. ostii, A. thaliana, and V. vinifera were BLASTP analyzed and mapped to genomes for the synteny and evolutionary relationships (Figure 3B,C). Among P. ostii, the syntenic blocks (gray line areas) were few, and the PoGAox gene family exhibited a surprising lack of collinear genes within its own genome (Figure 3B). However, when comparing the P. ostii genome to the genomes of A. thaliana and V. vinifera, a higher number of syntenic blocks were observed (Figure 3C). Specifically, three PoGAox genes (PoGA20ox1, PoGA2ox2, and PoGA2ox4) showed collinearity (blue lines) with GAoxs in both A. thaliana and V. vinifera; the other three PoGAox genes, PoGA3ox1.4, PoGA3ox1.2, and PoGA2ox9.3, presented only collinearity with VvGAoxs (Figure 3C; Table S6).

Gene duplication events play a critical role in the diversification of plant species and the expansion of gene families, and they can be categorized into four types: whole-genome (WGD)/segmental duplication, tandem duplication, and proximal and dispersed duplications. Proximal duplications may arise either from tandem duplications or from small-scale transpositions, or they may involve the insertion of additional genes. Dispersed duplications suggest genes that likely originated through transposition, such as “conservative transposition”, “nonreplicative transposition”, or “replicative transposition”. The current analysis revealed that among the PoGAox genes, most were classified as dispersed duplications, and only PoGA3ox1.1~1.7 were proximal duplications (Table S7). Notably, the study found no instances of WGD/segmental duplication and tandem duplication within this gene family.

In addition, strict criteria for identifying gene duplication, as established by Gu et al. (2002) [42], include a sequence similarity threshold and a requirement that the similar sequence length exceed 75% of the length of the longer sequence (Gu et al., 2002) [42]. Employing these conditions, only six PoGAoxs (forming seven pairs; colored lines and gene names) were pinpointed as products of gene duplications (Figure 3B; Table S8). In detailed, PoGA2ox6.1 and PoGA2ox6.2 came from one gene duplication and the four genes PoGA3ox1.1~1.4 were duplicated with each other. Further, the Ka/Ks values of these seven gene pairs were calculated and shown to be below zero, suggesting they likely had suffered purify selection (Table S8).

3.4. Cis-Element Prediction in PoGAox Promoters

To delve into the expression regulations of PoGAox, a detailed examination of the cis-elements within the promoter regions was performed. Promoter analyses spotlighted the presence of cis-elements that are responsive to hormonal signals (auxin (Aux), jasmonic acid/methyl jasmonate (MeJA), salicylic acid (SA), and gibberellin (GA)) and to various stresses including heat, drought, low-temperature, stress, and anoxic conditions. In detailed, more than half (64%) of PoGAoxs had MeJA- and SA-responsive cis-elements, and only three PoGAoxs, including PoGA20ox2, PoGA3ox1.1, and PoGA2ox6.3, had no STRE (stress-response element) in their promoters. These findings suggest that PoGAoxs may be involved in the regulation of hormonal signaling and stress, especially MeJA, SA, and stress related to the STRE. Surprisingly, not all GAox promoters contained cis-elements responsive to GA, and the GA-responsive elements in the GAox promoters were much fewer than those of SA and JA. In addition, cis-elements pertinent to specific tissues or processes—like light response, zein metabolism, flavonoid biosynthesis, cell cycle, circadian rhythm, endosperm development, meristem activity, and anaerobic response—were identified in promoters (Figure 4), indicating that PoGAoxs likely play a role in the mentioned biological processes. Furthermore, the MYB- and MYC-binding elements were the most conserved elements in the PoGAox promoters.

3.5. Expression and Regulation of PoGAox Genes in Flowers

The exploration of PoGAox expression patterns during flower senescence was facilitated by analyzing RNA-seq data (unpublished) from “Xueyingtaohua” flowers in stages S1~S5 (initial opening, half opening, full opening, initial senescence, and senescence, respectively). As displayed in Figure 5A, only PoGA2ox1 and PoGA2ox8.1, which both catalyze GA degradation, had relatively higher expression levels in gradually opening flowers. Further, the quantitative real-time PCR (qRT-PCR) results show that PoGA2ox1 had higher expression levels in S4~S5 compared with other periods, while PoGA2ox8.1 had relatively high levels in S2~S4, suggesting PoGA2ox1 and PoGA2ox8.1 are likely involved in flower senescence (Figure 5B).

To understand the regulatory mechanisms of PoGA2ox1 and PoGA2ox8.1 expressions, TFs with potential binding to their promoters were predicted in A. thaliana (Figure 5C; Table S9). It was found that three BBR/BPC family transcription factors, including AtBPC1, AtBPC5, and AtBPC6, all bind to promoters of PoGA2ox1 and PoGA2ox8.1 (proPoGA2ox1 and proPoGA2ox8.1), while three DREB family TFs (DREB1C, DREB1D, and DDF1) all specifically recognized the proPoGA2ox8.1. Moreover, another five TFs (SOC1, RGA1, TCX2, BBM, and CIL2 (CIB1-like protein)) were found to bind to proPoGA2ox1. To confirm the regulatory relationships, the expressions of peony genes homologous to these 11 TFs in the S1~S5 flowers were obtained from the RNA-seq data (Figure 5D). There were four peony TFs, including Pos.gene61094 (CIL2), Pos.gene61099 (BPC6), Pos.gene38359 (DDF1), and Pos.gene17639 (DREB1C), that had higher expressions in flowers in S1~S5. In the qRT-PCR results, Pos.gene61099 (BPC6) exhibited an increased expression in the S4 flowers, whereas Pos.gene61094 (CIL2) showed elevated expression levels in the S1~S3 flowers, indicating that Pos.gene61099 (BPC6) and Pos.gene61094 (CIL2) probably act as positive and negative regulators of PoGA2ox1, respectively (Figure 5B). Additionally, Pos.gene38359 (DDF1) and Pos.gene17639 (DREB1C), both of which were more abundantly expressed in the S3 flowers, may be positively correlated with the expression of PoGA2ox8.1 (Figure 5B).

4. Discussion

Tree peony is a popular ornamental plant worldwide distinguished by their sizable blooms, diverse color palette, and robust aroma (Luo et al., 2021; Yang et al., 2020) [2,3]. However, the short lifespan of its flowers limits its horticultural application, floral designs and bouquets, damaging their economic value (Liu et al., 2017; Liu et al., 2024) [4,5]. Gibberellins are instrumental across plant development stages and in mitigating stress, such as seed germination, shoot elongation, leaf development, and fruit and flower senescence (Liao et al., 2019; Toyomasu et al., 2015; Wang et al., 2016) [8,12,13]. In the GA pathway, GA20oxs, GA3oxs, and GA2oxs are the key enzymes that catalyze the synthesis and inactivation of GA intermediates and bioactive GAs, regulating the dynamic equilibrium of GA (Cho et al., 2016; Hedden and Thomas, 2012; Olszewski et al., 2002; Rieu et al., 2008; Wang et al., 2016; Yamaguchi, 2008) [13,18,19,21,22,23]. In the current study, a total of 3 PoGA20oxs, 10 PoGA3oxs, and 12 PoGA2oxs were identified in tree peony (Figure 1; Table S3), and an in-depth analysis of their protein characteristics, chromosomal locations, synteny relationships, and gene expression patterns was conducted to detect the key genes related to flower senescence.

The functional classification of genes is reflected in their subfamilies to some extent (Cheng et al., 2022; Yuan et al., 2021) [44,45]. On the basis of the phylogenetic tree, PoGA20oxs and PoGA3oxs were grouped into one subfamily, while PoGA2oxs were divided into two subgroups: GA2ox-A subfamily (GA2ox1~6) and GA2ox-B subfamily (GA2ox7~10) (Figure 1). In Arabidopsis thaliana, AtGA2ox1~6 hydroxylated C19-GA substrates, while AtGA2ox7~10 catalyzed C20-GA precursors (Hedden and Thomas, 2012; Thomas et al., 1999) [18,31] (Lange et al., 2020; Schomburg et al., 2003) [32,33], displaying the functional classification of members in the GA2ox-A and GA2ox-B subfamilies. It was also found in the FAD family that, for instance, the members of the FAD2 and FAD3 subfamilies were responsible for lipid desaturation in the endoplasmic reticulum, while FAD6, FAD7, and FAD8 members were involved in desaturation reactions in the chloroplasts (Cheng et al., 2022) [44]. The functional specialization and subfamily classification were also reflected in the diversity in the gene structure, domains and motifs (Cheng et al., 2022; Yuan et al., 2021) [44,45]. As shown in tree peony GA2oxs, all PoGA2ox1~6 had three exons, whereas the members of the GA2ox-B subfamily exhibited variety in their gene structures (Figure 2). The domains PLN02156 and PLN02984 were, respectively, conserved in the GA2ox-A and GA2ox-B subfamilies. For motif analysis, motifs 8 and 9 were specific to the GA2ox-A group, and motif 10 was conserved in the GA2ox-B subfamily (Figure 2). In addition, the gene structures and motifs varied in the different PoGAoxs, indicating that these differences might contribute to the functional diversity and specialization within or among different subfamilies of PoGAoxs. Furthermore, the similarities and differences in motif distributions and gene structures within subfamilies further support the accuracy of the evolutionary tree and provide bases for understanding the functional similarities and differences among subfamilies.

Duplication plays a vital role in the expansion of a genome and gene family, and there are four types of gene duplication events, including whole genome duplication, segmental duplication, tandem duplication, and transposition events (Paterson et al., 2010) [46], which can be identified and categorized into WGD/segmental duplication, tandem duplication, dispersed duplication, and proximal duplication in TBtools (Chen et al., 2020) [37]. Commonly, tandem duplication along with WGD/segmental duplication are pivotal in augmenting the diversity of gene families, and tandem duplication often leads to clusters of genes with similar sequences and functions (Panchy et al., 2016; Tyagi et al., 2021) [47,48]. However, there were a few collinear blocks, and only 5439 (7.69% of 70700) tandem and 1154 (1.63% of 70700) WGD or segmental genes were identified in P. ostii, and neither tandem nor WGD or segmental events were found in PoGAoxs in this study. PoGA3ox1.1~1.7 were regarded as proximal genes (arising from transposition or tandem duplication and insertion of other genes), while other PoGAoxs were designated as dispersed genes (arising from transposition). These results are contrary to the general rule that WGD and tandem duplication contribute a lot in genome and gene family expansion. Surprisingly, previous research revealed that P. ostii was likely to experience explosive transposition of the long terminal repeat (LTR), instead of pedigree-specific WGD, thus forming its gigagenome (Yuan et al., 2022) [35], which supports our findings in this study. In addition, there were more collinear pairs of PoGAoxs among P. ostii–V. vinifera than P. ostii–A. thaliana, which was consistent with the fact that P. ostii and V. vinifera have greater segments from an ancestral reconstructed karyotype (Yuan et al., 2022) [35].

The variable expression patterns of genes, to a certain degree, indicate functional disparities (Ye et al., 2019) [49]. Consequently, the cis-elements within promoters associated with transcriptional regulation could also reflect the gene function. Previous research has revealed that the JA/SA signaling pathways jointly downregulate genes associated with gibberellin, redirecting growth resources towards defense (Zhang et al., 2020) [50]. In V. vinifera, the transcript levels of VvGA2ox1a, VvGA2ox10 and VvGA2ox11 were upregulated significantly after MeJA treatment (Gao et al., 2022) [51]. Similarly, CaGA2ox1 was induced and upregulated expression by salicylic acid treatment in pepper (Lee et al., 2012) [52]. In A. thaliana, It was found that the declines in GA3ox1 transcriptions were rescued in ibuprofen-treated or JA signaling-impaired plants, revealing that GA biosynthesis in Ph. liquidambaris-inoculated roots was antagonized by JA (Zhang et al., 2021) [53]. In this study, 64% of PoGAoxs had MeJA- and SA-responsive cis-elements, and only three PoGAoxs had no STRE, which responds to stress, in their promoters, suggesting that PoGAoxs may be involved in the responses of MeJA, SA, and stress or in cross-talk between JA/SA and GA signals. Remarkably, not all GAox promoters harbor cis-elements responsive to GA, and the GA-responsive elements within GAox promoters are markedly fewer than those associated with SA and JA, suggesting the JA and SA might exert a more direct regulatory effect on GAoxs than GA itself, as mentioned above (Gao et al., 2022) (Lee et al., 2012) (Zhang et al., 2021) [51,52,53], or the modulation of GAoxs could potentially necessitate the involvement of specific transcription factors, such as BBR/BPC and DDF1 (Sun et al., 2018; Magome et al., 2008; Wang et al., 2016) [13,54,55]. Furthermore, the distinct temporal and spatial expression profiles of genes yield critical insights into the timing and location of their roles in growth and developmental processes (Tong et al., 2013) [56]. In the qRT-PCR results of this study, PoGA2ox1 and PoGA2ox8.1 were highly expressed in flowers in S4~S5 and S2~S4, respectively, suggesting that PoGA2ox1 and PoGA2ox8.1 are likely involved in flower senescence (Figure 5B).

In addition, the TF regulation prediction results for PoGA2ox1 and PoGA2ox8.1 and their expression patterns showed that Pos.gene61099 (BPC6) and Pos.gene61094 (CIL2) probably activated and decreased the expression of PoGA2ox1, respectively, and both Pos.gene38359 (DDF1) and Pos.gene17639 (DREB1C) likely promote the transcription of PoGA2ox8.1 (Figure 5B,C). In a previous study, basic pentacysteine (BPC) proteins have been identified to bind specifically to GA enrichment regions within gene promoters (Monfared et al., 2011) [57]. Furthermore, Class II BPC6 proteins are known to modulate developmental processes in plants, such as the flowering period (Hecker et al., 2015) [58]. Similarity, it was also predicted that BBR/BPC proteins directly bind to the CsGAox promoters in cucumber (Sun et al., 2018) [55]. These reports support our previous finding that Pos.gene61099 (BPC6) was likely a regulator of PoGA2ox1. A previous study also demonstrated that CIB (cryptochrome2 (CRY2)-interacting bHLH) proteins could form heterodimers and interact with CRY2 to directly activate the transcription of FT (Flowering Locus T), resulting in the regulation of flowering time in Arabidopsis (Liu et al., 2013) [59]. Moreover, CIB1 and PHR2 can form a complex to activate the transcription of the GA-receptor GID1B, which is vital in the integration of plant age and photoperiod signaling to harmonize signal perception and internal GA concentrations in Chrysanthemum morifolium (Zhao et al., 2023) [60]. However, the direct regulation of PoGA2ox1 via Pos.gene61094 (CIL2) should be further investigated. The candidate regulator Pos.gene38359 (DDF1) and Pos.gene17639 (DREB1C) of PoGA2ox8.1 were both DREB1/CBF transcription factors belonging to the ERF/AP2 family. In A. thaliana, DDF1 overexpression resulted in increased stress tolerance, dwarfism, and postponed flowering, and these phenotypic effects were reversed upon application of exogenous gibberellin, highlighting the role of GA pathway in DDF1-mediated regulation (Kang et al., 2011) [61]. Furthermore, DDF1 has been documented to interact with DRE motifs and directly initiate the transcription of GA2ox7, resulting in the suppression of growth (Magome et al., 2008; Wang et al., 2016) [13,54]. Additionally, constitutive expression of CBF1 could promote the transcription of GA2oxs and a decrease in GA content, resulting in the accumulation of the growth inhibitory protein DELLAs, which leads to dwarfism and delayed flowering (Achard et al., 2008) [62]. Therefore, Pos.gene38359 (DDF1) and Pos.gene17639 (DREB1C) are very probably involved in GA signal regulation by activation of the transcription of PoGA2ox8.1.

5. Conclusions

The study conducted a comprehensive examination of the GAox gene family in Paeonia ostii utilizing a suite of analytical approaches, including phylogenetic, gene structure, conserved domains and motifs, collinearity, duplication events, promoter cis-element, regulatory factor prediction, and gene expression pattern analyses. A total of 3 PoGA20ox, 10 PoGA3ox, and 12 PoGA2ox proteins in P. ostii were identified and renamed according to their similar genes in A. thaliana. The gene structure, conserved domains, and conserved motifs revealed structural differences between the subfamilies of GA20ox, GA3ox, GA2ox-A, and GA2ox-B, which is reflected in the functional classification of the genes to some extent. The synteny analysis revealed that only a few collinear blocks were exhibited among the P. ostii genome and no tandem or WGD/segmental events were found in PoGAoxs, which support the role of the explosive transposition of the long terminal repeat (LTR) in forming a gigagenome. The prediction of the TFs binding to promoters and the expression level analyses suggest that PoGA2ox1 and PoGA2ox8.1 are candidate genes related to flower senescence; Pos.gene61099 (BPC6) and Pos.gene61094 (CIL2) probably regulate PoGA2ox1 transcription positively and negatively, respectively; Pos.gene38359 (DDF1) and Pos.gene17639 (DREB1C) likely enhance the expression of PoGA2ox8.1. This research establishes the groundwork for the further elucidation of PoGAox functions and development of strategies to delay flower senescence in tree peony.