Allelic variants of OsSUB1A cause differential expression of transcription factor genes in response to submergence in rice

Background Flooding during seasonal monsoons affects millions of hectares of rice-cultivated areas across Asia. Submerged rice plants die within a week due to lack of oxygen, light and excessive elongation growth to escape the water. Submergence tolerance was first reported in an aus-type rice landrace, FR13A, and the ethylene-responsive transcription factor (TF) gene SUB1A-1 was identified as the major tolerance gene. Intolerant rice varieties generally lack the SUB1A gene but some intermediate tolerant varieties, such as IR64, carry the allelic variant SUB1A-2. Differential effects of the two alleles have so far not been addressed. As a first step, we have therefore quantified and compared the expression of nearly 2500 rice TF genes between IR64 and its derived tolerant near isogenic line IR64-Sub1, which carries the SUB1A-1 allele. Gene expression was studied in internodes, where the main difference in expression between the two alleles was previously shown. Results Nineteen and twenty-six TF genes were identified that responded to submergence in IR64 and IR64-Sub1, respectively. Only one gene was found to be submergence-responsive in both, suggesting different regulatory pathways under submergence in the two genotypes. These differentially expressed genes (DEGs) mainly included MYB, NAC, TIFY and Zn-finger TFs, and most genes were downregulated upon submergence. In IR64, but not in IR64-Sub1, SUB1B and SUB1C, which are also present in the Sub1 locus, were identified as submergence responsive. Four TFs were not submergence responsive but exhibited constitutive, genotype-specific differential expression. Most of the identified submergence responsive DEGs are associated with regulatory hormonal pathways, i.e. gibberellins (GA), abscisic acid (ABA), and jasmonic acid (JA), apart from ethylene. An in-silico promoter analysis of the two genotypes revealed the presence of allele-specific single nucleotide polymorphisms, giving rise to ABRE, DRE/CRT, CARE and Site II cis-elements, which can partly explain the observed differential TF gene expression. Conclusion This study identified new gene targets with the potential to further enhance submergence tolerance in rice and provides insights into novel aspects of SUB1A-mediated tolerance. Electronic supplementary material The online version of this article (doi: 10.1186/s12284-017-0192-z) contains supplementary material, which is available to authorized users.


Background
Rice (Oryza sativa L.) is an important food crop globally and the most advanced genetic monocot model amongst the cereal crops (Cantrell and Reeves, 2002). Biotic and abiotic stresses, such as drought and heat, are known to be detrimental to crop production. Rice production is additionally constraint by submergence stress during the rainy season with complete submergence and water stagnation affecting about 20 million hectares of rice fields in the tropics, causing significant yield and economic losses, and food insecurity (Septiningsih et al., 2012). Rice fields can be flooded with several meters of water for weeks and plants die within a few days from lack of oxygen and impaired photosynthesis (Gibbs and Greenway, 2003;Fukao and Bailey-Serres, 2004). Furthermore, in the attempt to escape the water, plants excessively elongate leaves thereby depleting their carbohydrate and energy reserves, and as water recedes, the plants lodge thus impairing any recovery growth (Singh, Singh, and Ram, 2001;Das, Sarkar, and Ismail, 2005;. Submergence-tolerant rice varieties were identified in the early 1960s and breeding efforts to introgress this trait into rice cultivars have been ongoing since the 1970s. Submergence tolerance was initially considered a complex trait involving many genes located in multiple quantitative trait loci (QTL) (Mackill, Ismail, Singh, Labios, and Paris, 2012). However, advances in molecular marker technologies and optimized phenotyping techniques aided in the identification of the major QTL Submergence tolerance 1 (Sub1) (Xu, Deb, and Mackill, 2004). Sub1 was identified from the aus-type rice landrace FR13A (Flood Resistance 13A) and it provided a breakthrough in the understanding of submergence tolerance mechanisms, enabling markerassisted breeding of submergence-tolerant rice Xu et al., 2006;Singh, Mackill, and Ismail, 2009;Mackill et al., 2012). Sub1 rice varieties are now widely grown across Asia enhancing plant survival in flooded fields with yield advantages of one ton per hectare or more (Septiningsih et al., 2009;Mackill et al., 2012).
Molecular and comparative sequence analysis of the Sub1 genomic region in submergence-tolerant and -intolerant rice varieties revealed the presence of a variable cluster of two to three ethylene-responsive TFs (ERFs) genes, namely SUB1A, SUB1B and SUB1C . In submergence-intolerant and moderately tolerant rice genotypes, the SUB1A-1 gene is either absent, due to an inversion-deletion, or present as an allelic variant, SUB1A-2 (Singh et al., 2010) (Additional file 1: Figure S1a).
Transgenic approaches identified SUB1A-1 as the major tolerance gene, since its constitutive expression conferred submergence tolerance to an intolerant rice variety (M202), which naturally lacks the SUB1A gene . Under complete submergence, the main phenotypic effect observed in SUB1A-1 overexpression and Sub1 near isogenic lines (NILs), is a significantly reduced elongation growth (Additional file 1: Figure S1b). The Sub1 plants assume a "quiescence" status which prevents excessive elongation growth thereby preserving carbohydrate reserves and preventing an energy crisis, as well as lodging once the water recedes Xu et al., 2006;Voesenek and Bailey-Serres, 2009).
In contrast to the SUB1A-mediated tolerance, deepwater rice responds to submergence with rapid GA-induced elongation growth, so that plants reach the water surface for access to light and oxygen Bailey-Serres and Voesenek, 2008;Singh et al., 2010). As is the case for SUB1A-1, rapid elongation growth in deepwater rice is also ethylene induced but regulated by different ERFs, namely SNORKEL1 (SK1) and SNORKEL2 (SK2) (Hattori et al., 2009;Nagai, Hattori, and Ashikari, 2010).
At the molecular level, SUB1A-1 functions by preventing ethylene-induced, GA-mediated elongation growth via enhancing the level of the GA-inhibitors SLENDER RICE 1 (SLR1) and SLENDER RICE LIKE 1 (SLRL1) (Fukao and Bailey-Serres, 2008). Further to this, presence of SUB1A-1 is associated with the negative regulation of genes involved in carbohydrate catabolism and cell elongation (via expansins) and positive regulation of genes involved in ethanolic fermentation, such as alcohol dehydrogenase Bailey-Serres and Voesenek, 2010). Recently some studies have revealed that SUB1A-1 is phosphorylated by MITOGEN-ACTI-VATED PROTEIN KINASE3 (MAPK3) and in turn physically interacts with MAPK3 and binds to its promoter (Singh and Sinha, 2016). Overall, the suppressed growth and metabolic adjustments during flooding confers a higher recovery rate to Sub1 plants compared to non-Sub1 plants as illustrated in the aerial photo of the IRRI demonstration plot (Additional file 1: Figure S1b).
The aforementioned molecular and physiological studies were conducted by comparing Sub1 NILs and SUB1A-1 overexpressing plants with their corresponding wild-type varieties (M202 and Liaogeng, respectively), which have the SUB1B and SUB1C genes but naturally lack the SUB1A gene Xu et al., 2006). However, an allelic variant of the SUB1A gene exists, SUB1A-2, which is present in a range of rice varieties Singh et al., 2010). Phenotyping of these SUB1A-2 varieties revealed a variable level of submergence tolerance ranging from about 4% to 40% plant survival after two weeks of submergence. SUB1A expression analysis subsequently showed that both, SUB1A-1 and SUB1A-2 are submergenceinducible and expressed at similar levels in leaves, stems and developing panicles. However, gene expression in internodes associated well with submergence tolerance, i.e., a higher level of SUB1A expression in internodes under submergence was observed in varieties with the Sub1A-1 allele compared to varieties with the SUB1A-2 allele (Singh et al., 2010) (Additional file 1: Figure S1c).
To further address differences between SUB1A-1 and SUB1A-2 allele effects, we analysed internode samples from IR64 (SUB1A-2) and the IR64-derived NIL IR64-Sub1, which carries the tolerant Sub1 locus including the SUB1A-1 allele (Singh et al., 2010). Quantitative real-time PCR (qPCR) targeting nearly 2500 rice TF genes based on the Nipponbare reference genome (Caldana, Scheible, Mueller-Roeber, and Ruzicic, 2007) was conducted with internode samples from plants submerged for 30 h and corresponding non-submerged control plants to (a) identify TFs that are differentially regulated upon submergence in IR64 and IR64-Sub1, (b) compare TFs in both genotypes to provide insights in the SUB1A regulatory pathway and (c) analyse the promoter region of SUB1A-1 and SUB1A-2 to identify single nucleotide polomorphisms (SNPs) that create putative cis-elements and could be binding sites for upstream TFs.

