HsGA20ox1, HsGA3ox1 and HsGA2ox1 are Involved in Endogenous Gibberellin Regulation within Heracleum sosnowskyi Ovaries after GA3 Treatment

1. Introduction

Gibberellins (GAs) are a large group of plant hormones, encompassing 136 chemically similar compounds found in plants, fungi, and bacteria [1,2]. They belong to diterpene carboxylic acids with a tetracyclic ent-gibberellane skeleton, composed of 19–20 carbon atoms [3]. Biologically active GAs act as endogenous plant growth and development regulators, transmitting internal and environmental signals within the organism [1]. They stimulate cell elongation, cell division, regulate seed germination, control the transition between juvenile and generative plant phases, induce flowering, influence sex determination, and play a role in fruit setting and growth processes [4,5]. However, only a small group are biologically active: gibberellin A₁ (GA₁), gibberellin A₃ (GA₃), gibberellin A₄ (GA₄), and gibberellin A₇ (GA₇) [1]. Interestingly, GA₁ and GA₄ are the dominant forms in plant tissues [6,7,8]. GA₃ and GA₇ metabolism in plants is not a common process, and these GA species are present in trace amounts, however, microorganisms can metabolize GA₂₀ to GA₃, or GA₇ directly from GA₄, involving GA desaturases [9,10]. For plants, it is crucial to regulate GA levels precisely; however, the chemical structure specificity of 1,2-unsaturated GAs (GA₃ and GA₇) obstructs their inactivation in plants by 2β-hydroxylation, which results in their prolonged activity and may explain their low amounts in plants [2,10].

GA biosynthesis has been studied in many plant species, and it occurs in three stages according to subcellular localization, where involved enzymes act [4,11]. The first stage of GA biosynthesis takes place in plastids, where geranylgeranyl diphosphate is converted to ent-kaurene through an intermediate, catalysed by terpene synthases: ent-copalyl diphosphate synthase (CPS) and ent-kaurene synthase (KS) [12,13]. Subsequently, in the outer plastid membrane and endoplasmic reticulum, cytochrome P450 monooxygenases catalyse oxidation reactions, during which ent-kaurene is converted through four intermediates to GA₁₂ [14,15]. The third stage of synthesis occurs in the cell cytoplasm, where the common precursor GA₁₂ is oxidised by two distinct soluble 2-oxoglutarate-dependent dioxygenases (2-ODDs) – GA 20-oxidase (GA20ox) and GA 3-oxidase (GA3ox) – leading to the biologically active GA₁ or GA₄ [11,16,17]. To keep optimal levels of active GAs species within tissues, the inactivation process is essential for plant growth and development. The main mechanism involves two groups of 2-ODDs enzymes, differentiated by their substrate specificity: C₁₉-GA 2-oxidase and C₂₀-GA 2-oxidase, both of which utilise 2β-hydroxylation to catabolize precursors and biologically active GAs into inactive products [2,18,19,20].

The 2-ODD superfamily contains four enzyme subfamilies (GA20ox, GA3ox, C₁₉-GA2ox and C₂₀-GA2ox), involved in GA biosynthesis and catabolism, which are encoded by small multigene families and are conserved between multiple species [19,21,22,23]. Genes of GAoxs subfamilies were first isolated in model plant species: GA20ox in pumpkin [24], GA3ox in Arabidopsis thaliana [25], GA2ox in runner bean (Phaseolus coccineus) and in A. thaliana [26]. To date, GAoxs genes have been identified in many commercially important plants such as: cucumber [18], watermelon [27], tomato [28,29,30], maize [31], rice [32,33], peach [34], grapevine [35], wild cherry, wild strawberry [36], breadfruit [37] and others. With the identification of GAoxs enzymes and their genes, the mechanism of fruit-set became better understood, which opened opportunities to manipulate fruit-set and induce parthenocarpy [28,29,38,39].

H. sosnowskyi is a noxious invasive weed species enlisted into Invasive Alien Species of Union concern and national lists of invasive species in many EU countries [40,41]. H. sosnowskyi propagates only by seeds and after bearing them, eventually dies [42,43]. Interestingly, Pastinaca sativa (Apiaceae) is adapted to cope with some pests by producing seedless fruits. This strategy enables it to regulate outbreaks of herbivores [44]. Previous studies have indicated an exogenous GA₃ effect in inducing seedlessness in H. sosnowskyi [45,46]. However, to date, analysis of endogenous GAs and related genes has not been available. With the recent article on the H. sosnowskyi genome, GAoxs and gene analysis have become accessible in this invasive species [47].

In this study, we identified 27 putative GAoxs enzymes in the H. sosnowskyi genome. Phylogenetic analysis clustered these proteins into three subgroups: GA20ox, GA3ox and C₁₉-GA2ox. We identified the effects of exogenous GA₃ on the expression of HsGA20ox1, HsGA3ox1 and HsGA2ox1, as well as on the endogenous GA profiles in gradually opening H. sosnowskyi flowers. These results open opportunities to further analyse the role of GAs in H. sosnowskyi fruit-set mechanism and develop invasion control strategies.

