Evolutionary Analysis of GH3 Genes in Six Oryza Species/Subspecies and Their Expression under Salinity Stress in Oryza sativa ssp. japonica

Glycoside Hydrolase 3 (GH3), a member of the Auxin-responsive gene family, is involved in plant growth, the plant developmental process, and various stress responses. The GH3 gene family has been well-studied in Arabidopsis thaliana and Zea mays. However, the evolution of the GH3 gene family in Oryza species remains unknown and the function of the GH3 gene family in Oryza sativa is not well-documented. Here, a systematic analysis was performed in six Oryza species/subspecies, including four wild rice species and two cultivated rice subspecies. A total of 13, 13, 13, 13, 12, and 12 members were identified in O. sativa ssp. japonica, O. sativa ssp. indica, Oryza rufipogon, Oryza nivara, Oryza punctata, and Oryza glumaepatula, respectively. Gene duplication events, structural features, conserved motifs, a phylogenetic analysis, chromosome locations, and Ka/Ks ratios of this important family were found to be strictly conservative across these six Oryza species/subspecies, suggesting that the expansion of the GH3 gene family in Oryza species might be attributed to duplication events, and this expansion could occur in the common ancestor of Oryza species, even in common ancestor of rice tribe (Oryzeae) (23.07~31.01 Mya). The RNA-seq results of different tissues displayed that OsGH3 genes had significantly different expression profiles. Remarkably, the qRT-PCR result after NaCl treatment indicated that the majority of OsGH3 genes play important roles in salinity stress, especially OsGH3-2 and OsGH3-8. This study provides important insights into the evolution of the GH3 gene family in Oryza species and will assist with further investigation of OsGH3 genes’ functions under salinity stress.


Phylogenetic Analysis
Multiple sequences alignments with Arabidopsis and O. sativa ssp. japonica were conducted by Clustal W, separately. A phylogenetic tree was generated by MEGA 6.0 using the Neighbor Joining (NJ) method with 1000 bootstrap replicates [27,28,30,31]. Subsequently, GH3 sequences of five Oryza species were systematically named based on the clustering results and names from previous studies [11,12,31].