Plant material and submergence treatment
Seeds of the rice variety IR64 (accession number 66970) and the IR64-Sub1 NIL (accession number IR8419422-139) were provided by the International Rice Germplasm Collection (IRGC) at IRRI, Philippines. Seeds were incubated at 55°C for 5 days to break dormancy before pregermination in petri dishes at 37°C in an incubator. Three-day-old seedlings were transplanted into pots filled with soil substituted with 3 g ammonium sulphate. Four pots with two plants for each accession were grown for 75 days until booting/heading stage, i.e., the onset of internode elongation. At that stage, plants were completely submerged in a concrete tank filled with about 2 m of tap water for 30 h (Singh et al., 2010). Nonsubmerged control plants were kept under natural conditions in an IRRI screenhouse and sampled in parallel. From each plant, multiple internodes from three tillers were sampled from a total of four submerged and two control plants per genotype. Samples were immediately frozen in liquid nitrogen and stored at −80°C.

RNA isolation and cDNA preparation
Total RNA was isolated from the internode tissues of the control and submerged plants using TRIzol® (Invitrogen) according to the manufacturer's instructions. The quantity and quality of RNA samples were assessed using a Nanodrop ND-100 (Thermo Scientific) and by agarose gel electrophoresis. DNA contamination was removed by treating RNA samples with RNase-free DNase I according to the protocol provided (Promega). cDNA was synthesized from 5 μg of total RNA using Super-script™ III reverse transcriptase (Invitrogen, Germany) as per the manufacturer's protocol.

