Genome-Wide Identification and Evolution Analysis of the Gibberellin Oxidase Gene Family in Six Gramineae Crops

The plant hormones gibberellins (GAs) regulate plant growth and development and are closely related to the yield of cash crops. The GA oxidases (GAoxs), including the GA2ox, GA3ox, and GA20ox subfamilies, play pivotal roles in GAs’ biosynthesis and metabolism, but their classification and evolutionary pattern in Gramineae crops remain unclear. We thus conducted a comparative genomic study of GAox genes in six Gramineae representative crops, namely, Setaria italica (Si), Zea mays (Zm), Sorghum bicolor (Sb), Hordeum vulgare (Hv), Brachypodium distachyon (Bd), and Oryza sativa (Os). A total of 105 GAox genes were identified in these six crop genomes, belonging to the C19-GA2ox, C20-GA2ox, GA3ox, and GA20ox subfamilies. Based on orthogroup (OG) analysis, GAox genes were divided into nine OGs and the number of GAox genes in each of the OGs was similar among all tested crops, which indicated that GAox genes may have completed their family differentiations before the species differentiations of the tested species. The motif composition of GAox proteins showed that motifs 1, 2, 4, and 5, forming the 2OG-FeII_Oxy domain, were conserved in all identified GAox protein sequences, while motifs 11, 14, and 15 existed specifically in the GA20ox, C19-GA2ox, and C20-GA2ox protein sequences. Subsequently, the results of gene duplication events suggested that GAox genes mainly expanded in the form of WGD/SD and underwent purification selection and that maize had more GAox genes than other species due to its recent duplication events. The cis-acting elements analysis indicated that GAox genes may respond to growth and development, stress, hormones, and light signals. Moreover, the expression profiles of rice and maize showed that GAox genes were predominantly expressed in the panicles of the above two plants and the expression of several GAox genes was significantly induced by salt or cold stresses. In conclusion, our results provided further insight into GAox genes’ evolutionary differences among six representative Gramineae and highlighted GAox genes that may play a role in abiotic stress.


Introduction
Gibberellins (GAs) form a large family of diterpenoid compounds and are involved in plants' development and growth in various aspects, including seed germination [1], stem elongation [2], leaf expansion [3], flowering [4,5], root growth, and fruit development [6][7][8]. The first discovery of gibberellin was in the pathogenic fungus Gibberella fujikuroi, and the fungus-infected plants showed an excessive elongation phenotype; this phenomenon is known as 'foolish-seedling' [9]. To date, more than 130 GAs have been discovered in plants, bacteria, and fungi, of which only a few GAs harbor biological activity, namely, GA 1 , GA 3 , Sorghum bicolor (Sb)) at the genome-wide scale. We systematically analyzed the phylogenetic relationships, gene structures, motif compositions, orthogroups (OGs), chromosomal locations, duplication events, selective forces, and cis-elements in the promoters of GAox genes. In addition, the RNA-seq data were used to analyze GAox genes' tissue expression patterns and their responses to cold or salt stress in rice and maize. Our results provided a basis for further study on the gene expansion, evolutionary patterns, and functional diversity of GAox gene families in Gramineae.

Plant Materials and Treatments
Nipponbare (Oryza. Sativa) rice seeds were cultivated with Youshida medium in the 26 • C greenhouse of Wuhan University, the light/dark photoperiod was 16/8 h, and relative humidity was 60% [34]. At three-leaf stage, some of the seedlings were transferred to Youshida medium with 125 mmol/L NaCl for salt treatment. The growth temperature of the other seedlings was changed to 4 • C for cold treatment. The roots were collected at 0, 3, and 24h after stress treatments, and three biological replicates for each treatment were conducted.

