To identify all the putative GATA gene members of rice, MSU rice genome annotation project (RGAP) release 7 database (http://rice.plantbiology.msu.edu/) was scanned with GATA pfam (http://pfam.xfam.org/) domain ID PF00320. List of genes retrieved from MSU rice genome database was further confirmed by BLASTP (protein BLAST) in three different databases: NCBI (https://blast.ncbi.nlm.nih.gov/), phytozyme V.11 (https://phytozome.jgi.doe.gov/pz/portal.html), and plant genome database release 187 (http://www.plantgdb.org/), using Arabidopsis GATA1 protein sequence as the reference. Redundant sequences were omitted manually. Functional domains in the full-length protein sequence were identified using pfam, Interpro (https://www.ebi.ac.uk/interpro/), and SMART (http://smart.embl-heidelberg.de/) databases. To avoid ambiguity because of multiple names, nomenclature of these GATA genes was retained same as that of Reyes et al. (2004). However, newly identified alternative spliced products were denoted as gene number extended with suffix “a” and “b” as suggested earlier by Pareek et al. (2006).

For locating the GATA members on rice chromosomes, CDS coordinates were retrieved from MSU RGAP database version 7 for each GATA gene and were placed on each of the rice chromosomes according to the physical location of the gene. Plant genome duplication database (http://chibba.agtec.uga.edu/) was used to search for the duplication events (segmental duplication and tandem duplication) of the OsGATA genes. Duplicated genes have been connected by dotted lines.

Multiple sequence alignment of protein sequences of only GATA domain was performed using ClustalW program available in MEGA 7. The Neighbor-Joining tree was generated based on the MUSCLE alignment of the full-length OsGATA protein sequences using Jones-Taylor-Thornton (JTT) model under default setting in the MEGA 7 program. To infer phylogeny, 1,000 bootstrap replicates were taken. Gene structure display server version 2.0 (http://gsds.cbi.pku.edu.cn) was used to analyze the gene structure and for calculating the number of exon and intron in the gene sequence. The cis-acting elements in the promoter region of OsGATA genes were deduced using PlantPan version 2 (http://plantpan2.itps.ncku.edu.tw/).

Seeds of Oryza sativa L. cv IR64 and landrace Pokkali were rinsed thoroughly in sterile water, germinated in hydroponic set-up and raised for 7 days on half-strength Yoshida medium, at 28 ± 2°C for 12 h light and dark cycles in a plant growth chamber. For stress treatment, 7 days old seedlings were transferred to half-strength Yoshida medium supplemented with either 200 mM NaCl (for salinity), 20% PEG (for drought) or 100 μM ABA and shoot samples were harvested after 4 and 24 h of stress application. Un-treated seedlings growing in half-strength Yoshida medium were taken as control. After harvesting, samples were frozen immediately in liquid nitrogen and stored at −80°C until further use.

Total RNA was isolated from the shoots of seedlings using TRIzol Reagent (Life Technologies, USA). For extraction, 100 mg tissue was homogenized to a fine powder with liquid nitrogen using pre-chilled mortar and pestle. RNA was extracted as described earlier by Soda et al. (2013). Purity and integrity of the total RNA was analyzed using spectrophotometry (Thermo Scientific, USA) and denaturing agarose gel electrophoresis respectively. The quality of RNA was checked by A₂₆₀/A₂₈₀ ratio and samples having A₂₆₀/A₂₈₀ > 1.8 were used for further analysis.

First-strand cDNA was synthesized using first strand cDNA synthesis kit (Fermentas Life Sciences, USA) as per manufacturer's instructions. For removing genomic DNA contamination, DNAse (Epicenter, USA) treatment was done before proceeding for the cDNA synthesis as described earlier (Soda et al., 2013).

Primers for OsGATA genes were designed from the region corresponding to the junction of unique 3′ UTR region and CDS sequence, using Primer Express 3.0 software (Applied Biosystems, USA). In the case of alternative splice variants, primers were designed from a unique region within the CDS sequence. The uniqueness of each primer pair to amplify a selected gene was confirmed by BLASTN using the RGAP database and NCBI databases.

The qRT-PCR analysis was performed with a Sequence Detection System ABI Prism 7500 (Applied Biosystems, USA). Reactions (final volume, 10 μl) were set up with the 2X SYBR Green PCR Master Mix (Applied Biosystems, USA), 3 μl cDNA sample and 0.5 mM of gene-specific forward and reverse primers. All the PCR reactions were performed under the following conditions: 2 min at 50°C, 10 min at 95°C and 40 cycles of 15 s at 95°C, 1 min at 58–62°C (annealing temperature range for different genes) and 30 s at 72°C in 96-well optical reaction plates (Applied Biosystems, USA). The specificity of the amplification was tested by dissociation curve analysis. Three technical replicates were analyzed for each sample and the data analysis was performed using SDS 1.4 software (Applied Biosystems, USA). For data normalization, the rice eukaryotic elongation factor 1 alpha (eEF-1α) was taken as internal control. Transcript abundance of the selected group of genes was analyzed by qRT-PCR using the 2^−ddCT and 2^−dCT method for the calculation of fold change and transcript abundance respectively (Livak and Schmittgen, 2001).

To analyze the expression of the OsGATA genes at different developmental stages of rice, publicly available microarray database (https://www.genevestigator.com/gv/) was scanned with locus ID listed in Table 1. Expressions of these GATA genes were analyzed at germination, seedling, tillering, stem elongation, booting, heading, flowering, milking, and dough stages of rice plant development.

BLASTP search in NCBI using full-length protein sequence from Arabidopsis GATA1 as query identified 35 sequences which contain at least one GATA zinc-finger domain (Table 1). Further, MSU RGAP database version 7 was scanned for putative numbers of GATA genes using GATA domain ID PF00320 retrieved from pfam. Domain search also yielded 35 putative GATA transcription factors encoded by 28 gene loci. Our analysis yielded additional 7 OsGATA transcripts which were not reported earlier. Protein BLAST searches in RGAP database also yielded similar results. These 28 GATA genes were named as OsGATA1-OsGATA28 as described earlier by Reyes et al. (2004) (Table 1). The alternative spliced forms were named as “a” and “b” along with the GATA gene number (Table 1).

All the 35 GATA proteins contain at least one conserved GATA domain with a typical CX₂CX₁₈₋₂₀CX₂ zinc-finger motif except OsGATA8b which has partially truncated zinc-finger loop. Protein domain analysis using pfam, SMART, INTERPRO databases confirmed that two of the GATA genes have more than one GATA domains in the encoded protein sequences. OsGATA26 has two, and OsGATA24 has three and one truncated zinc-finger loop in their encoded proteins. Protein sequences encoded by 15 GATA genes that contained accessory domains other than GATA might play additional roles in different physiological responses (Table 1). Among them, six members, OsGATA17, OsGATA18, OsGATA19, OsGATA20, OsGATA22, and OsGATA23 possess CX₂CX₂₀CX₂ zinc-finger loop in GATA domain. While rest of the GATA members contain a CX₂CX₁₈CX₂ type of domain structure. All the rice GATA members are listed in Table 1 along with the gene nomenclature, domain details, and amino acid length. However, we found many differences in the amino acid length as well as in exon/intron structure of GATA TFs from earlier reported information. Predicted amino acid length of OsGATA1, OsGATA2, OsGATA5, OsGATA6, OsGATA10, OsGATA13, OsGATA16, OsGATA20, OsGATA21, OsGATA22, OsGATA23, and OsGATA26 are 387, 431, 376, 386, 142, 225, 390, 292, 450, 732, 742, and 383 respectively but earlier reports by Reyes et al. (2004) showed the amino acid length as 386, 387, 390, 387, 140, 155, 348, 332, 303, 778, 786, and 415 respectively. In addition to amino acid length, we found differences in the number of exons. Predicted exon numbers in the OsGATA5, OsGATA7, OsGATA10, OsGATA13, OsGATA16, OsGATA19, OsGATA20, OsGATA22, and OsGATA23 are 2, 2, 2, 5, 3, 8, 6, 5, and 5 respectively which were previously reported as 3, 3, 3, 3, 4, 9, 8, 7, and 6. This variation could be because of rapidly evolving genomic data and availability of refined annotation tools in the rice genome database version 7.

The GATA TFs vary in amino acid length from 101 to 742 with a predicted isoelectric point (pI) ranging from 4.56 to 10.03 and molecular weight ranging from 10.98 to 84.43 kDa. OsGATA8b was found to be the smallest protein having amino acid length 101 and molecular weight 10.98 kDa. On the other hand, OsGATA23a was the largest protein with an amino acid length of 742 and molecular weight of 84.43 kDa.

The OsGATA family members are randomly distributed on all the rice chromosomes, except VIII and IX (Figure 1). Maximum GATA genes i.e., six have been found to be present on chromosome III. On the other hand, only one each GATA gene has been annotated each on chromosomes VII and XI. The number of GATA genes vary from two to four on other rice chromosomes. OsGATA5 and OsGATA21 were clustered on chromosome IV between 27 and 27.2 Mb segments. OsGATA14 and OsGATA15 were present on chromosome V between 28.2 and 28.9 Mb region (Figure 1). Gene duplication has always been one of the well-known basis for the expansion of a gene family. Duplication can be either tandem; if duplicated genes are located on the same chromosome and closely linked or segmental; if duplicated genes are located on different chromosomes. We have observed eight segmental duplication events between OsGATA gene members (Figure 1) and one tandem duplication between OsGATA18 and OsGATA19 located between 27.5 and 30 Mb region of chromosome III (Figure 1).

Further, to compute the evolutionary distance between the genes, a Neighbour-joining tree was constructed in the Mega 7 program using the Jones Taylor Thornton (JTT) model. In this analysis, proteins with similar kind of domains got clustered in one group (Figure 2). In the course of setting up a new structural classification criteria for monocots, we have re-categorized the GATA proteins on the basis of their gene structure, the number of GATA domains, the position of GATA domain, and accessory domains (Figure 2). On the basis of homology in the GATA domain as well as the presence of accessory domain other than GATA, all the GATA genes have been subdivided into seven subfamilies (Figure 2). Typical domain structures of these TFs belonging to diverse subfamilies are presented in Figure 3. Subfamily-I has eight gene members including OsGATA1, OsGATA2, OsGATA3, OsGATA4, OsGATA5, OsGATA6, OsGATA7, and OsGATA25 (Figure 2). All these GATA proteins carry a single GATA domain at the C-terminal end (Figure 3). Among them, OsGATA1, OsGATA3, OsGATA6, OsGATA7, and OsGATA25 show the highest homology within the GATA domain (Figure 5). Subfamily-II is the largest and comprises of nine GATA genes, OsGATA8, OsGATA9, OsGATA10, OsGATA11, OsGATA12, OsGATA13, OsGATA14, OsGATA15, and OsGATA16 with the GATA domain being centrally located. OsGATA8b, one of the alternative splice variant of OsGATA8, contains a partially truncated GATA domain. Members of the subfamily—II are complex in terms of their domain architecture (Figure 3) and have been functionally sub-categorized as B-class of GATA transcription factors in rice and Arabidopsis (Behringer et al., 2014; Behringer and Schwechheimer, 2015). OsGATA9, OsGATA14, and OsGATA15 possess a HAN (HANABA TARANU) domain at the N-terminal of the protein (Figure 3). On the other hand OsGATA8, OsGATA10, OsGATA11, OsGATA12, OsGATA13, and OsGATA16 possess a highly conserved LLM (leucine-leucine-methionine) domain at the C-terminal of the protein (Figure 3). Members of the subfamily-III include OsGATA27 and OsGATA28 (Figure 2) which contain an N-terminal located GATA domain. Though they possess similar domain structure but differ in their gene structure (Figure 4). These are intronless genes (Figure 4). The lone member of subfamily-IV, OsGATA21 came out as an outlier with an extreme N-terminal GATA domain (Figure 3). OsGATA21 possesses six exons and five introns (Figure 4). The two members of subfamily-V, OsGATA26 and OsGATA24 contain unique 2 and 3 & 1 truncated GATA domains respectively, in the encoded protein. Subfamily-VI comprises of GATA genes which encode for GATA protein having GATA domain along with two accessory domains namely TIFY and CCT (Figure 2). OsGATA17, OsGATA18, OsGATA19, and OsGATA20 belong to subfamily-VI. The alternative splice variants of the genes OsGATA17 (a and b), OsGATA18 (a and b), and OsGATA19 (a and b) possess typical GATA domain with CX₂CX₂₀CX₂ zinc-finger loop. Gene structure of the members of this subfamily is complex, having 6–7 exons in the coding sequence (Figure 4). The GATA subfamily-VII comprises of only two members, OsGATA22 and OsGATA23 (Figure 2). The coding sequence for these genes are interrupted by four introns and hence possess five exons (Figure 4). Typical domain structure of this subfamily includes GATA, FAR1, MULE, and SWIM domain (Figure 3). Both FAR1 (FAR Red Impaired Response1) and MULE (Mutator-like transposases) domains show sequence homology and possess C₂H₂ zinc-finger-like motif and SWIM domain (found in SWI2/SNF and MuDR transposases). Further, to analyze the conserved amino acid residues in the GATA zinc-finger loop; we have carried out multiple sequence alignment of the GATA domains from all the peptide sequences (Figure 5). In the case of OsGATA26, both the GATA domains were kept in analysis and numbered as OsGATA26_1 and OsGATA26_2. Similarly, all the three GATA domains of OsGATA24; OsGATA24_1, OsGATA24_2, and OsGATA24_3 were aligned along with the other GATA domains. Partially truncated GATA domains of OsGATA8b and OsGATA24_4 were not included in this analysis. Close inspection of the aligned protein sequences revealed that apart from the conserved Cys residues at Cys-1, Cys-4, Cys-25, and Cys-28 in the zinc-finger loop, few amino acid residues in between the Cys-4 and Cys-28 are also conserved. The residues Thr-11, Pro-12, Gly-17, Pro-18, Lys-24, Asn-26, and Ala-27 (Figure 5), contribute to the formation of α-helix in the zinc-finger loop; suggesting their role in maintaining the structural integrity of the domain.

To comment on the possible roles of OsGATA family members in abiotic stress response, transcript abundance of OsGATA genes were analyzed in two contrasting rice genotypes, IR64 and Pokkali in response to distinct abiotic stresses such as salinity, drought, and stress responsive phytohormone ABA (Figure 6). To check the expression of all the 35 transcripts, unique primer combinations (Table S2), from the junction of 3′ UTR and CDS sequence were designed. However, OsGATA6, OsGATA7, OsGATA9, OsGATA14, OsGATA15, OsGATA19b, OsGATA21a, OsGATA24, OsGATA27, and OsGATA28 could not be amplified from any of the cDNA used in this analysis; therefore these were kept out of the expression analysis. Expression data of remaining OsGATA genes has been presented in the form of heat map (Figures 6A–F).

The qRT-PCR analysis revealed unique findings for the OsGATA family expression. Differential accumulation of the OsGATA transcripts was observed in IR64 and Pokkali (Figures 6A–I). Basal level expression of OsGATA2b, OsGATA8b, OsGATA11, OsGATA16, OsGATA17b, OsGATA20, OsGATA22, OsGATA23b, and OsGATA25 were comparatively higher in both the rice genotypes (Figures 6A–F). However, expression of some of the genes was found to be genotype specific. GATA members such as OsGATA17a, OsGATA18a, OsGATA18b, OsGATA19a, and OsGATA21b were highly expressed in IR64 under control conditions. On the other hand, expression of OsGATA1 and OsGATA10 were higher only in Pokkali genotype under the control conditions (Figures 6D–F). Also, the transcript levels of the OsGATA2a and OsGATA13 in IR64 and Pokkali were not affected by any of the aforesaid stresses (Figures 6A–F).

OsGATA3 was upregulated in response to exogenous ABA in both the genotypes (Figures 6C,F). However, induction was more pronounced in Pokkali. OsGATA26 accumulated in response to salinity and drought in both the genotypes, but accumulation pattern varied with duration of stress given (Figures 6A,B,D,E). In Pokkali, the gene was upregulated after 4 h of salinity and drought stress. However, in IR64, under drought stress, transcripts accumulated at 24 h of stress while under salinity stress, induction is marked at 4 h (Figures 6A,B).

In IR64, members of subfamily-VI such as OsGATA17a, OsGATA17b, OsGATA18a, OsGATA18b, OsGATA19a, and OsGATA20 as well as OsGATA22, a member of subfamily-VII, maintained higher transcript level at 24 h of all the applied stresses (Figures 6A–F). OsGATA23a also maintained higher transcript level at 24 h but only under drought in IR64. Contrarily, in Pokkali, the expression pattern of the gene members from subfamily VI and VII varied with the duration of stress (Figures 6D–F). Though induction was observed with the onset of stress (4 h) in response to salinity, drought as well as ABA (Figures 6D–F); but the transcripts level declined as the stress continued for 24 h, in salinity and ABA (Figures 6D,F); finally maintaining a level higher than the control. Under drought, transcript levels maintained a similar accumulation at the end of 24 h as that of 4 h level (Figure 6).

In this study, differential regulation of the alternative spliced variants of some of the OsGATA genes was also identified (Figures 6A–I). The transcript level of OsGATA2a, a spliced variant of OsGATA2, remained unchanged in all the applied stresses for both the rice genotypes (Figures 6A–F). Interestingly, expression of OsGATA2b was induced in response to salinity, drought, and ABA in both the genotypes at early stress duration (4 h). However, responses varied from one genotype to another during the later duration of stresses. In case of IR64, OsGATA2b showed higher transcript levels at 24 h post salinity and drought while the expression was downregulated in the ABA treated samples (Figures 6A–C). On the other hand, downregulation was seen in the Pokkali samples for salinity as well as ABA stress treatment at the end of 24 h while a higher level was seen in case of drought stress imposed samples for the same duration of time (Figures 6C–F).

Expression of OsGATA8b, the alternative spliced variant of the gene OsGATA8, was higher than the OsGATA8a in both the genotypes under all the three conditions tested (Figures 6A–F). Similarly, splice variants of OsGATA23 also varied in their expression pattern in both the genotypes in response to abiotic stresses (Figures 6A–F). In IR64 as well as in Pokkali, OsGATA23a showed an induction at 4 h stress duration and continued to maintain higher transcript till the end of 24 h of salinity and drought but it was downregulated in response to ABA at 24 h in both the genotypes (Figures 6A–F). On the other hand, transcripts of OsGATA23b were relatively low as compared to the OsGATA23a in both the genotypes. It can be said that, though the gene showed induction at the onset of stress but ultimately its transcript levels declined as the stress prolonged for 24 h in both the genotypes (Figures 6A–F).

In terms of fold change, maximum induction i.e., 40-folds was observed for OsGATA8a and OsGATA26 in IR64 under salinity stress (Figure 6G). Furthermore, in Pokkali salinity induced transcript levels of OsGATA8b, OsGATA18a, and OsGATA23a by more than 100-folds (Figure 6G). Interestingly in IR64, drought stress and ABA application lead to induction of OsGATA23a by more than 150- and 35-folds respectively (Figures 6H,I).

Similarly in Pokkali, it was seen that drought and ABA modulate the expression of a set of genes that included OsGATA4, OsGATA8a, and OsGATA21b upregulating them by more than 100- and 35-folds with respect to control under drought and ABA treatment respectively (Figures 6H,I). Similarly, OsGATA11 and OsGATA26 show more than 200- and 100-fold change respectively in transcript level under the drought stress. On the other hand, OsGATA18a showed ABA-dependent upregulation up to 300-folds and OsGATA23a showed a 250-folds change in response to ABA (Figure 6I).

The expression of the OsGATA23a was also found to be higher as compared to the OsGATA23b. In our study, OsGATA23a was identified as multi-stress responsive gene as it showed maximum upregulation for salinity, drought as well as ABA stress. In both the genotypes, this gene was induced at 4 h of stress. In the case of IR64, induction was 40-fold while in Pokkali more than 60-fold induction was observed. Since ABA is a stress hormone, it not only governs stomatal opening but also acts as a master regulator of abiotic stress signaling. Our data suggested that some of the OsGATA genes are highly responsive to ABA, hence these OsGATA genes might also be interplaying an important role at the junction of stress signaling cascade.

To understand the regulation of these GATA transcription factors, we analyzed the cis-acting elements in the promoter region of OsGATA genes (Table S1). Binding sites of various stress responsive transcription factors such as AP2 (Apetala2), ERF (Ethylene Response Factors), bHLH (basic helix loop helix), bZIP (basic leucine zipper), MADF (myb/SANT-like domain in Adf-1), Myb, NAC (NAM, ATAF1/2 and CUC1/2), WRKY, and MADS box were found. Interestingly, we found that promoter region of almost all the OsGATA genes possess GATA binding sites indicating that the expressions of OsGATA genes might be regulated by GATA transcription factors themselves.

To comment on the role of OsGATA transcription factors in rice developmental processes, we analyzed the expression data from publicly available microarray database, Genevestigator (Table S3). Transcript abundance of distinct members of the OsGATA gene subfamilies at various developmental stages like germination, seedling, tillering, stem elongation, booting, heading, flowering, milking, and dough stage was checked (Figure 7). Interestingly, transcripts of members of the subfamily IV and VI were found to be comparatively abundant throughout all the developmental stages of the rice plant. On the other hand, members of the subfamily II showed huge variation in terms of fold change at the different stages from seedling to maturity. The members of subfamily I and II showed mixed expression pattern. Our analysis showed that the OsGATA12 gene was induced the most during seedling stage amongst all OsGATA genes. Expression of OsGATA22, members of subfamily VII, was observed to be low as compared to other GATA gene members (Figure 7). On the other hand OsGATA23 falls under higher expression group (Figure 7).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.