Data source for transcription factor sequences and primers
A qRT-PCR platform was used for expression profiling of the rice TF genes (Caldana et al., 2007). The initial source dataset used for the annotation and primer design was version 2.0 of the Rice Genome Annotation Project (RGAP) (http://rice.plantbiology.msu.edu/) and encompassed a total of 2508 (2487 unique) TF genes. TF sequences were extracted from version 1 of the Plant Transcription Factor Database (Riano-Pachon, Ruzicic, Dreyer, and Mueller-Roeber, 2007;Perez-Rodriguez et al., 2010). The TF genes were initially synchronized with version 5.0 of the RGAP Pseudomolecule and genome annotations resulting in 2221 unique genes with appropriate annotation in version 5.0. Twenty of these were represented in two distinct positions in the qPCR platform as internal controls. Based on version 5.0 of the RGAP Pseudomolecules, 266 unique genes out of 2487 missed annotations (Additional file 2: Table S1). The TF genes were further revised based on the most recent RGAP version 7.0 annotation and details are given in Additional file 2: Table S2. For the genes with missing annotations in Release 7, BLAST searches were performed taking the sequences of the PCR amplicons as queries against the RGAP database.

Transcription factor profiling
The qPCR platform used for this study had the primer pairs specific to the TF genes described above and were distributed on a total of seven 396 well-plates (Caldana et al., 2007). On each plate, twelve wells were reserved for three reference genes (RGs; represented in quadruplicate) and four wells were for negative water control. For the current study, the RGs 9631.m04973, 9629.m05807 and 9633.m03388 were used (Additional file 2: Table S3 and S4). Of these, 9633.m03388 (LOC_Os05g36290, actin) was selected as the best suitable RG for all subsequent calculations and normalizations. There were three technical replicates for each of the four samples (IR64 C , IR64 S , IR64-Sub1 C and IR64-Sub1 S , where C denoted control and S submergence treated plants). Real-time PCR analysis was based on SYBR Green (Applied Biosystems, Germany) and was conducted in a total volume of 5 μl (2.5 μl of SYBR Green master mix, 0.5 μl of cDNA (1.25 ng/μl) and 200 nM forward and reverse primers). PCR amplifications were carried out as described in (Caldana et al., 2007) in an ABI PRISM 7900 HT sequence detection system (Applied Biosystems, Germany).

Data processing and statistical analysis
Based on the obtained amplification curves, the Ct values (fractional cycle number at threshold), the value of PCR efficiency and the corresponding coefficient of determination (R 2 ) were calculated using LinRegPCR software (Ruijter et al., 2009) for all four samples (IR64 C , IR64 S , IR64-Sub1 C and IR64-Sub1 S ). Ct values for all reactions were log 2 transformed (Additional file 2: Table  S3 and S4). The average Ct values for RGs (Ct RG ) were calculated from Ct values of quadruplicates for each 396 well-plate separately for control and submerged conditions. The RG selected for all subsequent calculations was actin (LOC_Os05g36290). Based on the average Ct RG , the ΔCt for each selected gene (SG) was calculated as ΔCt = Ct SG -Ct RG . All reactions with R 2 < 0.995, reflecting low-quality amplifications and those designated as "undetermined" (did not deliver a Ct < 40) were excluded from further analyses. The obtained sets of ΔCt values for control and submerged conditions for each gene model (in three replicates) were examined for significant differential expression using LIMMA's moderated t-test (Smyth, 2005). LIMMA computes moderated t-statistics and log-odds of differential expression by empirical Bayes shrinkage of the standard errors towards a common value. False Discovery Rate (FDR) corrections of p-values were carried out using Benjamini and Hochberg method (Benjamini and Hochberg, 1995). A gene was considered differentially expressed if the moderated t-test resulted in a corrected p-value <0.01. Hence, the list of differentially expressed TF genes and corresponding log 2 FC (fold change) values under submergence stress were generated from both, the IR64 and the IR64-Sub1 submerged plants compared to the respective control plants. In addition, we examined the data for genotype-specific expression differences. For this, we used the contrast matrix referred from the LIMMA manual (https://www.bioconductor.org/packages/devel/ bioc/vignettes/limma/inst/doc/ usersguide.pdf ) as follows: Contrast.matrix <− makeContrasts (SubmergenceIR64 = IR64 S -IR64 C , SubmergenceIR64Sub1 = IR64Sub1 S -IR64Sub1 C , Genotype = (IR64Sub1 S + IR64Sub1 C ) -(IR64 S + IR64 C ) levels = design).

Results
In order to investigate the difference between submergence tolerance as mediated by the strong allele SUB1A-1 (present in IR64-Sub1) and the weak allele SUB1A-2 (present in IR64), expression of 2487 rice TF genes was quantified by qPCR in internodes of submerged and non-submerged control plants. Internodes were chosen for this study since differences in the expression of the SUB1A-1 and SUB1A-2 alleles were most pronounced in this tissue (Singh et al., 2010).

Update of TF gene annotations
The 2487 TF genes given in Additional file 2: Table S1 were synchronised to the latest version 7.0 of the RGAP Pseudomolecules and annotation of the rice genome (Kawahara et al., 2013). For 72 of the 266 genes with missing annotations, we could obtain annotations from version 7.0 (Additional file 2: Table S2) and one gene was marked as obsolete. Therefore, only 193 genes remained as not being annotated in version 7.0. Nine genes previously annotated in version 5 miss annotations in version 7, eight of these are now considered obsolete and one is without annotation in versions 6.1 and 7.0 of RGAP annotations. Hence, a total of 2284 TF gene models were updated with version 7.0 annotations (Table 1).

Differential gene expression analysis in IR64 and IR64-Sub1
The expression patterns of the TF genes in four samples (IR64 C , IR64 S , IR64-Sub1 C and IR64-Sub1 S ) were analysed by qPCR. Statistical comparisons for the identification of DEGs were performed only for those genes which possessed at least two calculated ΔCt values in both genotypes and both growth conditions (control and submerged). For performing differential expression analysis an intersect of 1823 TF genes, 1944 genes of the IR64 dataset and 1933 of the IR64-Sub1 dataset, were taken as input for LIMMA's moderated t-test comparisons (Fig. 1a). The resulting p-values, FDR corrected p-values and log 2 FC values for all gene models are provided in Additional file 2: Table S6.

Submergence-responsive TFs in IR64-Sub1
In IR64-Sub1, 19 TF genes were found to be submergence-responsive, of which thirteen genes were upregulated and six genes were downregulated (Fig. 1b, Table 2). The group of upregulated genes include: 6 MYB factors, 3 Zn finger TFs, 1 Heat Shock Factor (HSF), 1 bHLH, 1 bZIP, and 1 HOX TFs. And the seven (6 + 1) downregulated genes include: 2 AP2-type TFs, 2 MYBs, a WRKY, a TIFY and a RING-H2 finger gene. One of the seven downregulated genes, the MYB factor LOC_01g19330, was downregulated in response to 30 h of submergence in both IR64 and IR64-Sub1 genotypes and is placed in the category of a common gene (Fig. 2, Table 2). Studies of these TFs in relation to abiotic stress and hormonal response are detailed below and an overview of the data is provided in Fig. 3.
Only two of the above-mentioned genes have previously been reported to be upregulated under submergence and dehydration in IR64-Sub1, namely the bZIP gene (LOC_Os01g36220) and the Zn-finger gene (LOC_Os07g43740) (Fukao et al., 2011). The latter is also upregulated under anoxia, heat and chilling (Ham et al., 2013;Shankar et al., 2016).

Downregulated genes
In this category, two AP2-type TFs (LOC_Os07g22730, LOC_Os05g27930) have previously been reported as submergence-responsive (Jung et al., 2010). LOC_Os07g22730 belongs to the AP2 group IIIb and was slightly upregulated in M202-Sub1 after 6d of submergence but downregulated under anoxia. This gene is also upregulated under desiccation and salt stress (Shankar et al., 2016) and known to be highly expressed in roots and embryos (Rashid, He, Yang, Hussain, and Yan, 2012). The other AP2 factor (LOC_Os05g27930) codes for OsDREB2b, which was reported to be slightly downregulated under submergence in both, M202 and M202-Sub1 (Jung et al., 2010). OsDREB2b is a transcriptional activator and is upregulated in the nodes of droughtstressed plants (Todaka et al., 2012). It is also responsive to high and low-temperature stress and over-expressing lines exhibit increased plant survival under drought (Matsukura et al., 2010).
The two downregulated MYB TF genes (LOC_Os02g40530, LOC_01g19330) encode putative calmodulin-binding proteins and are known to play important roles in signal transduction regulating plant development and adaptation responses to different stress conditions (Chantarachot, Buaboocha, Gu, and Chadchawan, 2012). Both genes were reported to be upregulated under drought, desiccation and salt stress whereas downregulated when treated with GA (Chantarachot et al., 2012;Katiyar et al., 2012;Shankar et al., 2016;Xiang et al., 2017).
The other downregulated genes in this category code for OsWRKY76 (LOC_Os09g25060), OsTIFY11a/OsJAZ9 (LOC_Os03g08310) and a RING-H2 finger gene (LOC_Os02g35329). OsWRKY76 is constitutively expressed in vegetative and reproductive tissues and upregulated under desiccation, salinity, auxin, GA and MeJA treatments (Ramamoorthy et al., 2008;Shankar et al., 2016). In contrast, OsTIFY11a/JAZ9 is specifically expressed in young panicles (Ye et al., 2009). OsTI-FY11a/JAZ9 is widely stress-responsive and induced by JA, wounding, drought, desiccation, salt and cold stress but downregulated by ABA treatment (Shankar et al., 2016). It was shown that overexpression of TIFY11a/ JAZ9 enhances shoot growth in rice under salt and mannitol treatment (Ye et al., 2009). OsTIFY11a/JAZ9 interacts with OsCOI1 suggesting its involvement in JA signalling as a transcriptional regulator for salt stress tolerance (Wu, Ye, Yao, Zhang, and Xiong, 2015). No information was found for the RING-H2 finger gene.

Submergence-responsive TFs in IR64
There were 26 submergence-responsive TF genes specific to IR64, excluding the above-mentioned MYB Fig. 4 Differentially expressed transcription factor genes in IR64. An overview of the transcription factor genes and their related pathways is provided for all genes that are specifically submergence-responsive in IR64 which carries the SUB1A-2 allele. The figure is based on a literature review, see text for details factor LOC_01g19330 common to both genotypes. The genes encode: 5 AP2 domain TFs, 4 NAC, 4 WRKY, 4 MYB, 1 B3 DNA binding protein and a number of different Zn-finger proteins. All of these genes were downregulated under submergence ( Fig. 1b; Table 2) and an overview of the data is provided in Fig. 4. Documented studies of these TF genes are summarised below based on their response to submergence, hormone signalling pathways and other abiotic and biotic stresses.
The most responsive gene in IR64 was LOC_Os07g07690, an uncharacterized gene with a putative Jas-domain (pham16135), plant homeodomain (PHD) and Zn-finger domain (cd15539, pfam00628, smart00249). The second most downregulated gene, interestingly, was SUB1B (LOC_Os09g11480), one of the three ERF factors within the Sub1 QTL   (Table 2, Fig. 2 and Additional file 1: Figure S1a). SUB1B is reportedly not responsive to GA or ethylene Xu et al., 2006). The other AP2 TF genes include SUB1C (LOC_Os09g11460) and three additional ERFs (LOC_Os03g22170, LOC_Os07g47790 and OsERF922: LOC_Os01g54890). The SUB1C gene was reported to be induced by GA and ethylene and upregulated under submergence in the intolerant japonica variety M202, which naturally lacks the SUB1A gene . The same study also suggested that SUB1A-1 negatively regulates SUB1C expression. Our results are in agreement with that finding and suggest that SUB1A-2 also negatively regulates SUB1C under submergence (Fig. 2). However, downregulation of SUB1C was not observed in IR64-Sub1 suggesting an alternate regulatory pathway. One of the recent studies revealed that desiccation stress could also downregulate SUB1C and SUB1B in IR64 plants (Shankar et al., 2016). The three other AP2/ERF genes have also been reported from studies on submergence tolerance. OsERF067 (LOC_Os07g47790) was shown to be upregulated in M202 under submergence and in IR64 under dehydration. This gene was seen to be upregulated by SUB1A during submergence and dehydration in M202-Sub1 and IR64-Sub1 suggesting a downstream position in the SUB1A pathway (Fukao et al., 2011). Contrastingly, another study reported downregulation of OsERF067 in M202 under submergence but upregulation in M202-Sub1 (Jung et al., 2010). The same study showed a similar expression pattern for OsERF066 (LOC_Os03g22170). Both of these AP2 factors are downregulated under GA treatment (Xiang et al., 2017), involved in cytokinin response (Lara et al., 2004) and are positively responsive to anaerobic germination in Nipponbare (Magneschi and Perata, 2009). The gene OsERF922 was shown to be upregulated in M202-Sub1 after 6 days of submergence but downregulated in M202 after 1 day of submergence (Jung et al., 2010). OsERF922 is strongly induced by ABA, salt treatment, and biotic stress. A study on the crosstalk between abiotic and abiotic stress signalling through modulation of ABA levels further suggests that OsERF922 functions as a transcriptional activator (Liu, Chen, Liu, Ye, and Guo, 2012).
Of the 4 NAC genes downregulated in IR64, two have been reported to respond to various biotic and abiotic stresses, including submergence. OsANAC34 (LOC_Os03g04070) is upregulated under submergence, desiccation and by rice stripe virus infection but downregulated under GA (Nuruzzaman et al., 2010;Shankar et al., 2016;Xiang et al., 2017). The OsNAC6/SNAC2 (LOC_Os01g66120) gene contains ABRE and CBE elements in its promoter and is upregulated under drought, cold and submergence stress (Nuruzzaman et al., 2010). These two genes are responsive in Nipponbare, suggesting they are not under the control of SUB1A, which is absent in this rice variety. OsNAC6/SNAC2 was also reported to be upregulated under desiccation and salinity stress in a drought-tolerant rice (N22) mutant (Shankar et al., 2016).
The C3HC4 Zn finger (LOC_Os01g11460) and B-Box Zn finger OsBBX1 (LOC_Os01g10580) TF genes were upregulated under submergence in the intolerant variety M202, the former also in M202-Sub1, and in IR64 upon dehydration (Fukao et al., 2011). OsBBX1 is highly expressed in seeds, roots and the endosperm, and is induced by auxin, GA and cytokinin (Huang et al., 2012). The bZIP TF factor LOC_Os01g64020 was found to be downregulated under submergence in IR64 plants in our study and was reported to be nitrogen-responsive in studies in rice and Arabidopsis (Obertello, Shrivastava, Katari, and Coruzzi, 2015).
The induction of ethylene in response to submergence is known to contribute to a decline in ABA and an increase in GA thereby regulating elongation growth (see above). ABA in-turn is involved in the regulation of the crosstalk between JA and SA signalling pathways (Liu et al., 2012). A complex interplay and fine-tuning of hormonal responses is therefore likely to determine the level of submergence tolerance.
The other TF genes have not been reported to be ABA-or GA-responsive and/or responsive to submergence. OsBBI1 (LOC_Os06g03580) and OsONAC7 (LOC_Os11g05614) have been associated with biotic stress response. OsBBI1 encodes a RING finger protein with E3 ligase activity, which modifies cell wall defence and confers broad-spectrum disease resistance to blast fungus (Magnaporthe oryzae). Expression of OsBBI1 is also induced by the rice blast fungus, benzothiadiazole and SA (Li et al., 2011). OsONAC7 is expressed and highly upregulated in rice stripe virus (RSV) and tungro virus (RTSV) infected plants (Nuruzzaman et al., 2010). The NAC TF gene OsSNAC (LOC_Os07g37920) was shown to be downregulated under drought stress in the intolerant variety Nipponbare (Nuruzzaman et al., 2010) but upregulated under desiccation stress in the rainfed varieties N22 and Pokkali (Shankar et al., 2016).
Limited information is available for the remaining genes downregulated in IR64. An in-silico study proposed B3 domain factor LOC_Os08g06120 to colocalize with a QTL for seed width (Peng and Weselake, 2013). The PHD Zn finger gene OsMYB7 (LOC_Os08g43550) and the CCT motif orphan gene (LOC_Os02g05470) have no defined function but were both downregulated under desiccation stress in IR64 (Shankar et al., 2016). OsWRKY19 (LOC_Os05g49620) was shown to be constitutively low expressed across different tissues and was upregulated by IAA and SA with no clear involvement in any stress response (Ramamoorthy et al., 2008).