Analysis of Gene Structure and Conserved Motifs
The exon/intron structure of GH3 genes was analyzed by comparing the coding DNA sequences (CDS) and the genomic sequences using the GSDS 2.0 (http://gsds.cbi.pku.edu.cn/). The Multiple Expectation Maximization for motif Elicitation (MEME, http://meme-suite.org/tools/meme) tool was used to predict conserved motifs of GH3 proteins with these parameters: the number of motifs (20) and other parameters (default values) [27,29,30]. Gene structure and conserved motifs were visualized using TBtools [32].

Expression Analysis and Co-Expression Network Analysis of OsGH3 Genes
A co-expression network in different tissues was constructed based on RNA-seq datasets using the Comparative Co-Expression Network Construction and Visualization tool (CoExpNetViz, http://bioinformatics.psb.ugent.be/webtools/coexpr/) with previously reported parameters [27]. The Co-expression network was visualized using Cytoscape V.3.1.0. The correlation coefficient >0.50 or < −0.50 was limited.

Expression Analysis of OsGH3 Genes under Salinity Stress by qRT-PCR
Primers of OsGH3 genes were designed by Primer 5.0 in specific regions or 3'-UTR regions (Primers in Table S1) [29]. The qRT-PCR reaction (10 µL) was formulated using ChamQ™ SYBR ® Color qPCR Master Mix (Vazyme, Shanghai, China). qRT-PCR was carried out in 96-well plates on a CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Ubi (LOC_Os03g13170, encodes the ubiquitin fusion protein) was used as an internal control. The average threshold cycle (Ct) from three biological replicates was used to determine the fold change of OsGH3 gene expression by the 2 −∆∆CT method [29].

Identification and Classification of the GH3 Gene Family
A total of 13, 13, 13, 12, and 12 members were identified in O. sativa ssp. indica, O. rufipogon, O. nivara, O. punctata, and O. glumaepatula, respectively. Based on the homologous sequence cluster result with Arabidopsis and O. sativa ssp. japonica, all GH3 genes were grouped into two groups (Table 1 and Figure 1): group1 and group2. In O. sativa ssp. indica, O. rufipogon, and O. nivara, four members of GH3s belonged to group1, and nine members of GH3s belonged to group2. In O. punctata and O. glumaepatula, four members of GH3s belonged to group1, and eight members of GH3s belonged to group2 (Table 1 and Figure 1). We found that the classification, chromosome locations, and number of GH3 genes were strictly conservative among these six Oryza species, and small differences in the number of GH3 genes could be produced from assembly or sequencing errors because these two incomplete sequences were also identified, namely OpGH10 (ID: OPUNC07G18660.1) in O. punctata and OgGH7 (ID: OGLUM06G14160.1) in O. glumaepatula. Considering this reason, these two incomplete sequences were removed in further research. Interestingly, no members of group III were found in five Oryza species, while 10 members of group III were found in Arabidopsis (Figure 1).

Gene Structure and Conserved Motif Analysis
Earlier studies have shown that gene structure diversity can provide the primary power for the evolution of multigene families [7,[27][28][29]36]. Thus, an exon/intron analysis was performed to obtain more insights into the structural diversity of GH3s in Oryza species ( Figure 2). The analysis results showed that the intron number of GH3s in six Oryza species/subspecies ranged from 1 to 7, GH3-4 contained the fewest introns (1), and GH3-13 contained the most introns (3-7) ( Figure 2). Furthermore, conserved motifs of 63 GH3 proteins were analyzed by MEME. As a result, 20 conserved motifs were identified, and the 63 GH3 proteins showed a similar conserved motifs arrangement. Notably, we found that GH3s from the same group showed differences in the number and the length of exons/introns. These results suggest that the gene function from the similar group (B) A phylogenetic tree of GH3 protein sequences from six Oryza species/subspecies. Clustal W is used for multiple sequence alignment. MEGA 6.0 is adopted for phylogenetic reconstruction by using the Neighbor Joining (NJ) clustering method. Bootstrap numbers (1000 replicates) are shown. Different colors of circles represent different subfamilies. The different species are indicated by different shaped markers.

Gene Structure and Conserved Motif Analysis
Earlier studies have shown that gene structure diversity can provide the primary power for the evolution of multigene families [7,[27][28][29]36]. Thus, an exon/intron analysis was performed to obtain more insights into the structural diversity of GH3s in Oryza species (Figure 2). The analysis results showed that the intron number of GH3s in six Oryza species/subspecies ranged from 1 to 7, GH3-4 contained the fewest introns (1), and GH3-13 contained the most introns (3-7) ( Figure 2). Furthermore, conserved motifs of 63 GH3 proteins were analyzed by MEME. As a result, 20 conserved motifs were identified, and the 63 GH3 proteins showed a similar conserved motifs arrangement. Notably, we found that GH3s from the same group showed differences in the number and the length of exons/introns. These results suggest that the gene function from the similar group has diversified. In short, the gene structure and conserved motif analysis of GH3s strongly supports the reliability of the group classification of GH3s in Figure 1. , and exon/intron structure (C) of the GH3 genes in six Oryza species/subspecies. A: Sequence alignments and the NJ-Phylogenetic trees were made using ClustalW and MEGA 6.0, respectively. A bootstrap number (1000 replicates) is adopted. The red and green colors in the phylogenetic tree represent group1 and group2, respectively. B, C: The widths of the grey bars represent the relative lengths of genes and proteins. The green boxes and grey lines display exons and introns, respectively.

Chromosome Locations, Duplication Events, Selection Pressure, and Microsynteny Analysis
The chromosome locations result showed that 13 GH3 genes were unevenly mapped on 12 chromosomes among the six species and subspecies. Four GH3 genes (30.77%) mapped on chromosome (Chr) 7, three GH3 genes (23.08%) mapped on Chr1 and Chr5, two GH3 genes (15.38%) mapped on Chr11, and one GH3 gene (7.69%) mapped on Chr6. No GH3 gene was found on Chr2, Chr3, Chr4, Chr8, Chr10, and Chr12 ( Figure 3A-F). These uneven distribution patterns of the GH3 gene family have also been observed in Arabidopsis, maize, tomato, and potato [5,6,[8][9][10]. Moreover, we discovered six pairs of segmental duplication events and two pairs of tandem duplication events in the six Oryza species/subspecies ( Figure 3). Interestingly, GH3-1 and GH3-4 segmental duplication events were found in every species/subspecies ( Figure 3). However, GH3-9 and GH3-10 tandem duplication events were only found in O. sativa ssp. indica and O. nivara (Figure 3B,E,G). Our findings suggest that duplication events were major factors determining the expansion of the GH3 gene family. Next, Ka/Ks values of duplicated GH3 gene pairs were calculated to evaluate the driving force underlying the GH3 gene's evolution. Ka/Ks >1, <1, and =1 mean a positive selection, a negative selection, and a neutral selection, respectively. The Ka/Ks results showed that the Ka/Ks values of eight duplicated GH3 genes ranged from 0.1245 to 0.2070 and were less than 1 ( Table 2). These results demonstrated that the duplicated GH3 genes were under a strong negative selection during the evolution process. The segmental duplication events of these six gene pairs were estimated to occur between 23.20 and 31.01 Mya (Table 2). Besides this, to further understand the evolutionary process of the GH3 genes in Oryza species, a microsynteny analysis was conducted among the six Oryza species/subspecies. In total, 169 collinear gene pairs were identified (Figure 4, Table S2). findings suggest that duplication events were major factors determining the expansion of the GH3 gene family. Next, Ka/Ks values of duplicated GH3 gene pairs were calculated to evaluate the driving force underlying the GH3 gene's evolution. Ka/Ks >1, <1, and =1 mean a positive selection, a negative selection, and a neutral selection, respectively. The Ka/Ks results showed that the Ka/Ks values of eight duplicated GH3 genes ranged from 0.1245 to 0.2070 and were less than 1 ( Table 2). These results demonstrated that the duplicated GH3 genes were under a strong negative selection during the evolution process. The segmental duplication events of these six gene pairs were estimated to occur between 23.20 and 31.01 Mya ( Table 2). Besides this, to further understand the evolutionary process of the GH3 genes in Oryza species, a microsynteny analysis was conducted among the six Oryza species/subspecies. In total, 169 collinear gene pairs were identified ( Figure 4, Table S2).

Analysis of Cis-Elements in OsGH3 Genes
The Cis-elements present in the stress-responsive gene promoters can provide an important insight into the stress response of plants [28,41,42]. Thus, the PlantCare database was used to identify the Cis-elements present in the promoter regions in the OsGH3 genes. Our results showed a high frequency of occurrence of Cis-elements in OsGH3 genes ( Figure 5) and that 47 Cis-elements were identified. These Cis-elements can be divided into 3 primary categories, including 20 secondary categories based on the previously described functional categories [43]. The cis-elements in the growth and development primary category showed a higher frequency of occurrence than that in the stress response and phytohormone response primary category ( Figure 5A,B). In the growth and development primary category, the number of light responsive/responsiveness secondary category elements exceeded 160 and contained 23 types of elements ( Figure 5A). In the phytohormone response primary category, the MeJA-responsiveness secondary category was the top secondary category (64), including the CGTCA-motif and the TGACG-motif, followed by the abscisic acid responsiveness secondary category (55), including the ABRE element ( Figure 5A). In the phytohormone response primary category, the top three secondary categories were the anaerobic induction, anoxic specific inducibility, and drought inducibility secondary categories, respectively ( Figure 5A). The result of the Cis-element position analysis showed that Cis-elements are unevenly distributed on all genes, and several Cis-elements were preferentially present on individual genes; for instance, the OsGH3-9 promoter regions had a lot of MeJA-responsive Cis-elements ( Figure 5C). Thus, we proposed that OsGH3 genes have the potential to improve stress tolerances because the OsGH3 genes contain several biotic/abiotic stress motifs in their promoter regions. It could also be further speculated that OsGH3-9 plays an important role in MeJA-related processes.
abscisic acid responsiveness secondary category (55), including the ABRE element ( Figure 5A). In the phytohormone response primary category, the top three secondary categories were the anaerobic induction, anoxic specific inducibility, and drought inducibility secondary categories, respectively ( Figure 5A). The result of the Cis-element position analysis showed that Cis-elements are unevenly distributed on all genes, and several Cis-elements were preferentially present on individual genes; for instance, the OsGH3-9 promoter regions had a lot of MeJA-responsive Cis-elements ( Figure 5C). Thus, we proposed that OsGH3 genes have the potential to improve stress tolerances because the OsGH3 genes contain several biotic/abiotic stress motifs in their promoter regions. It could also be further speculated that OsGH3-9 plays an important role in MeJA-related processes.

Functional Annotations Analysis of the OsGH3 Proteins
The results of the GO annotation and KEGG pathway analysis provide us with a better understanding of the biological functions of the different OsGH3 proteins. In this study, the GO annotation result revealed that 13 OsGH3 proteins were divided into 17 specific classes represented under the functional domains molecular functions, cellular components, and biological processes. Catalytic activity, metabolic process, and response to stimulus were predominant among the above-specified classes ( Figure 6). KEGG pathway results revealed that nine OsGH3 proteins (eight OsGH3 proteins belonged to group2, except for OsGH3-6) were enriched on the Auxin pathway, while three OsGH3 proteins (belonging to group 1) were enriched on the JA pathway. These results support previous reports that GH3 proteins from group I, with JA and/or SA-amido synthetase activity, use JA or SA as a substrate. GH3 proteins from group II, with IAA-amido synthetase activity, have Auxin-inducible expression profiles [5,9].
OsGH3 proteins belonged to group2, except for OsGH3-6) were enriched on the Auxin pathway, while three OsGH3 proteins (belonging to group 1) were enriched on the JA pathway. These results support previous reports that GH3 proteins from group I, with JA and/or SA-amido synthetase activity, use JA or SA as a substrate. GH3 proteins from group II, with IAA-amido synthetase activity, have Auxin-inducible expression profiles [5,9].

Expression Analysis of OsGH3 Genes in Different Tissues and under Salinity Stress
The expression analysis results revealed that OsGH3 genes showed different expression patterns in different tissues (Figure 7, Table S3). All OsGH3 genes clustered into two major groups based on expression levels ( Figure 7A). The genes in Group I (including eight members) showed low expression levels in all tissues, whereas the genes in Group II (including five members) had relatively high expression levels in some tissues. For example, OsGH3-5 showed a high expression

Expression Analysis of OsGH3 Genes in Different Tissues and under Salinity Stress
The expression analysis results revealed that OsGH3 genes showed different expression patterns in different tissues (Figure 7, Table S3). All OsGH3 genes clustered into two major groups based on expression levels ( Figure 7A). The genes in Group I (including eight members) showed low expression levels in all tissues, whereas the genes in Group II (including five members) had relatively high expression levels in some tissues. For example, OsGH3-5 showed a high expression level in anther and shoots. OsGH3-4 displayed relatively high expression levels in seed-5 DAP. OsGH3-2 exhibited a high expression level in anther and post-emergence inflorescence. OsGH3-8 exhibited a high expression level in anther, post-emergence inflorescence, and seed-5 DAP ( Figure 7A). Besides this, the co-expression network results indicated that there was a strong co-expression relationship network within the GH3 family's genes and that OsGH3-11, OsGH3-2, and OsGH3-4 were Hub genes in this network. level in anther and shoots. OsGH3-4 displayed relatively high expression levels in seed-5 DAP. OsGH3-2 exhibited a high expression level in anther and post-emergence inflorescence. OsGH3-8 exhibited a high expression level in anther, post-emergence inflorescence, and seed-5 DAP ( Figure  7A). Besides this, the co-expression network results indicated that there was a strong co-expression relationship network within the GH3 family's genes and that OsGH3-11, OsGH3-2, and OsGH3-4 were Hub genes in this network. It is common knowledge that salinity stress is a serious threat to crop yield worldwide [44]. The expression pattern under salinity stress can provide crucial clues to help us identify OsGH3 genes' functions. Hence, a qRT-PCR analysis at different time points after NaCl treatment was carried out. After NaCl treatment, the expression levels of all OsGH3 genes showed significant changes. The expression patterns of all OsGH3 genes in roots were different from that in leaves (Figure 8). In leaves, the expression levels of OsGH3-1, OsGH3-2, OsGH3-8, and OsGH3-10 were upregulated at all of the tested points and reached the highest at 6 h, 6 h, 24 h, and 24 h, respectively, whereas OsGH3-3, OsGH3-4, OsGH3-7, OsGH3-9, OsGH3-12, and OsGH3-13 were upregulated at two or three time points ( Figure 8A). In contrast, the expression levels of OsGH3-5 and OsGH3-11 were downregulated at two or three time points in leaves ( Figure 8A). The expression level of OsGH3-6 showed downregulation at all of the tested points in leaves ( Figure 8A). In roots, the expression levels of five genes, namely OsGH3-2, OsGH3-4, OsGH3-8, OsGH3-9, and OsGH3-12, were upregulated at all of the tested points and reached the highest at different points, while OsGH3-3, OsGH3-5, OsGH3-7, and OsGH3-10 were upregulated at two or three time points ( Figure 8B). Interestingly, OsGH3-1 and OsGH3-13 were upregulated at only one time point, respectively 12 h and 24 h, in roots ( Figure 8B). Yet, the expression levels of OsGH3-11 and OsGH3-13 showed downregulation at two or three time points ( Figure 8B) in roots. Specifically, we observed that OsGH3-2 and OsGH3-8 were upregulated at all tested points, namely 3 h, 6 h, 12 h, and 24 h in leaves and 3 h, 6 h, 12 h, and 24 h in roots, while OsGH3-6 was downregulated at 3 h, 6 h, 12 h, and 24 h in leaves and at 3 h, 6 h, and 12 h in roots ( Figure 8). Overall, the expression profiles in leaves and roots of these genes were different, indicating that OsGH3 genes have different roles under salinity stress. Correlation from weak to strong is shown by dotted line to solid line. Connectivity from weak to strong is shown from green to red. DAP, days after pollination.
It is common knowledge that salinity stress is a serious threat to crop yield worldwide [44]. The expression pattern under salinity stress can provide crucial clues to help us identify OsGH3 genes' functions. Hence, a qRT-PCR analysis at different time points after NaCl treatment was carried out. After NaCl treatment, the expression levels of all OsGH3 genes showed significant changes. The expression patterns of all OsGH3 genes in roots were different from that in leaves (Figure 8). In leaves, the expression levels of OsGH3-1, OsGH3-2, OsGH3-8, and OsGH3-10 were upregulated at all of the tested points and reached the highest at 6 h, 6 h, 24 h, and 24 h, respectively, whereas OsGH3-3, OsGH3-4, OsGH3-7, OsGH3-9, OsGH3-12, and OsGH3-13 were upregulated at two or three time points ( Figure 8A). In contrast, the expression levels of OsGH3-5 and OsGH3-11 were downregulated at two or three time points in leaves ( Figure 8A). The expression level of OsGH3-6 showed downregulation at all of the tested points in leaves ( Figure 8A). In roots, the expression levels of five genes, namely OsGH3-2, OsGH3-4, OsGH3-8, OsGH3-9, and OsGH3-12, were upregulated at all of the tested points and reached the highest at different points, while OsGH3-3, OsGH3-5, OsGH3-7, and OsGH3-10 were upregulated at two or three time points ( Figure 8B). Interestingly, OsGH3-1 and OsGH3-13 were upregulated at only one time point, respectively 12 h and 24 h, in roots ( Figure 8B). Yet, the expression levels of OsGH3-11 and OsGH3-13 showed downregulation at two or three time points ( Figure 8B) in roots. Specifically, we observed that OsGH3-2 and OsGH3-8 were upregulated at all tested points, namely 3 h, 6 h, 12 h, and 24 h in leaves and 3 h, 6 h, 12 h, and 24 h in roots, while OsGH3-6 was downregulated at 3 h, 6 h, 12 h, and 24 h in leaves and at 3 h, 6 h, and 12 h in roots ( Figure 8). Overall, the expression profiles in leaves and roots of these genes were different, indicating that OsGH3 genes have different roles under salinity stress.

Discussion
The GH3 gene family has been identified in several plants, such as Arabidopsis [9], O. sativa ssp. japonica [11][12][13], Z. mays [7], S. lycopersicum [6], P. patens [8], and S. moellendorffii [8], and it has important effects on plant growth, the plant developmental process, and various stress responses. Unfortunately, the evolutionary analysis of the GH3 gene family in Oryza species has not been well-studied to date. In recent years, the constantly released genomes of various wild rice species provide us with an opportunity to conduct a comprehensive analysis of this important gene family, including their gene structure, conserved motifs, a phylogenetic analysis, chromosome locations, gene duplication events, Ka/Ks ratios, and expression patterns.

Discussion
The GH3 gene family has been identified in several plants, such as Arabidopsis [9], O. sativa ssp. japonica [11][12][13], Z. mays [7], S. lycopersicum [6], P. patens [8], and S. moellendorffii [8], and it has important effects on plant growth, the plant developmental process, and various stress responses. Unfortunately, the evolutionary analysis of the GH3 gene family in Oryza species has not been well-studied to date. In recent years, the constantly released genomes of various wild rice species provide us with an opportunity to conduct a comprehensive analysis of this important gene family, including their gene structure, conserved motifs, a phylogenetic analysis, chromosome locations, gene duplication events, Ka/Ks ratios, and expression patterns.
The present study demonstrated that the number (13 members) and the gene structure of GH3 genes are strictly conservative across six Oryza species/species. The Oryza species has more GH3 genes than Marchantia polymorpha L. (two, mosses) [45], Physcomitrella patens (two, mosses) [46], Picea abies (seven, gymnosperms), Amborella trichopoda (six, angiosperms), Prunus persica (seven, eudicots), Capsicum annuum (11, eudicots), Vitis vinifera (nine, eudicots), Hordeum vulgare (five, monocots), and Phalaenopsis equestris (six, monocots) [8]. The number of GH3 genes in Oryza species is lower than that in Selaginella moellendorffii (18, ferns), A. thaliana (17, eudicots), Brassica rapa (38, eudicots), S. lycopersicum (17, eudicots), Glycine max (24, eudicots), Elaeis guineensis (16, monocots), and Musa acuminata (18, monocots) [8]. These results revealed that the GH3 genes have been expanded to different degrees in various plants. Additionally, there was no positive correlation between the number of GH3 genes and the size of the specie genome. For instance, the genome size of A. trichopoda is nearly twice that of O. sativa ssp. japonica, while O. sativa ssp. japonica (13) has a larger number of GH3 genes as compared with A. trichopoda (6). Therefore, we speculated that the difference in the number of OsGH3 genes is not related to the size of the genome. Besides this, we also observed an interesting phenomenon in which group III did not exist in mosses, ferns, gymnosperms, and angiosperms, whereas group III existed in eudicots and a few monocots [8], such as Brassicaceae plants. Considering these results, it can be deduced that group III might be the youngest group and originated from group I or group II and that group III could be related to the adaptation of plants to specific environments because group III has a species-specific expansion in various plants.
In the six Oryza species/subspecies, the same duplication events (GH3-1 and GH3-4) were found and duplication events of these six gene pairs were estimated to occur between 23.20 and 31.01 Mya. These results suggest that the expansion of the GH3 gene family might be attributed to duplication events and this expansion could occur in the common ancestors of Oryza species, resulting in similar structures and characteristics among the existing Oryza species. The calculated divergence time (23.20-31.01 Mya) of these duplication events is earlier than the differentiation time (~14 Mya) of Oryza species, while it is close to the origin time (~23.9 Mya) of the rice tribe (Oryzeae) [47][48][49][50]. Thus, it can be further deduced that the GH3 gene family produced these duplication events in a common ancestor of Oryzeae. Interestingly, tandem duplication events (GH3-9 and GH3-10) were only found in O. sativa ssp. japonica and O. nivara. This may be evidence to support the double domestication hypotheses of Chinese-cultivated rice subspecies [49,[51][52][53] that O. sativa ssp. japonica originated from O. rufipogon and O. nivara and that O. sativa ssp. indica originated from O. rufipogon [47].
The expression results of different tissues were consistent with previous findings on other species that GH3 genes displayed tissue-specific expression profiles [10,54]. For example, OsGH3-5 showed a high expression level in anther and shoots, implying that OsGH3-5 is involved in anther and shoot development. OsGH3-3 displayed relatively high expression levels in seed-5 DAP, which suggests that OsGH3-3 may be associated with seed development and growth. In addition, OsGH3-2, OsGH3-4, and OsGH3-8 had a high expression level in post-emergence inflorescence, indicating that these three genes may work together on post-emergence inflorescence development. Previous studies have reported that several GH3 genes play crucial roles in biotic and abiotic stress response [5][6][7]. For example, the SbGH3 gene was expressed at a low level under a normal condition, whereas it was substantially enhanced under salt and drought stress [55]. In chickpea CaGH3-1 and -7, and in Medicago MtGH3-7, -8, and -9, expression levels were significantly enhanced under drought or/and salinity stress [7]. In this study, we also found that the expression levels of OsGH3 genes were different under salinity stress. The majority of OsGH3 genes showed upregulation at different time points after NaCl treatment, indicating that OsGH3 genes play important roles in salinity stress response. However, some OsGH3 genes showed different expression profiles in leaves and roots. It can be inferred that the OsGH3 genes have different roles in roots and leaves. In addition, several OsGH3 genes formed a co-expression relationship network. Coincidentally, several genes showed a similar expression trend in the same tissues after NaCl treatment, such as OsGH3-5 and OsGH3-6, in leaves. These results indicated that OsGH3 genes may collaborate with each other in the response to salinity stress. Importantly, OsGH3-2 and OsGH3-8 were significantly upregulated at all the tested points in leaves and roots. OsGH3-1 was upregulated at all the tested points in leaves, while it was upregulated at only one time point in roots. Conversely, OsGH3-9 and OsGH3-12 were upregulated at all the tested points in roots, while they were upregulated at two time points in leaves. Based on these results, we speculate that OsGH3-2 and OsGH3-8 play important roles in leaves and roots and that OsGH3-1 plays a greater role in leaves than in roots under salinity stress, while OsGH3-9 and OsGH3-12 play a greater role in roots than in leaves under salinity stress.
In summary, a systematic analysis of GH3 in six Oryza species/subspecies was performed. The results revealed that the gene structure, conserved motifs, phylogenetic analysis, chromosome location, gene duplication events, and Ka/Ks ratios of the GH3 family were strictly conservative among Oryza species. The expansion of the GH3 family might be attributed to segmental duplication and tandem duplication, and this expansion could have occurred in the common ancestors of Oryza species and can be traced back to the origin time (~23.9 Mya) of the rice tribe (Oryzeae) [47]. The tandem duplication events (GH3-9 and GH3-10) were only found in O. sativa ssp. japonica and O. nivara. This may be evidence to support the double domestication hypotheses of Chinese-cultivated rice subspecies [49]. Similar to previous reports [10,54], OsGH3 genes showed tissue-specific expression. In addition, the qRT-PCR result indicated that OsGH3 genes play vital roles under salinity stress.