Phylogenetic Tree Construction of GAox Genes
For phylogenetic analysis, the final candidate GAox proteins of each species were aligned by MUSCLE v3.8.1551 [38], trimAl v1.4 was used to filter gap and nonconservative sites [39], and then the phylogeny trees were generated by iqtree v2.0.3 with 1000 bootstrap replicates [40]. To ensure the accuracy of phylogenetic tree, the optimal model identified by iqtree was adopted to rebuild the tree by RaxML v0.9.0 [41]. Then, the Evolview online site (https://www.evolgenius.info/evolview/, accessed on 10 January 2021) was used to visualize the phylogenetic tree [42].

Homology Analysis of GAox Genes
The lineal ortholog gene groups of GAoxs in six plants were analyzed by OrthoFinder v2.4.0 software [43]. An all-vs.-all BlastP search of all OsGAox protein sequences was performed by BlastP through diamond v0.9.24.125 software with default parameters as the input file for OrthoFinder software [44]. Finally, the candidate GAox genes were named according to their homology in rice [32]. In addition, DnaSP 5.0 (http://www.ub.edu/ dnasp/, accessed on 25 January 2021) software was used to calculate the value of Tajima's D for each OG [45].
To explore the expansion modes and differentiation times of GAox genes in six Gramineae, the duplication events were analyzed by MCScanX software [46]. First, we used Diamond software to perform BlastP among the GAox genes in each species [44], and then the 'duplicate_gene_classifier' script in the MCScanX software was adopted to identify the paralogous genes with the default parameters. The nonsynonymous/synonymous (Ka/Ks) substitution rates of the identified gene pairs were calculated by DnaSP 5.0 [45], and the divergence times were estimated by T = Ks/(2 × 9.1 × 10 −9 ) × 10 −6 million years ago (Mya) [47].

Distribution and Structure of GAox Genes
The chromosomal positions of GAox genes were obtained through GFF3 files of six species and visualized through Tbtools v1.089 [48]. All putative GAox protein sequences were submitted to the MEME Suite 5.3.3 (http://meme-suite.org/tools/meme, accessed on 12 February 2021) with the maximum number of motif sets at 15 and the optional width of motifs from 6 to 150 to check conservative motifs [49]. Both gene structure and conservative motifs were visualized by Tbtools [48].

RNA Extraction and qRT-PCR Analysis
Total RNA from the samples were extracted by TRIzol reagent. Purified total RNA were reversed to first-strand cDNA using HiScript III 1st Strand cDNA Synthesis Kit (Vazyme, Shanghai, China). The qTR-PCR analysis was conducted by 2 × SYBR Green qPCR Master Mix (US Everbright ® Inc., Suzhou, China) according to the manufacturer's instructions. The qRT-PCR reactions were detected with CFX96 TouchTM Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The actin gene (UBI) was used as an internal reference control, and the primers in this experiment were designed by Primer Premier 5.0 (Table S1). The relative expression level was calculated based on three biological replicates using 2 − CT method [35].

Identification and OG Analysis of GAoxs in Six Gramineae
A total of 105 nonredundant GAox proteins were identified in six Gramineae crops: O. sativa (17), Z. mays (20), H. vulgare (17), S. bicolor (17), S. italica (15), and B. distachyon (19) (Table S2, Figure 1). The number of GAox genes was similar in the six species mentioned above, with the largest number of ZmGAoxs and the smallest number of SiGAoxs. In order to explore the evolutionary model of GAox genes among six Gramineae, both iqtree and RaxML were used to construct a phylogenetic tree based on an edited multiple-sequence alignment (MSA) file (Supplementary File S1-S3), and the results of the two software were consistent, which proved the reliability of the phylogenetic tree ( Figure S1). GAox genes could be clearly divided into four groups, as the GA20ox, GA3ox, C 19 -GA2ox, and C 20 -GA2ox subfamilies ( Figure 1A) [54]. Subsequently, we identified nine orthologous groups using OrthoFinder (Table S3), and GA3ox genes corresponded to OG-7, GA20ox genes were classified into OG-8,9, and the remaining six OGs all belonged to GA2ox genes ( Figure 1A,B). In each OG group, the number of orthologous genes among the tested species was nearly equal. For example, each tested species contained two GAox genes belonging to OG-4 or OG-7 ( Figure 1B). These results strongly indicated that the evolution of GAox families occurred before the separation of the major grass families.  Figure 1B), and this result emphasized that all OGs were subjected to varying degrees of purification selection. Among them, GA3ox genes (OG-7) had the highest Tajima's D-value, and we speculated that owing to the function of GA3ox genes being indispensable to plants and that there were only two members in each species, it faced strong selection pressure. This part of the results implied that GA3ox genes were relatively conservative, while in GA2ox and GA20ox genes, gene loss/gain may have occurred in the course of the evolution of six Gramineae after the separation of the major Gramineae families.

Chromosome Localization, Gene Duplication Events, and Selective Forces Analysis
In order to explore the expansion patterns of paralogous genes, we separately executed the gene positions and duplication events of GAox genes in six species. GAoxs were unevenly distributed on chromosomes (Chrs) (Figure 2). For example, five and six GAox genes were located on chromosomes 1 and 5 in rice, respectively. However, the other chromosomes only harbored 0-2 genes (Figure 2A). Interestingly, we noticed that GA20ox-1 genes from six species were commonly located on the edge of the chromosome, and the positions of GA2ox-4 and GA2ox-8 were in proximity in five plants except for Si. Thus, we speculated that the physical positions of some GAox genes were conservative during the Gramineae plant evolution process. Gene duplication events were some of the main drivers of genomics evolution [55]. Based on our results, gene duplication events were found in all tested species (Table 1) and also detected in the GA2ox, GA3ox, and GA20ox subfamilies, respectively. Their expansions were mainly through whole-genome duplication or segmental duplication events (WGD/SD), while tandem duplications (TD) only existed in Hv, Bd, and Si ( Figure 2B,D,E). Notably, some paralogous pairs were simultaneously identified in several plants, such as GA2ox-3. GA2ox-4 as well as GA2ox-6 gene pairs existed in Os, Hv, Zm, Si, and Sb, and GA2ox-9 gene pairs occurred in Os, Zm, Si, Bd, and Sb, implying that these paralogous genes had been duplicated prior to the ancestral divergence ( Figure 2).  ). Maize harbored the most GAox genes, which was due to the recent doubling of several pairs of genes, such as ZmGA2ox-7a-ZmGA2ox-7b, ZmGA2ox-3a-ZmGA2ox-3b, and ZmGA20ox-2a-ZmGA20ox-2b.

Gene Structure and Conserved Motif Analysis
The GAox gene structures were visualized based on the GFF3 files of six species ( Figure 3B). The lengths of coding sequences and the numbers of exons in all tested GAox genes were various, of which most GAox genes consisted of 1 to 3 exons. However, even in the same OGs, some genes contained long introns, while others had only one exon, such as OG-1. Therefore, it was speculated that in GAox genes, sequence loss and exon fusion might occur during evolution. We executed the MEME online tool and searched 15 motifs in 105 GAox protein sequences ( Figure S2, Table S4). The composition and arrangement of motifs in four subgroups were similar ( Figure 3A). The conserved domain of 2OG-FeII_Oxy consisted of motifs 1, 2, 4 and 5, which were identified in all GAox protein sequences ( Figure 3A). It was worth noting that the tails of GAox protein sequences contained specific motifs that could distinguish different subfamilies. For example, motif 11, motif 14 and motif 15 existed only in GA20ox, C 19 -GA2ox, and C 20 -GA2ox, separately. Furthermore, motif 7 existed specifically in the GA3ox and GA20ox genes. GA3ox and GA20ox proteins were responsible for the synthesis of gibberellin, while GA2ox proteins degraded gibberellin [24], so motif 7 may have been related to the substrate recognition of GAox gene families. GA2ox proteins were classified as C 20 -GA2ox or C 19 -GA2ox genes based on their preference for C 20 -GA or C 19 -GA substrates [20], and their differences in conserved domains were reflected in motifs 8, 14, and 15. In general, our results revealed that GAox families formed unique conserved motifs in the process of differentiation.

Cis-Elements Analysis of GAoxs
We identified a total of 31 cis-elements distributed unevenly on the promoter region which were involved in growth and development, stress response, light response, and phytohormone ( Figure 4A). Growth-and development-related cis-elements accounted for 8% ( Figure 4B), involving endosperm development, circadian rhythm, meristem, and metabolism. Thirteen percent of the cis-elements were relevant in the stress response, which was associated with drought, low temperature, and anoxic specific inducibility. A total of 16.3% of cis-elements were predicted to respond to a light signal, such as G-box, GATA-motif, MRE elements, and so on. Moreover, phytohormone-related cis-elements (43%) included auxin responsiveness (AuxRR-core), gibberellin responsiveness (TATC-box, P-box, and GARE-motif), and MeJA-responsiveness (CGTCA-motif and TGACG-motif). Previous studies had also successively proved the various functions of GAox genes [56,57]. These results indicated that GAox genes may execute multiple biological functions and play a key role in plant growth, development, and stress.

Tissue-Specific Expression of GAox Genes
We retrieved the expression data of rice and maize from the public database and our previous study ( Figure 5A,B; Table S5). GAox genes were widely expressed in rice or maize tissues. In rice, OsGA2ox-10 and OsGA20ox-1 were highly expressed in most tissues ( Figure 5A). On the contrary, the expression levels of OsGA2ox-2, OsGA2ox-4, OsGA2ox-5, OsGA2ox-6, and OsGA2ox-8 were low in most tissues ( Figure 5A). The majority of OsGAox genes were predominantly expressed in rice panicles, indicating that OsGAox genes were vital in panicle development, which was consistent with previous research results [58]. Moreover, the expression profiles of paralogous gene pairs were various ( Figure 5A). For example, OsGA2ox-10 was expressed highly in most tissues, while its paralogous genes OsGA2ox-7 were predominantly expressed in flowers and panicles, which revealed that these gene pairs underwent functional differentiation after duplication events. There were many similarities between the GAox genes' expression patterns in rice and maize. In maize, ZmGAox genes were predominantly expressed in ear primordium, which was consistent with that in rice. This result showed that GAox genes played a key role in panicle development of both maize and rice. Interestingly, GA2ox-10 and GA20ox-1 were clustered together in rice or maize tissue expression profiles, and two pairs of genes showed similar expression patterns, indicating that their expression in Gramineae seemed to be very conservative and played similar biological functions.

GAox Genes Responded to Cold or Salt Stress
In plants, GAox genes participated in plenty of physiological processes, especially in resisting abiotic stress, which played a significant role. Many researchers have investigated the performance of GAox genes under abiotic stress [20,59]. Under salt stress, there was little difference in the expression patterns of GAoxs between Chao2R (indica) and RPYgeng (japonica). For example, OsGA20ox-2 and OsGA20ox-3 were up-regulated and OsGA2ox-3 and OsGA2ox-8 were down-regulated in both varieties ( Figure 5C, Table S6). For cold stress, we obtained the transcriptome data of GAox families from TNG67 (japonica) under low temperature for 3 h and 24 h ( Figure 5C). Most members of OsGAox genes were sensitive to cold stress, and the expression of OsGA2ox-3 altered significantly. Moreover, OsGA3ox-1 responded differently to cold stress in various tissues, which was down-regulated in roots but up-regulated in shoots.
In our study, ZmGAox genes mostly responded to salt or cold stress by down-regulating their expression, while ZmGA2ox-10 and ZmGA2ox-6 were up-regulated under both stresses ( Figure 5D). There were several ZmGAox genes which showed the opposite expression pattern under salt or cold stress, such as ZmGA2ox-7a, ZmGA2ox-9, ZmGA2ox-11, and ZmGA20ox-1. It was worth noting that the great differences appeared in the expression patterns of homologous gene pairs in rice and maize. The exception was that ZmGA2ox-1a and ZmGA2ox-1b clustered together in the maize stress expression profile, and they may have been functionally redundant in response to stress. In order to verify the correctness of the RNA-seq data, we tested the expression levels of GAoxs in the roots of japonica rice (Nipponbare) treated with cold and salt for 3h and 24h, respectively. Five representative OsGAox genes that were significantly differentially expressed under stress were selected to perform qRT-PCR, and the test results showed that these genes were significantly induced by stress ( Figure 5E). This work was consistent with the analysis of the cis-elements of promoters. Most GAoxs were induced by cold or salt stress, and some GAox genes may have been functionally redundant.

Discussion
Bioactive gibberellins (GAs) are diterpene plant hormones that are biosynthesized through complex pathways and control diverse aspects of growth and development [10]. In addition, in the 1960s, the regulation of GAs on plant height had a significant influence on agriculture, and the cultivation of semi-dwarf wheat and rice contributed to a substantial increase in food production, which was called the Green Revolution [60]. Owing to the great significance of gibberellin oxidase in plant growth and development, we herein systematically identified 105 GAox genes in six Gramineous plant genomes. The number of GAox genes in each tested species was similar. The smallest genome (genome size: 490 Mbp) of S. italica contributed the fewest number of GAox genes. Although maize had the most GAox genes among six plants, it is obvious that ZmGAox genes did not double with their special genome-wide duplication [61]. Moreover, maize had the most GAox genes, and gene duplication events showed that they had experienced doubling recently, which meant that their functions may have been redundant, and our expression profile results also provided relevant proof, such as that ZmGA2ox-1a and ZmGA2ox-1b showed similar expressions in stress response.
The final stage of bioactive GAs synthesis is catalyzed by two GAox families, called the GA20ox and GA3ox genes [12,24]. Moreover, in the pathway of GA degradation, biologically active GAs or their immediate precursors are inactivated by two GA2ox families, including the C 19 -GA2ox and C 20 -GA2ox genes [10]. The above four GAox families all belong to the 2-ODDs superfamily and share an extremely conservative domain [62,63]. However, three families perform completely different functions [64,65]. In order to explore the sequence differences between them, motif identification was performed and revealed that some common motifs were shared by the four GAox subfamilies, while some unique motifs existed only in specific families. In our results, the GAox families (GA3ox and GA20ox genes) related to GA biosynthesis specifically contained motif 7 compared to the families (GA2ox genes) responsible for degradation, and motif 11 could further distinguish the GA3ox and GA20ox genes. The difference in motif composition between the C 19 -GA2ox and C 20 -GA2ox families depended on the existence of motif 14 ( Figure 3A). In summary, the four GAox subfamilies all harbored unique motifs that could distinguish each other, which may have been the result of the differentiation of the 2-ODDs superfamily. Whether these motifs are related to substrate binding is worthy of further exploration.
In the previous study, Han and Zhu identified eight GA20oxs, while we only identified four OsGA20ox genes in our Blastp results by using Arabidopsis' GAox genes [21]. We speculated that the previous incomplete genome, various parameters, or softs led to the different results. Liu et al. conducted functional analysis of OsGA20ox7, there was no significant change of GA content in the leaves of two OsGA20ox7 mutants generated by CRISPR/Cas9, and they speculated the role for OsGA20ox7 related to salicylic acid (SA) metabolism [66]. Another study found that OsGA20ox-5 and OsGA20ox-8 clustered in the CsGA7ox subfamily, which was a type of multifunctional enzyme with 7-oxidase and 3β-and 15 α-hydroxylase activity. This result indicated that OsGA20ox5 and OsGA20ox8 were homologs of CsGA7ox1 and CsGA7ox2 [67]. Furthermore, our conservative motif analysis also confirmed that OsGA20ox1-4 were significantly different from OsGA20ox5-8 ( Figure S3). In general, OsGA20ox5-8 may not belong to the GA20ox subfamily, and the functional study of OsGA20ox5-8 will further confirm it.
Here, we confirmed that the evolution of GAox families took place before their ancestors divided from several aspects. First, orthogroups analysis classified GAox genes into nine OGs, and the number of GAox genes in each OG was similar in each tested species ( Figure 1B), which indicated that these OGs were generated prior to six species' divergence. In addition, the gene duplication event is one of the driving factors of species evolution [55]. We identified 33 pairs of homologous genes with duplication events, and the tested species commonly shared some of gene pairs (Table 1), showing that the duplication of these gene pairs also occurred before the differentiation of six plants. Furthermore, the GA20ox-1 genes of all six species were located in the tails of chromosomes (Figure 2), indicating that they may have been located in the tails of chromosomes before the differentiation of the tested species. Thus, based on the above evidence, we speculated that the evolution of GAox families occurred before the ancestral separation.
The results of cis-acting elements analysis showed GAox genes may respond to a variety of biological processes, including growth and development, stress response, light response, and phytohormone. Tissue expression profiles showed that GAox genes were highly expressed in both rice and maize, which indicated GAox genes may regulate plant panicle development. A previous study demonstrated that OsGA20ox1 contributed to increasing grain number and yield by regulating cytokine activity in rice panicle meristems [58]. In addition, the expression data of GAox genes under cold or salt stress were retrieved and analyzed. The results showed that the majority of GAox genes responded to cold or salt stress and the expression of OsGA2ox-3, OsGA20ox-2, ZmGA2ox-7a, ZmGA2ox-6, and ZmGA2ox-10 changed significantly. Previous studies had proved that the application of GAox genes can improve plant stress resistance [20,68,69]. For example, the ectopic expression of GA2ox6 mutants can result in reduced plant height, increased yield, and improved abiotic and biotic stress tolerance in transgenic rice [68]. Therefore, future research on the function of GAox genes may provide valuable gene resources for stress-resistant breeding in Gramineae.

Conclusions
In this study, 105 gibberellin oxidase genes were identified in Gramineous crops, which can be subdivided into four subfamilies and nine OG groups. The similar number of GAox genes in the six species suggested that they may have evolved prior to the isolation of six Gramineous crop, and mainly expanded with WGD/SD expansion methods as well as undergone purification selection. The cis−acting elements analysis indicated that GAox genes may respond to growth and development, stress, hormones, and light signals. Moreover, the expression profiles of rice and maize showed that GAox genes were pre-dominantly expressed in the panicles of the above two plants and the expression of sever-al GAox genes was significantly induced by salt or cold stresses. In summary, this study provides a comprehensive insight into the evolution of GAox genes in six Gramineous crops.  Figure S1. The figure of phylogenetic tree constructed by RaxML. Figure S2. Characteristics of 15 conserved motifs of GAox genes in Gramineous crops. Figure S3. The evolutionary tree and conserved motifs of eight OsGA20oxs. Table S1. The primers of this experiment. Table S2. Summary of information on GAox genes in six species. Table S3. The members of OG groups. Table S4. Annotation of conserved motifs in all GAox proteins. Table S5. The expression values of GAox genes in rice and maize tissues. Table S6. The expression values of GAox genes in rice and maize tissues under salt or cold stress.
Author Contributions: C.Z. and Y.L. conceived and designed the experiments. C.Z. and X.N. performed all of the experiments, analyzed the data, prepared the figures and tables, and wrote the paper. X.N., T.S., and X.L. performed parts of the experiments, figures, and tables. W.K. and X.D. provided comments during the writing. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All data generated or analyzed during this study are included in this published article.

Conflicts of Interest:
The authors declare no conflict of interest.