SUB1A promoter analysis
A comparative sequence analysis of the promoter region of the SUB1A-1 and SUB1A-2 alleles (Additional file 2: Table S5) revealed the presence of 13 SNPs. Out of these, 5 SNPs gave rise to allele-specific putative cisregulatory elements that could be targets of upstream TFs. (Fig. 5; Additional file 3: Figure S2).
In the SUB1A-1 allele, SNP4, SNP5, and SNP13 create putative motifs (Fig. 5). SNP4 constitutes a Site I (SITEIOSPCNA) motif that resembles the G-box and was discovered in the rice PROLIFERATING CELL NUCLEAR ANTIGEN (PCNA) promoter. Promoter deletion studies had shown that the Site I element is a transcriptional activation site of PCNA and that it is essential for meristematic tissue-specific expression (Kosugi, Suzuka, and Ohashi, 1995). SNP5 resulted in two distinct cis-elements in the two SUB1A alleles (see below) creating a CAREOS-REP1 motif in SUB1A-1. This motif was shown to confer GA-mediated increased expression of a proteinase gene in rice seeds (Sutoh and Yamauchi, 2003). SNP13 creates the other constituent of the Site I bipartite motif, i.e., the Site II element (SITEIIATCYTC). This motif was also formed by SNP9 in the SUB1A-2 allele (see below). In addition to PCNA, Site II elements are present in the promoters of many Arabidopsis cell cycle, mitochondrial, ribosomal and cytochrome genes (Tremousaygue et al., 2003;Welchen and Gonzalez, 2005). They are also overrepresented in the promoters of genes encoding nuclear components of the oxidative phosphorylation complex in Arabidopsis and rice (Welchen and Gonzalez, 2005;Welchen and Gonzalez, 2006). Promoter::GUS studies showed that the Site II motif is sufficient to mediate gene expression in shoot and root meristems, as well as in anthers (Kosugi and Ohashi, 1997). It is recognised by TCP and bHLH TFs, which play crucial roles in hormone response pathways, light and abiotic stress responses (Danisman, 2016). Importantly, the position of the Site II element specific to the SUB1A-1 allele is closer to the transcription start site (TSS), which was shown to be the more active position, thus makes this element more relevant in SUB1A-1 allele compared to SUB1A-2, where it is located further upstream (Tremousaygue et al., 2003;Welchen and Gonzalez, 2006).
A third Site II element is present in both alleles (Fig. 5). In addition to Site I and Site II elements, a telo-box (5'-AAACCCTAA-3′) was identified in the PCNA promoters and shown to act as a transcriptional activator (Tremousaygue et al., 2003). In the SUB1A alleles, no perfect match for the telo-box could be found, however, two truncated telo-box motifs (5'-ACCCTA-3′ and 5'-AAACCCT-3′) are present near the SUB1A-2 and SUB1A-1 specific Site II elements, respectively (Fig. 5).

Promoter analysis of the DEGs
The analysis of the 2Kb upstream regions of the DEGs identified 32 distinct motifs present in at least one of the genes (Table 3). These motifs can be categorized into four groups based on their relatedness to hypoxia & light, energy & sugar metabolism, transcription & cell division, and hormones & stress.
Under submergence, plants are exposed to hypoxia and low light conditions and in agreement with this, four anaerobic consensus motifs (ANAERO1-4) and the light response motif GT1CONC4:C36SENSUS were found in most DEGs. The anaerobic consensus motifs were previously identified in genes involved in the anaerobic fermentative pathway (Mohanty, Krishnan, Swarup, and Bajic, 2005). The GT1 motif has been found in light-regulated genes and is known to stabilize the TATA box complex (Terzaghi and Cashmore, 1995;Villain, Fig. 5 Comparison between the SUB1A-1 and SUB1A-2 upstream region. The SUB1A-1 and SUB1A-2 upstream regions were aligned and analysed to identify the promoter region and allele-specific putative cis-elements constituted by the thirteen single nucleotide polymorphisms (SNPs). Annotation of the putative cis-elements is indicated in the figure and explained in detail in the text. Other relevant promoter elements are additionally indicated. TSS, transcription start site. The complete alignment is provided in Additional file 1: Figure S1 Mache, and Zhou, 1996;Buchel, Brederode, Bol, and Linthorst, 1999;Le Gourrierec, Li, and Zhou, 1999;Zhou, 1999).
Sixteen of the identified motifs were related to energy and sugar metabolism, and studies associated with some of these motifs are detailed below. RAMY motifs, identified in alpha-amylase genes, were highly represented in all four groups of DEGs (up and downregulated DEG in IR64 and IR64-Sub1). The RAMY3D motif (TATCCAY) was found in the rice RAmy3D alpha-amylase gene promoter and is responsible for sugar repression (Toyofuku, T-a, and Yamaguchi, 1998). The AMY3 motif (CGACGO) is present in the GC-rich regions of amylase genes and functions as a coupling element for the G-box element (Hwang, Karrer, Thomas, Chen, and Rodriguez, 1998). The other AMY motif (TATCCA) is a binding site for OsMYBS1, OsMYBS2 and OsMYBS3 TFs, which mediate sugar and hormone regulation of alpha-amylase gene expression (Lu, Lim, and Yu, 1998;Lu, Ho, Ho, and Yu, 2002;Chen, Chiang, Tseng, and Yu, 2006). GARE and CARE motifs are necessary for GA regulation of alphaamylases genes (Sutoh and Yamauchi, 2003) and the BP5OSWX motif acts as a transcriptional activator in the waxy gene (Zhu, Cai, Wang, and Hong, 2003). The GARE motif was present in about 50% of the IR64-Sub1 upregulated genes and less frequent (15%) in the downregulated genes in both genotypes (Table 3). The GARE and RAMY1A (PYRIMIDINEBOX) motifs participate in sugar repression and are present in the promoter of the rice alpha-amylase gene RAmy1A and the barley alphaamylase gene Amy2/32b, both of which are GA inducible (Morita et al., 1998;Mena, Cejudo, Isabel-Lamoneda, and Carbonero, 2002). The four identified GLUB1 motifs have been associated with the rice glutelin gene OsGluB-1.
They are related to energy metabolism and are proposed to control endosperm-specific expression (Washida et al., 1999;Wu, Washida, Onodera, Harada, and Takaiwa, 2000;Qu, Xing, Liu, Xu, and Song, 2008). Glutelins are primary energy storage proteins and the relevance of these motifs remains to be assessed with respect to submergence response. Site II motifs as detailed above, were also found in most DEGs. Seven of the identified motifs have been related to transcription (E2F, EF1, and TATA Box), meristems (LEAFYA-TAG and WUSATAg) and homeotic development (BIHD1OS). The E2F binding site was present in 80% of the IR64-Sub1 upregulated genes and less frequent in the downregulated genes (IR64-50% and IR64-Sub1-10%). E2F TFs are regulators of the cell cycle and can act as transcriptional activators or repressors (Vandepoele et al., 2005). The meristem motifs were rare compared to BIH-D1OS, which was present in most DEGs. The BIHD1OS motif has been reported to be the binding site for OsBIHD1, a rice BELL homeodomain TF that integrates ethylene and BR hormonal responses (Liu et al., 2017).
Four of the DEG motifs were associated with hormone and stress responses. The drought-related DRE/CRT and ABA-responsive ABRE motifs were most frequent (Dubouzet et al., 2003;Qin et al., 2004;Diaz-Martin et al., 2005;Skinner et al., 2005;Suzuki et al., 2005), whereas the Fe starvation related motif IRO2 (Ogo et al., 2006) was present in only about 10% of the genes. Amongst the groups of DEGs, DRE/CRT and ABRE motifs were most common in the IR64-Sub1 downregulated genes (60-80%) compared to the other two groups (Table 3).

Discussion
All SUB1A-related studies so far have been conducted using contrasting genotypes with and without the SUB1A gene and, therefore, limited information is available on the differential effect of the two alleles. As a first step, the objective of this study was therefore to assess if different sets of TF genes are regulated by the strong SUB1A-1 allele, as present in IR64-Sub1, and the weaker SUB1A-2 allele, naturally present in the wild-type IR64. An earlier study on a range of rice genotypes had shown that both alleles are highly expressed in nodes of submerged plants, however, differences in the expression level in internodes were most pronounced and correlated best with submergence tolerance (Singh et al., 2010). Therefore, internode tissue was used for this comparative study.
The qPCR analysis of the 2487 TF genes revealed a relatively small number (46) of TF genes that were significantly differentially expressed under submergence. Interestingly, the predominant transcriptional response was reduced expression and only thirteen genes, all in IR64-Sub1, were upregulated under submergence. This, per se, might at least partly explain the higher tolerance level of IR64-Sub1. Further to this, with the exception of an uncharacterized MYB factor, distinct sets of TFs were identified in IR64 and IR64-Sub1, suggesting that the two genotypes employ different regulatory pathways leading to different levels of tolerance.

Association between submergence responsive DEGs and plant hormones
Nineteen TF genes were specifically submergenceresponsive in IR64-Sub1, making them candidates to explain the higher tolerance of IR64-Sub1 compared to IR64. The upregulated genes were previously reported in studies on drought, salt and other abiotic stresses, but not with respect to submergence stress (Additional file 2: Table S7). In agreement with our data, the OsbZIP4 and the Zn-finger gene LOC_Os07g43740 were found to be upregulated under anoxia (Ham et al., 2013;Shankar et al., 2016).
Available data on the hormonal response of some of these genes (OsHsfC2a, OsBBX11, OsBBX24, OsWRKY76 and OsTIFY11a) suggest that they are part of the ABA/ GA regulatory pathway (Fig. 4). For instance, one of the six upregulated MYB factors (LOC_Os08g06110) was reported to reduce seed dormancy in a rice mutant through ABA-GA antagonism . OsHsfC2a is downregulated by ABA, while, GA causes downregulation of OsBBX11 and OsBBX24 (Huang et al., 2012). Additionally, OsBBX11 is regulated by light, which plays a critical role in the submergence response (Das et al., 2009).
Among the genes downregulated in IR64-Sub1 were two AP2-type TFs (OsERF136, OsDREB2b), both of which are submergence-responsive in M202 and M202-Sub1 (Jung et al., 2010). Since the wild-type genotype M202 naturally lacks SUB1A, the two AP2 genes might act in a different pathway, independent of SUB1A. DREB2b has also been reported to function as a transcriptional activator under low-and high-temperature stress and it enhances drought tolerance when overexpressed (Matsukura et al., 2010;Todaka et al., 2012). Similarly, SUB1A has been implicated with enhanced recovery growth after severe drought (Fukao et al., 2011) and submerged plants generally survive longer in cooler water (Das et al., 2009). These similarities between SUB1A and DREB2b make the latter an interesting candidate gene for integrating these signal responses.
In contrast to IR64-Sub1, the DEGs in IR64 were all downregulated (Fig. 5). The most highly downregulated gene encodes an uncharacterized protein, with putative Jas-and PHD Zn-finger domains, suggesting that it may be a part of the JA signalling pathway and/or a regulator of development. The second most highly responsive gene was OsSUB1B, followed by WRKY and additional ERF/AP2 factors, including OsSUB1C. Interestingly, both OsSUB1B and OsSUB1C are, despite also in IR64-Sub1, only differentially expressed (downregulated) in IR64.
OsSUB1B is non-responsive to GA and ethylene and no tolerant-specific allele has been identified for this gene . A clear downregulation of OsSUB1B under submergence in IR64 justifies more in-depth analysis of its function and regulation of expression. In contrast, tolerance-specific alleles for OsSUB1C are known and its induction by GA and ethylene specifically in M202-Sub1 has been documented Xu et al., 2006). These studies also indicated that SUB1A-1 might be a negative regulator of OsSUB1C. Downregulation of OsSUB1C in IR64 suggests that the OsSUB1A-2 allele can likewise act as a negative regulator of OsSUB1C, directly or indirectly.
The three AP2/ERF genes downregulated in IR64 (OsERF067, OsERF066, OsERF922) had previously been shown to be upregulated in M202-Sub1 under submergence (Jung et al., 2010). This was not observed in IR64-Sub1 in this study, however, their downregulation in IR64 might at least partly explain a lower level of submergence tolerance. OsERF922 may be of particular interest here since this gene acts as a transcriptional activator and modulator of ABA levels (Liu et al., 2012).
The BBX1 gene had previously been shown to be submergence induced in the variety M202, which lacks the SUB1A gene (Fukao et al., 2011). This is in contrast to the observed downregulation of BBX1 in IR64 in our study, suggesting that BBX1 might be under the control of the SUB1A-2 allele contributing to the higher tolerance level of IR64 compared with M202.
As was the case for IR64-Sub1, a number of the IR64 downregulated DEGs have been associated with GA and ABA regulation. Among the GA-responsive genes, the MYB TF (LOC_Os03g04900) might be of particular interest as it is similar to OsMPS, which acts as a regulator of plant growth by suppressing expansin genes and is downregulated by GA and upregulated by ABA treatment (Schmidt et al., 2013). Suppression of expansins has previously been shown to be a critical component of the SUB1A-mediated quiescence response and it will be interesting to assess if this MYB factor functions in a way similar to OsMPS.
Genes upregulated by ABA (OsMYB2, OsWRKY90, OsERF922 and OsMYB48-1) activate the ABA synthesis pathway and OsMYB48-1, specifically, has been reported to induce other ABA-responsive genes under drought (Xiong et al., 2014). On the other hand, four genes (OsSUB1C, OsBBX1, OsBBX30, OsWRKY62) have antagonistic effects to ABA and respond positively to GA (Ramamoorthy et al., 2008;Huang et al., 2012;Shankar et al., 2016). The fact that all of these genes are downregulated in IR64 under submergence suggests that at least part of the GA signalling pathway is inactivated, as required for submergence tolerance. That these genes have not identified as differentially expressed in IR64-Sub1 is not necessarily a contradiction but instead suggests that submergence tolerance might be related to low or neutral expression of these genes.