2. Results

2.1. Genome-Wide Identification and Analysis of GA-Oxidase Genes in Heracleum sosnowskyi

To identify GAox genes from GA20ox, GA2ox and GA3ox subfamilies in Heracleum sosnowskyi, the amino acid sequences of GA20ox, GA2ox and GA3ox from Arabidopsis thaliana and Oryza sativa were used as a reference in BLASTp program against the H. sosnowskyi genome. A total of 32 putative protein homologues were obtained. However, 5 protein candidates did not possess 2OG-Fe(II) oxygenase (PF03171) or non-haem dioxygenase N-terminal (PF14226) domains, as confirmed by InterPro Pfam (available online: http:// https://www.ebi.ac.uk/interpro) (Table S1). Thus, we concluded that H. sosnowskyi contains 27 GAox family members. We identified and renamed GA20ox, GA2ox and GA3ox genes in the H. sosnowskyi genome: 9 HsGA20ox1-9, 13 HsGA2ox1-13 and 5 HsGA3ox1-5 (Table 1).

The CDS lengths of the HsGA20ox, HsGA2ox and HsGA3ox genes varied from 978 (HsGA20ox7) to 1221 nt (HsGA20ox4), from 945 (HsGA2ox8) to 1665 nt (HsGA2ox2) and from 1020 (HsGA3ox5) to 1068 nt (HsGA3ox2), respectively. In HsGAox families, amino acid residue ranges slightly differed: HsGA20ox (325-406 aa), HsGA2ox (314-554 aa) and HsGA3ox (339-355 aa). The theoretical molecular weight ranged from 36.81 to 46.14 kDa (HsGA20ox), for HsGA2ox from 35.39 to 61.92 kDa and for HsGA3ox from 38.00 to 40.02 kDa (Table 1). Using Plant-mPLoc software, we identified the subcellular localization of the HsGAox proteins. Analysis revealed that all members from HsGA20ox, HsGA2ox and HsGA3ox are localized in the cytoplasm. Interestingly, only one member, HsGA2ox2 protein, displayed localization in both the cytoplasm and nucleus (Table 1).

To determine the exact genomic positions of HsGAox genes, chromosome localization was carried out. A total of 27 HsGAox genes were widely distributed across eleven H. sosnowskyi chromosomes (Figure 1). HsGA20ox gene subfamily members were identified to be distributed in chromosomes 1, 2, 3, 7 and 8. The vast majority of gene loci are localized in the lower arm of the chromosomes; however, HsGA20ox7 and HsGA20ox8 are situated in the upper arm and in the middle of chromosome 3, respectively. Thirteen HsGA2ox genes are scattered across six chromosomes. Eight HsGA2ox1, 2, 5, 6, 9, 10, 11, 12 loci are located in the lower arm of the chromosomes 2, 4, 5, 6, 10 and 11. The rest of the genes localized in the upper arm of the chromosomes 6, 10, and 11. The H. sosnowskyi GA3ox subfamily genes are scattered among chromosomes 5, 9, 10 and 11. HsGA3ox1, 4 and 5 are located in the lower arm of chromosomes 11 and 5, respectively. The HsGA3ox3 locus was located in the end of the upper arm of chromosome 9. Interestingly, chromosomes 1-6 are considered longer than others, but HsGAox genes are distributed across all chromosomes and do not show a dominant correlation to chromosome length.

2.2. Phylogenetic Analysis of the HsGAox Genes

To explore the phylogenetic relationship of HsGAox with homologues from other taxa, we involved two model species, Arabidopsis thaliana and Oryza sativa, in the analysis. A maximum-likelihood (ML) phylogenetic tree was built based on alignments of the complete protein sequences from 27 Heracleum sosnowskyi, 16 A. thaliana, and 21 O. sativa found in NCBI and Phytozome databases. Four subfamilies (GA20ox, C₁₉-GA2ox, C₂₀-GA2ox and GA3ox) were identified. However, H. sosnowskyi proteins clustered into 3 different subfamilies: GA20ox, C₁₉-GA2ox and GA3ox. Interestingly, among thirteen HsGA2ox members, none belonged to the C₂₀-GA2ox subfamily. Moreover, HsGA20ox1 and -2, HsGA20ox5 and -6 appear to be duplicate proteins. These results allow us to make speculations on the functions of the HsGA20ox, HsGA2ox, and HsGA3ox according to the phylogenetic classification mentioned above (Figure 2).

2.3. Conserved Motif and Gene Structure Analysis of the HsGAox Genes

To support the phylogenetic analysis, we predicted conserved motifs of the HsGAox protein family using the Multiple Em for Motif Elicitation (MEME) online tool. 10 conserved motifs for HsGA20ox, C₁₉-HsGA2ox and HsGA3ox subfamily proteins were obtained and ranged from 21 to 50 amino acid residues (Figure 3B). The assessment of motifs 1–10 with InterPro Pfam revealed that motif 1 exhibits 2OG-Fe(II) oxygenase and motif 2 exhibits non-haem dioxygenase N-terminal domains. Except for HsGA3ox, all HsGA20oxs and HsGA2oxs were endowed with 8 motifs in the same pattern. In contrast, HsGA3oxs contained 7 motifs (HsGA3ox4 had 6 motifs), moreover, motif composition of the HsGA3oxs was distinct from HsGA20oxs and HsGA2oxs by the C-terminal end, where motif 6 or motif 8 was absent. In more detail, HsGA20oxs subfamily proteins are distinct from other HsGAoxs, containing a family specific motif 6 at N-terminal end of the sequence. Furthermore, motif 10 is unique for HsGA20ox and HsGA3ox subfamily members; in contrast, at the same location of HsGA2ox sequences, motif 5 is observed.

We analysed the intron and exon compositions, as well as dispositions in the gene structure of HsGAoxs (Figure 3C). Analysis revealed that the majority of HsGA20ox members had a consistent number of 3 exons and 2 introns, whereas HsGA20ox4 and -7 number of exons varied from 4 to 5, respectively. HsGA2ox subfamily genes contained 3-5 exons: 10 members (HsGA2ox1,3,5-12) had 3 exons, HsGA2ox4 – 4 exons and HsGA2ox2,13 – 5 exons. The HsGA3ox subfamily (HsGA3ox1-5) showed a relatively simple gene structure, consisting of 2 exons and 1 intron. In summary, the length of introns varied within three HsGAoxs gene subfamilies, contributing to the overall size differences among gene members in these subfamilies. To sum up, gene structure and motif analyses support the phylogenetic findings, indicating the conserved evolution of these gene subfamilies in all three species.

2.4. Analysis of HsGAoxs Expression in Response to GA₃ Within Heracleum sosnowskyi Developing Ovaries

To evaluate the exogenous GA₃ effect on the expression of GA biosynthesis enzyme genes, the Real-time quantitative PCR (qPCR) method was used to determine the quantity of gene transcripts in the terminal umbel ovaries, which were collected from the central and lateral parts at different stages of development. Specific primers for the putative HsGA20ox1, HsGA3ox1, and HsGA2ox1 genes were designed and used in the analysis. It was found that the HsGA3ox1 gene was most intensely expressed at 3 and 10 days after the start of the experiment, both in the central and lateral parts of the terminal umbel (Figure 4B). The quantity of HsGA20ox1 transcripts in the terminal umbel ovaries peaked 3 days after the start of the experiment, but the expression of this gene decreased after 10 days from the onset of flower spreading (Figure 4B). Due to the effect of GA₃, the relative abundance of HsGA20ox1 gene transcripts decreased by almost half on the third day of the experiment, reaching 8.9 % in the central part (p = 0.04) and 8.7 % in the lateral part (p = 0.85). After 10 days, the expression of this gene was even weaker in both the central (2.4 %) and lateral (2.2 %) parts. Analysis of HsGA3ox1 gene expression revealed that 3 days after the GA₃ treatment, the expression of this gene decreased from 41.5 % to 7.18 % (p = 0.8) in the central part of the umbel and from 45.3 % to 7.27 % (p = 0.72) in the lateral part. Ten days after the GA₃ treatment, the relative abundance of HsGA3ox1 gene transcripts significantly increased, reaching 235 % (p < 0.001) in the central part and 168 % (p = 0.0019) in the lateral part of the umbel.

2.5. Changes in Endogenous GA levels in Heracleum sosnowskyi Ovaries After Treatment with Exogenous GA₃

The dynamics of endogenous GAs were also analysed in another experiment in which we aimed to determine the distribution of metabolites of the C-13-hydroxylated GA biosynthesis pathway in the central and lateral parts of the terminal umbel after GA₃ treatment (Figure 4C). The obtained results revealed that the level of biologically active GA₁ was extremely low in both parts of the umbel at 0 and 3 days after application (DAA): 963 pg/g FW in the central part and 264 pg/g FW in the lateral part. However, at a later stage of ovary development (10 DAA), the hormone level in the central part of the umbel significantly increased to 16.3 ng/g FW (p = 0.011) and to 43.6 ng/g FW (p < 0.001) in the lateral part. Interestingly, the GA₁ content in the ovaries of the central and lateral parts differed significantly, by almost a twofold difference (p < 0.001). The level of the precursor GA₂₀ reached its peak in the central (72.47 ng/g FW) and lateral (46.64 ng/g FW) parts of the umbel during the 3 DAA period. However, during the 10 DAA period, the level of this metabolite significantly decreased to 4.37 ng/g FW (p < 0.001) and 2.87 ng/g FW (p = 0.005) in the ovaries of the central and lateral parts, respectively. The level of the GA₂₀ catabolite (GA₂₉) remained low throughout the experiment. The level of the GA₁ catabolic form (GA₈) increased in the central part (12.8 ng/g FW) and in the lateral part (7.3 ng/g FW) during the 3 DAA period, and then significantly decreased to 1.19 ng/g FW (p < 0.001) and 0.89 ng/g FW (p = 0.059) in each part of the umbel, respectively, during the later period (10 DAA). The highest levels of early hormone precursors GA₄₄ and GA₁₉ in the terminal ovaries were detected only at the 3 DAA developmental stage. Interestingly, in the early stage of ovary development, GA₄₄ was found four times more abundant in the central part of ovary tissues (32.2 ng/g FW) than in the lateral part (7.35 ng/g FW), and the difference between them was significant (p < 0.001). After spraying the terminal umbels with exogenous GA₃, the GA₁ content increased insignificantly in the ovaries of the central part from 0.963 to 5.15 ng/g FW (p = 0.95), and significantly in the lateral part from 0.264 to 20.9 ng/g FW (p < 0.001) after 3 DAA. 10 days after application with GA₃ solution, an uneven distribution of GA₁ hormone levels was observed depending on the umbel area: significantly decreased in the central part to 388 pg/g FW (p < 0.001) and significantly decreased in the lateral part to 26.86 ng/g FW (p < 0.01). Due to GA₃ treatment, the level of catabolite GA₂₉ increased dramatically in the ovaries of the lateral part of the umbel: 6.66 ng/g FW (p < 0.001) and 6.28 ng/g FW (p < 0.01) in the 3 DAA and 10 DAA periods, respectively.

2.6. GA₃ Impact on the Phenotype of Heracleum sosnowskyi Ovaries

Morphometric analysis of the samples revealed that the development of mericarps is not uniform in the central and lateral parts of the umbel: after 10 days, the mericarps in the central part of the umbel were 11 % shorter (U = 2000.5, p < 0.001), 10 % narrower (U = 2385, p < 0.01), and 9 % lighter (U = 2185.5, p < 0.001) than the mericarps from the lateral part of the umbel (Figure 4B). Three days after GA₃ 150 mg/l treatment, there was a statistically significant decrease in mericarp length (U = 1459, p < 0.001) and weight (U = 2004, p < 0.01) in the central part, and width (U = 2203.5, p < 0.001) in the lateral part. Ten days after application, there was a statistically significant decrease in mericarp length (21 %) (U = 1656, p < 0.001), width (30 %) (U = 1020, p < 0.001), and weight (59 %) (U = 1008.5, p < 0.001) in the central part; and a decrease in length (20 %) (U = 1337, p < 0.001), width (31 %) (U = 712.5, p < 0.001), and weight (41 %) (U = 703, p < 0.001) in the lateral part (Figure 4B).

3. Discussion

Gibberellins play an important role in the processes of fruit set, development, and ripening [48,49]. In horticulture, plants are often sprayed with bioactive GAs to obtain high-quality seedless fruits, but this can have a negative impact on the histological and morphological structure of the fruit [50,51,52,53]. The 2-ODD protein superfamily is the second-largest enzyme family in plants. Moreover, 2-ODDs have a crucial role in leading oxygenation or hydroxylation in various plant metabolic events [19]. Members of 2-ODDs have been identified in many model plant species: Arabidopsis thaliana, rice, cucumber, tomato and others [18,26,54,55,56,57,58,59]. However, information about 2-ODDs in invasive Heracleum sosnowskyi is very limited. A recent study by [47] made the genome of H. sosnowskyi available for genome-wide analysis. According to GAoxs protein sequences from A. thaliana and Oryza sativa [60], we identified and named 9 HsGA20ox1-9, 13 HsGA2ox1-13 and 5 HsGA3ox1-5 protein homologues from H. sosnowskyi (Table 1). All putative GAoxs contained 2OG-Fe(II) oxygenase (PF03171) and non-haem dioxygenase N-terminal (PF14226) domains, which are known to be specific for 2-ODDs [61,62]. In our study, we predicted the 3D protein structures of several HsGAoxs and identified that both previously named domains are located in the reaction crevice of the proteins (Figure 3A). This aligns with the previous models of suggested GAoxs in rice [63]. Interestingly, HsGA20ox1 and HsGA20ox2, seem to be duplicates, which have one amino acid substitution caused by single nucleotide polymorphism (Figure 3A). Gene duplication is generally seen as a major driver of genes developing different functions, a process steered by natural selection, which is important to consider in further HsGAox analysis [64,65]. Phylogenetic analysis suggested that HsGA20oxs belong to three distinct clades, similarly, HsGA3oxs seem to belong to two different paraphyletic groups. Interestingly, we have not found any members of C₂₀-GA2oxs (Figure 2). However, the functional characterization using C₂₀- and C₁₉-GA substrates of HsGA2oxs might improve our classification. Previous studies indicate that GA2oxs are composed of three groups: C₁₉-GA2ox-I, C₁₉-GA2ox-II, and C₂₀-GA2ox-I [11,58,66,67]. In our study, we observed that H. sosnowskyi C₁₉-2ox members cluster together with O. sativa and A. thaliana members into two monophyletic groups: C₁₉-GA2ox-I and C₁₉-GA2ox-II (Figure 2). These findings suggest that C₁₉-HsGA2oxs might have functional differences among the groups. Furthermore, it is known that different GAox genes display significant variations in their expression across different tissues [39,68,69,70,71]. Comparison of the phylogenetic tree with the gene structure and motif analysis, it is apparent that the most closely related members within subfamilies share similar gene structure and motif composition (Figure 2; Figure 3). Interestingly, HsGA20ox and HsGA2ox are found in greater numbers than HsGA3oxs, similarly to A. thaliana, O. sativa, Cucumis sativus, Salix matsudana [18,63,72,73]. To delve into the evolutionary relationships of GAox genes, we analysed motif composition and gene structure (Fig. 3). However, GA20ox and GA3ox subfamilies were distinct from GA2ox by the LPWKET motif, which is characteristic of GA20oxs and for some GA3oxs [39,63,72]. Whereas motif 9 was observed in all C₁₉-HsGA2oxs, it is distinct from the rest of the 2-ODD subfamilies [63]. Gene structure revealed that all HsGA20oxs contain three exons and two introns except for HsGA20ox4 and -7, which are similar to the structure of the GA20ox subfamily in A. thaliana [58]. HsGA3oxs have a simple gene structure of one intron and two exons, consistent with those in A. thaliana, C. sativus and O. sativa [58,60]. Our findings reveal that closely related members within GAox subfamilies possess similar structures and motifs, implying shared functions.

A more detailed analysis of terminal umbel ovaries in the early stages of fruit development revealed a differential distribution of GA and metabolites between the central and lateral regions. This was particularly reflected in a significant increase in the levels of metabolites GA₄₄ and GA₂₀ in the control samples of the central umbel region at 0 and 3 days after application (DAA) (Figure 4D). This inter-regional distribution may be due to the complex umbel architecture, as the flowers in the lateral part of the umbel open first [74]. In addition, the literature indicates that GA is involved in the sex determination of male flowers [75,76,77], with high levels of GA precursors and the hormone GA₄ accumulating in the floral parts of stamens [78]. As is known, male flowers in the umbels of some Apiaceae family representatives, such as Zizia aurea and Thaspium barbinode, are distributed centripetally [79]. However, the terminal umbel of Heracleum mantegazzianum is usually composed only of hermaphrodite flowers [74]. There is no literature data on the distribution of hermaphrodite and male flowers in the terminal umbel of H. sosnowskyi but based on our endogenous GA profiles (Figure 4D), it can be assumed that there should be more of them in the centre of the H. sosnowskyi umbel. Three days after exogenous GA₃ treatment, the accumulation of the metabolite GA₂₀ was significantly reduced in the ovaries of both parts of the umbel, but the level of the hormone GA₁ in the ovary tissues increased (Figure 4D). Ten days after application, the synthesis of the hormone GA₁ in the ovary tissues of the central part of the umbel was completely suppressed, but only partially in the lateral part of the umbel. The inhibitory effect of exogenous GA₃ on GA biosynthesis can be identified by a significant decrease in the precursor GA₂₀ and a marked increase in the amount of the catabolite GA₂₉ in the ovaries of both parts of the umbel. However, it remains unclear why exogenous GA₃ forms a GA₁ gradient in different parts of the umbel (Figure 4D). It is known that an increase in the level of biologically active GA forms in the early stages of fruit development leads to the development of parthenocarpic fruit set [38,49,51,80,81]. It should be noted that in this experiment, flower pollination was not controlled, so the GA₁ content at 10 DAA could have increased due to natural fertilization (Figure 4D).

Before analysing the abundance of GA biosynthesis gene (HsGA20ox1, HsGA3ox1, and HsGA2ox1) transcripts, these genes were cloned from H. sosnowskyi umbel tissues collected 10 days after application (DAA) and specific primers were designed for them. It was observed that in the early stage of fruit development (3 DAA), HsGA20ox1 was more highly expressed, but its expression decreased in later stages (10 DAA). On the other hand, the expression of the metabolic HsGA3ox1 gene was much more intense, with no differences detected between umbel regions. Other studies have indicated that the increase in GA20ox gene expression and the decrease in GA2ox are essential factors in regulating GA biosynthesis during fruit set [48,51,82,83]. It should also be noted that in other plant species (e.g., white clover, tomato, cultivated rice), more than one GA oxidase gene is involved in fruit set control [28,49,78,80,84]. However, our results suggest that intense fruit development is dependent on endogenous GA, which is supported by the obtained levels of endogenous GA (Figure 4E). Additionally, the results of phenotypic analysis of ovaries (Figure 4B) correspond to the peak of gibberellin during the cell expansion phase observed in the model plant tomato [85,86]. However, the effect of exogenous GA₃ stimulated HsGA3ox1 expression in the central and lateral parts of the umbel ovaries, while the inhibition of HsGA2ox1 gene expression, as speculated, was not detected (Figure 4E). Furthermore, the obtained expression profile of the HsGA20ox1 gene at 10 DAA, both in the control and after GA₃ treatment, suggests that fruit set had already occurred, as seen in some species such as pear, where after a similar treatment with exogenous biologically active GA, the expression of genes PbGA20ox1, PbGA20ox2, and PbGA20ox3 is suppressed [39]. Our results align with the suppressed expression of paralogous genes of the PbGA20ox (PbGA20ox1, PbGA20ox3) in pollinated pear ovaries [39].

4. Materials and Methods

4.1. Research Object and Growth Conditions

An invasive habitat of Heracleum sosnowskyi was selected for field trials, and an experimental field was established. The test site was located in an anthropogenic H. sosnowskyi habitat with an area of 6.48 ha in an abandoned area in Vilnius, Lithuania (54.739749°N, 25.258871°E WGS80). The experiment field was fenced with a “stop” line and marked with signs informing about the plant research being conducted. The studies were conducted on H. sosnowskyi plants during the flowering and fruit formation stages from 2021. Samples were collected from 3 plants. The terminal umbel was divided into central and lateral parts following the guidelines [87].

4.2. Application of Plants with GA₃

Terminal umbels of Heracleum sosnowskyi were sprayed with the phytohormone gibberellic acid (GA₃) (SERVA, Heidelberg, Germany) at a 150 mg/l concentration. Certain order umbels were sprayed twice. Distilled water was used as a control spray. Manual sprayers “Venus” (KWAZAR, Budy-Grzybek, Poland) were used for spraying. Inflorescences were sprayed from a distance of 15-20 cm, evenly applying 18.5 ml onto each umbel.

4.3. Sample Harvesting

Samples for morphological and anatomical analysis of Heracleum sosnowskyi flowers and mericarps were collected in August-September from terminal umbels. To study ovary development dynamics in terminal umbels, ovary samples were collected from the central and lateral parts of the terminal umbel at 0, 3 and 10 days after application (at the beginning of flowering). For biochemical and molecular studies, flowers of different developmental stages were collected, weighed, flash-frozen in liquid nitrogen and stored in a low-temperature freezer (Skadi Green line, ES) at -80 °C until analysis.

4.4. Bioinformatic Analysis of Heracleum sosnowskyi GAoxs Subfamilies

In order to perform a comprehensive phylogenomic analysis of the GAox subfamilies in Heracleum sosnowskyi, a prerequisite was the in silico extraction of all potential GAox homolog sequences from publicly available whole-genome sequence data. We extracted 16 Arabidopsis thaliana and 21 rice GAox sequences from NCBI Phytozome and Rice Genome Annotation Project databases. These sequences were grouped according to protein subfamilies and used as query sequences against H. sosnowskyi genome (NCBI: PRJNA928505) in BLASTp [88]. Overall, 32 amino acid homologues were obtained and revised. Candidates were accepted if they shared at least 40% identity and had an expected threshold e-value ≤1.0 × e^-10. Subsequently, HsGAoxs protein candidates were submitted to InterPro Pfam database to identify the highly conserved and characteristic 2OG-FeII_Oxy (PF03171) and DIOX_N (PF14226) domains of the 2-ODDs superfamily [61,62]. Conserved motif analysis of HsGAoxs was performed using Multiple Em for Motif Elicitation (MEME) 5.5.7 software with default parameters [89]. The exon-intron gene structure of HsGAoxs genes was analysed using GSDS 2.0 online software [90]. Based on the CDS region range found in the NCBI database annotation of each HsGAox protein, we visualized gene locations on H. sosnowskyi chromosomes using MapGene2Chromosome v2 online tool [91]. In a subsequent analysis, the molecular weights of all HsGAox proteins were calculated using the Protein Molecular Weight tool from the bioinformatics.org website. Additionally, the Plant-mPLoc 2.0 software was applied to predict the subcellular localizations of HsGAox proteins [92]. Protein structures of HsGAoxs were predicted using Alphafold 3 [93], and the predicted model structure quality was assessed by the VoroMQA online tool using default parameters [94].

4.5. Multiple Sequence Alignment and Phylogenetic Analysis

Multiple sequence alignment of the GAoxs protein sequences from Heracleum sosnowskyi, Arabidopsis thaliana and Oryza sativa was conducted using MEGA X software through the application of MUSCLE (Multiple Sequence Comparison by Log-Expectation) algorithm with default parameters [95].

Optimal model fit for aligned protein sequences was generated, and a phylogenetic tree was constructed using Maximum Likelihood (ML) and Le-Gascuel model by MEGA 12 software [96]. The robustness of the phylogenetic inference was evaluated using 500 bootstrap replicates.

4.6. Extraction of Endogenous GAs

Endogenous GAs were extracted from fresh mericarp tissues that had been frozen in liquid nitrogen (N₂) according to the specified methodology [22]. GA extraction consisted of three stages. The first stage entailed sample extraction with 80 % methanol (purity ≥ 99.95 %, Carl Roth GmbH + Co. KG, Karlsruhe, Germany), introduction of a deuterated 17,17-d₂-GA standard mixture (internal standard) into the sample (purity ≥ 90 %, OlChemim, Olomouc, Czech Republic), and liquid-liquid extraction. The second stage involved solid-phase extraction using anion exchange “Bond Elut DEA” (Agilent, Santa Clara, USA) and reverse-phase “Sep-pak C18” (Waters, Wexford, Ireland) columns. The third stage was derivatization, during which the extracted samples were methylated with freshly prepared diazomethane according to De Boer and Backer [97] and trimethylsilylated with N-methyl-N-(trimethylsilyl)trifluoroacetamide (MSTFA) (purity ≥ 95 %, Macherey-Nagel, Düren, Germany). After extraction, the samples were evaporated under a stream of N₂ gas and stored at –20 °C until analysis by GC-MS method.

4.7. Endogenous GA Analysis by Gas Chromatography-Mass Spectrometry

Qualitative analysis of endogenous Gas was performed by GC-MS using a “GC/MS-Q2010 PLUS” (Shimadzu, Kyoto, Japan) gas chromatography system coupled with single quadrupole “GC-MS-QP2010 ULTRA” (Shimadzu, Kyoto, Japan) mass spectrometer. 1 μl of the sample was manually injected into the injector of the gas chromatography system coupled to a non-polar “BPX5” column (30 m × 0.25 mm × 0.25 μm) (SGE Analytical Science, Kiln Farm Milton Keynes, United Kingdom). Helium was used as the carrier gas. The split ratio was 30:1 for sample introduction into the column. The initial column temperature was 60 °C and was immediately raised at a rate of 45 °C/min to 220 °C, followed by a gradual increase at 4 °C/min intervals up to 300 °C. The interface temperature was 240 °C. The analysis time for one sample was 25 min. Data were recorded in SIM (selected ion monitoring) mode after 5 min from injection. Target ions for GA qualitative analysis were used for the corresponding GAs and d₂-GA: 300 and 302 (GA₁₂), 239 and 241 (GA₁₅), 314 and 316 (GA₂₄), 270 and 272 (GA₉), 284 and 286 (GA₄), 506 and 508 (GA₃₄), 284 and 286 (GA₅₁), 207 and 209 (GA₅₃), 207 and 209 (GA₄₄), 374 and 376 (GA₁₉), 418 and 420 (GA₂₀), 506 and 508 (GA₁), 594 and 596 (GA₈), 506 and 508 (GA₂₉). Compound identification was based on their retention times and the presence of additional characteristic ions in the obtained mass spectra, by comparing them with the known mass spectra of labelled and unlabelled GA methyl esters [98]. Endogenous GA levels were calculated based on peak areas, with corrections for the abundance of natural isotopes in the samples and the presence of unlabelled GAs in the internal GA standards [99].

4.8. Identification of GA Biosynthesis Gene Fragments

To amplify putative Heracleum sosnowskyi GA 20-oxidase (GA20ox), GA 3-oxidase (GA3ox) and GA 2-oxidase (GA2ox) gene fragments, a PCR based method was applied [22,100]. Based on the conserved amino acid sequences and nucleotide sequences of the GA20ox and GA3ox enzyme subfamilies in wild carrot (Daucus carota) (Apiaceae), degenerate primers were designed to amplify these genes. For GA2ox, a degenerated primer pair from cucumber (Cucumis sativus) (Cucurbitacea) was used. cDNA molecules were synthesized from frozen fresh mericarp samples of H. sosnowskyi. The resulting cDNA was used as a template in a PCR reaction, which was performed using the “Phusion High-Fidelity PCR Master Mix” (Thermo Scientific, Vilnius, Lithuania) kit, according to the manufacturer’s protocol. PCR products were analysed by electrophoresis on a 2 % agarose gel, stained with “Midori Green Advance” (Nippon Genetics, Duren, Germany) dye and visualized with a “FastGeneR GelPic LED Box” device (Nippon Genetics, Duren, Germany). The expected size of the DcGA20ox PCR product was 813 base pairs (bp), the expected size of the DcGA3ox product was 497 bp, and the expected size of the CsGA2ox product was 300 bp. PCR product bands were cut from the agarose gel and samples were purified using the “Gene JET^TM PCR Purification Kit” according to the manufacturer’s instructions (Thermo Scientific, Vilnius, Lithuania). Samples were prepared for sequencing according to the recommendations of “Microsynth SEQLAB” (Germany). Sequencing results allowed for the identification of regions in the H. sosnowskyi genome based on which specific primer pairs were designed. Using “NCBI PrimerBLAST” and “Primer3” software, forward (FW) and reverse (REV) primers were generated for the putative H. sosnowskyi HsGA20ox1, HsGA3ox1, HsGA2ox1 genes (Table S1).

4.9. Real-Time Quantitative PCR

Real-time quantitative PCR (qPCR) was performed to assess the transcript levels of late gibberellin biosynthesis enzyme (GA20ox, GA3ox, GA2ox) genes in Heracleum sosnowskyi ovary samples. Specific primers were used in the qPCR analysis to evaluate the expression of the genes encoding these proteins (Table S1). 1 µl of cDNA diluted 5 times with double distilled water was used as the template for the qPCR reaction. The total reaction mixture volume was 15 µl. Reactions were performed in an “AZURE CIELO^TM Real-time PCR System” (Azure Biosystems, Dublin, USA) thermocycler using the two-step cycle protocol of the “Maxima SYBR Green noROX qPCR Master Mix (2X)” (Thermo Scientific, Vilnius, Lithuania) kit. Cycling conditions consisted of one cycle at 95 °C for 10 min and 40 cycles at 95 °C for 15 s followed by 60 °C for 1 min. The cucumber (Cucumis sativus) β-actin gene (AB010922) was used as the internal standard. The relative abundance of each gene was calculated as the average of two technical replicates using the “Real-time PCR Miner” bioinformatics software and expressed as a percentage [101]. Technical replicates with a Ct difference ≥ 0.5 from the mean were discarded. All experiments were repeated twice with three biological replicates.

4.10. Statistical Analysis

Using descriptive statistical methods, the obtained results were presented as arithmetic means with standard errors (mean ± SE). The normal distribution of data was assessed using the Shapiro-Wilk test. If the assumption of data normality was met, parametric methods were employed, and the homoscedasticity of the data was assessed using Levene’s test with Brown-Forsythe correction. If the data were not normally distributed, statistically significant differences between group medians were assessed using the Kruskal-Wallis H-test. For pairwise comparisons between experimental groups, the Mann-Whitney U test or Dunn’s z test was applied. Statistically significant differences were indicated when p < 0.05. Data statistical analysis and graphical representation were performed using “PAST 4.16” and “R” version 4.0.2 in “RStudio” version 1.3.1093 with “ggplot2” package [102,103].

5. Conclusions

A total of 27 HsGAox genes were widely distributed across eleven Heracleum sosnowskyi chromosomes. Phylogenetic classification allows us to make speculations about the functions of the HsGA20ox HsGA2ox, and HsGA3ox. The analysed H. sosnowskyi proteins clustered into three different subfamilies. Among thirteen HsGA2ox members, none belonged to the C₂₀-GA2ox subfamily, contrary to the two subfamilies of GA2oxs in other plant species. Gene structure and motif analyses support the phylogenetic findings, indicating the conserved evolution of these gene subfamilies across three species. Intense fruit development is dependent on endogenous GA, as supported by the obtained levels of endogenous GA. The phenotypic analysis of ovaries corresponds to the peak of gibberellin during the cell expansion phase. However, the effect of exogenous GA₃ stimulated HsGA3ox1 expression in the central and lateral parts of the umbel ovaries. Overall, these results confirm that GA homeostasis is a tightly regulated, complex system including multiple HsGAox genes. It maintains a stable hormone level in the early stages of fruit development. A high or low level of the hormone in a particular umbel part can disrupt the development of the fruit. These results open opportunities to further analyse the role of GAs in H. sosnowskyi fruit-set mechanism and develop invasion control strategies.

In this study, 27 HsGAox genes were identified in H. sosnowskyi, distributed across eleven chromosomes. Phylogenetic analysis revealed three distinct HsGAox subfamilies, indicating the functional roles of HsGA20ox, HsGA2ox, and HsGA3ox. The HsGA2ox group, which is distinct from other plant species and lacks C₂₀-GA2ox. The conserved motif and gene structure patterns across species support the phylogenetic classification. It was found that endogenous gibberellin (GA) levels are crucial for robust fruit development, and ovary phenotypes coincide with the highest GA levels during cell expansion. External GA₃ application upregulated HsGA3ox1 in the ovary regions of the umbels. Complex GA homeostasis, regulated by multiple HsGAox genes, is crucial for stable early fruit development; GA imbalances in the umbel parts negatively affect fruit formation. These findings provide a basis for investigating the role of GA in H. sosnowskyi fruit formation and developing strategies to control the invasive spread.