Role of iron uptake and JA in submergence tolerance
Iron solubility and availability in plants is altered based on soil pH and submergence conditions and can result in excess iron availability. Difficulties in iron uptake cause chlorosis whereas iron overload or increased uptake leads to oxidative stress and permanent cell and tissue damage (Zheng 2010). One of the four TFs that showed a constitutive difference between the two genotypes, OsIRO2, was constitutively lower expressed in IR64-Sub1 (Fig. 2a). IRO2 has so far only been described in relation to Fe and it was shown to be specifically induced by Fe deficiency (Ogo et al., 2006;Kobayashi et al., 2007). Although the specific role during submergence remains to be determined, the higher expression of the IRO2 in IR64 might suggest a higher Fe uptake in this cultivar resulting in oxidative stress and cell damage.
As expression of IRO2 is not induced by submergence, a specific role for IRO2 is perhaps unlikely, however, it might be possible that the SUB1A-2 allele specifically acts as an enhancer of OsIRO2 thereby causing a secondary, negative effect under submergence. It will, therefore, be interesting to compare the performance of IR64 and IR64-Sub1 under Fe deficiency conditions and quantify plant iron content.
Five of the other TF genes identified in this study (OsTIFY11e, OsTIFY11d, OsWRKY76, OsTIFY11a, OsWRKY71) were reported to be upregulated by JA and downregulated by ABA (Ye et al., 2009;Ranjan et al., 2015;Yang et al., 2015;Shankar et al., 2016). Additionally, WRKY76 and WRKY71 are also upregulated by GA (Ramamoorthy et al., 2008;Shankar et al., 2016). This is well in support of the finding that fine-tuning of GA and ABA levels under submergence is a major component of regulating elongation growth in rice Fukao and Bailey-Serres, 2008) and studies on Arabidopsis have further shown that JA plays a crucial role in protecting plants during post-submergence reoxygenation (Yuan et al., 2017). The genes identified in this study might therefore play a role in integrating various hormonal signals and analysing their expression in IR64 and IR64-Sub1 plants upon de-submergence will provide more insights into the role of SUB1A-2 during recovery growth.
A significant number of DEGs controlled by submergence have also been described as cold stress responsive. In addition, it has been shown that plant survival in cooler submergence water is higher (Das et al., 2009). Therefore, higher expression of these cold-induced genes in IR64-Sub1 may contribute to its higher tolerance level. For the genes listed in Additional file 2: Table S7, cold responsiveness was also observed in M202 and M202-Sub1 (Jung et al., 2010).
Putative cis-regulatory elements The diversity of the putative regulatory cis-elements identified in the promoters of the DEGs (Table 3) is well in agreement with the diverse environmental stimuli associated with submergence requiring plants to adjust to low oxygen, light and temperature, alter their energy and sugar metabolism, cell divisions, as well as hormone-regulated growth. For instance, the light-responsive motif (GT1CONC4:C36SENSUS) occurred in all genes except one, suggesting that low-light is an integral signal that can trigger both, up-and downregulation of genes. Likewise, the majority of genes appears to be responsive to hypoxia, as indicated by the presence of the anaerobic consensus motifs, ANAERO1CONSENSUS and ANAERO2CON-SENSUS, as well as the PYRIMIDINEBOXOSRAMY1A, found in alpha-amylase genes. The AMY motif (TATCCA) appears to be more frequent in IR64-Sub1 downregulated genes, whereas the AMY3 motifs (CGACG), which functions as a G-box coupling element, is more prevalent in upregulated genes. The former is the binding site of OsMYBs, which mediate sugar and hormone regulation of alpha-amylase genes (Lu et al., 1998;Lu et al., 2002;Chen et al., 2006). Since maintenance of a higher level of soluble sugars has been implicated with tolerance and better recovery growth (Singh et al., 2001;Das et al., 2005; these TF genes might be part of the regulatory pathway controlling sugar consumption under submergence. Despite the importance of ABA in submergence tolerance, surprisingly few genes had ABRE elements. In contrast, the DRE/CRT element was more frequent, which, in conjunction with the identified differential expression of DREB2b (Matsukura et al., 2010;Todaka et al., 2012), reinforces a possible relationship between submergence and drought tolerance (Fukao et al., 2011).
Among the most interesting elements identified in the promoter analysis is the bi-partite Site II element. This motif was frequently present in the DEGs and also distinguished the SUB1A-1 and SUB1A-2 alleles at the SNPs (Fig. 5). Since this motif has been implicated with meristem-specific expression (Tremousaygue et al., 2003;Welchen and Gonzalez, 2006) its presence might indicate that control of cell division, in addition to elongation growth, is important for tolerance. This is supported by promoter::GUS reporter studies that showed specific SUB1A-1 expression in the base of growing leaves (Singh et al., 2010). It will be interesting to study this motif in more detail and to determine if the position closer to the TTS of the SUB1A-1-specific Site II element enhances the capacity of tolerant plants to suppress cell division under submergence. In relation to this, it is relevant that the two SUB1A alleles are also distinct in their coding region for a MAPK3 phosphorylation target site, which is specifically present in SUB1A-1 (Singh and Sinha, 2016). The same study also showed that both, SUB1A-1 and SUB1A-2, physically interact with MAPK3 but that the interaction with SUB1A-2 was weaker. Thus, SUB1A-2 would still be a component of this regulatory pathway but might be less efficient. It will indeed be interesting to investigate the role of SUB1A in regulating cell division in more detail and to establish whether SUB1A-2 binds to the MAPK3 promoter, as was shown for SUB1A-1 (Singh and Sinha, 2016).

Conclusion
This study identified distinct sets of submergenceresponsive transcription factor genes in IR64 and IR64-Sub1 which might at least partly explain differences in the level of tolerance between the two genotypes. These genes should now be assessed for putative additive effects with SUB1A and their potential to enhance tolerance. The promoter analysis identified allele-specific site II elements, which should be studied in more detail to establish if SUB1A directly or indirectly suppresses cell division, in addition to its documented role in suppressing cell elongation.

Additional files
Additional file 1: Figure S1. The submergence tolerance locus Sub1 in rice. The Sub1 locus is located on rice chromosome 9 with a variable number of SUB1 genes. In different rice genotypes, SUB1A is either absent or present as different allele. (a). Introgression of the Sub1 locus into the rice variety IR64 (IR64-Sub1) enhanced survival after complete submergence. (b). The photo shows the IRRI demonstration field plot in the Philippines. Phenotyping for submergence tolerance can also be conducted by submerging plants grown in trays (inlay). The SUB1A-1 and SUB1A-2 alleles are both highly expressed in nodes of submerged plants but SUB1A-1 expression is higher in internodes. (c). Schematic illustration based on (Singh et al., 2010). (PPTX 3066 kb) Additional file 2: Table S1. List of genes represented in Rice TF Primer Platform (Caldana et al., 2007). Complete list of 2508 TF genes, their primer sequences and their corresponding PCR efficiencies. Genes were presented with version 2.0 and version 5.0 RGAP genome annotations. Table S2. Transcription factor genes with and without annotations in version 5.0 and 7.0 of RGAP Pseudomolecule. Table S3. Expression profile of transcription factor genes in IR64 control and stress-treated plants upon submergence. Recorded Ct, R 2 , Efficiency and ΔCt values for all reactions. The list of used reference genes and their distribution on the plates is also given. Table S4. Expression profile of transcription factor genes in IR64-Sub1 control and stress-treated plants upon submergence. Worksheets record Ct, R 2 , Efficiency and ΔCt values for all reactions. The list of used reference genes and their distribution on the plates is also given. Table S5. Promoter sequences of SUB1A-1 and SUB1A-2. Table S6. Moderated t-tests results from LIMMA package for identifying differentially expressed genes. P-values were recorded and FDR corrections were also performed. Table S7. Details about the list of differentially expressed genes and their putative functions. (DOCX 29 kb) Additional file 3; Figure S2. Sequence alignment of the SUB1A-1 and SUB1A-2 upstream promoter regions. 2 kb upstream of the start site (promoter regions) were analysed for SUB1A-1 and SUB1A-2 alleles and putative cis-regulatory elements were identified. (XLSX 3.76 kb)