Six rice genotypes were used in this study: (1) The japonica variety Nipponbare; (2) The indica variety IR64; (3) Two recombinant inbred lines (RILs; F291 and F274-2a) derived from a cross between Nipponbare and IR64, and; (4) Two accessions that originate from southeast Asia: 8391 from Laos (IRGC 94599) and 8753 from Indonesia (IRGC 54313).

All of the seeds used in this study were freshly harvested and stored at 4°C prior to experiments; sterilized seeds were germinated in capped glass tubes in a growth chamber at 25°C for 7 days under a photoperiod comprising a cycle of 16 h of light (150 μmole m⁻² s⁻¹) and 8 h of dark, designated the “control” condition. For “submergence” treatments, sterilized seeds were germinated in 5 cm of water for 7 days.

Coleoptile length was measured using an ordinary ruler; however, considering that coleoptiles in the control group grown in air were barely elongated, the entire aboveground shoot (i.e., tissue including the coleoptile and primary leaf) of 7-day-old seedlings that showed a consistent response to both air and submergence was harvested and stored at −80°C. Three independent experiments including 15 plants were performed for each accession, and the plant tissues from three independent biological replicates were pooled for RNA isolation.

Total RNA was extracted using TRI Reagent® (Invitrogen, MA, USA) and a Direct-zol™ RNA MiniPrep Kit (ZymoResearch, CA, USA). RNA samples with an RQI (RNA quality indicator) greater than 7.5 were used for library preparation. Poly-A RNA containing mRNA was purified using poly-T oligo-attached magnetic beads and fragmented, and complementary DNA (cDNA) was synthesized using random hexamer primers, followed by purification, end-repairing, poly-A tailing, and adaptor ligation.

Twelve cDNA libraries comprising unique barcodes were pooled and sequenced in one run using an Illumina HiSeq™ 2500 to generate single-end reads, each 100 base pairs (bp) in length.

Raw sequenced short reads were aligned to the rice reference genome (MSU7.0) using the CLC Genomics Workbench 6.5 software (CLCbio-Qiagen, Aarhus, Denmark). Expression levels of each gene were quantified by normalizing total exon read counts with effective library size, and tests for pairwise differential expression were performed using the R software package DESeq in Bioconductor (Anders and Huber, 2010). Genes with P-values less than 0.05 were considered differentially expressed (DE). These were then used as the basis for biological pathway annotation in the MapMan software (ver. 3.5.1R2) (Usadel et al., 2005). Osa_MSU_v7 mapping files were downloaded from the MapManStore server (http://mapman.gabipd.org/web/guest/mapmanstore).

Applying multi-factorial linear modeling, we were then able to test three null hypotheses of effects on gene expression: (1) whether the expression of each gene was significantly regulated by submergence treatment; (2) whether gene expression was affected by genotype, and; (3) whether gene expression was affected by submergence in a genotype-dependent manner. We therefore fitted our experimental data into four different linear models: (1) FM_trt: Y = τ + ε; (2) FM_geno: Y = τ + ε; (3) FM_add: Y = τ + γ + ε, and; (4) FM_full: Y = τ + γ + τ:γ + ε. In each of these formulae, Y is the expression value of each gene, τ is the effect of submergence treatment, Y is the effect of different genotypes, and ε is random error. Comparing FM_geno and FM_trt to FM_add separately, we then tested whether the expression of each gene was regulated by submergence, and whether there was a significant genotypic effect. Comparisons of FM_full and FM_add allowed us to test whether gene expression was affected by submergence in a genotype-dependent manner.

In all cases, expression values for genes were standardized using the expression z=(x-x-)/sx for cross-genotype comparisons.

Total RNA from 7-day-old aboveground shoot tissues was extracted using TRI Reagent® (Invitrogen, MA, USA) with a Direct-zol™ RNA MiniPrep Kit (ZymoResearch, CA, USA). We included three independent biological replicates in this experiment, and designed gene-specific primers using the NCBI primer BLAST. The complete primer sets used in this study are listed in Table S1, and quantitative RT-PCR was performed using the Rotor-Gene® SYBR® Green RT-PCR kit (Qiagen, Hilden, Germany), following the manufacturer's protocol for one-step RT-PCR. The PCR conditions used in this study encompassed 10 min at 55°C for reverse transcription, 5 min of pre-denaturation at 95°C, 40 cycles of 5 s at 95°C and 10 s at 60°C, followed by melting curve generation. We selected the OsACT1 gene (forward primer: 5′-ATGAAGATCAAGGTGGTCGC-3′; reverse primer: 5′-GTACTCAGCCTTGGCAATCC-3′) as the internal reference for relative quantification, and used the delta-delta C_T method to calculate relative expression level fold-changes between submerged and control samples (Pfaffl, 2001).

Log₂-transformed relative expression levels of different samples were used to evaluate correlation in expression patterns determined by RNA-Seq and qRT-PCR for each selected gene.

We performed GO enrichment analysis on DE genes using the web-based toolkit “AgriGO” (Du et al., 2010). We used Fisher's exact test and the Benjamini-Hochberg's false discovery rate (FDR) adjustment to control for multiple comparisons (Benjamini and Hochberg, 1995).

All the sequence data generated in this research was deposited in the Experiment ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/) at EMBL-EBI (European Molecular Biology Laboratory-European Bioinformatics Institute) under accession number: E-MTAB-3834.

We examined the length of 7-day-old rice seedling coleoptiles across the six rice genotypes (i.e., Nipponbare, IR64, F291, F274-2a, 8391, and 8753) either under normal conditions or submerged (Figure 1, Table 1). The results of Fisher's least significant difference (LSD) tests demonstrated significant variation in coleoptile length when submerged seedlings of different genotypes were compared. Although variation in coleoptile growth between the rice genotypes was also observed under normal conditions, this was not as significant as in submerged examples. Thus, we used coleoptile length differences between submerged and control samples to evaluate the ability of seedlings to tolerate anaerobic stress induced by waterlogging. Results show that the anaerobic response of Nipponbare genotype seedlings was significantly better than the response of the IR64 genotype, and that two Nipponbare/IR64-derived RILs, F291 and F274-2a, both outperformed the former. We have previously confirmed the presence of this significant transgressive variation in an earlier mapping study (Hsu and Tung, 2015), suggesting that epistatic interactions between parental genomes could contribute to anaerobic tolerance. We have also previously identified two extremely tolerant rice landrace genotypes, 8391 and 8753; under submergence the coleoptile length in both these landraces was twice the length of that seen in genotype IR64 seedlings. According to their coleoptile elongation ability, we designated IR64 as sensitive genotype, Nipponbare and others are moderate and extreme tolerant genotypes (Figure 1).

We constructed 12 cDNA libraries using total RNA purified from the tissue of six genotypes of 7-day-old control (i.e., germinated in air) and submerged aboveground seedlings tissue. We generated a total of 403.65 million short reads (100 bp) from 12 barcoded cDNA libraries in one run for this study using a HiSeq 2500 machine. The detailed sequencing statistics for these cDNA libraries are presented in Table S2, and quality-trimmed reads were aligned to the rice reference genome (i.e., Nipponbare, MSU7.0 annotation). More than 90% of these reads were mapped, and between 87 and 90% were uniquely aligned (Table S3). Total counts of exon-aligned reads were normalized on the basis of effective library size for each sample, using the DESeq software (Anders and Huber, 2010), and the expression signal of each gene was calculated (see above). To further validate transcription expression generated by RNA-seq and DESeq, we selected eight genes based on their read counts and examined their expression levels using RT-PCR. The results of this analysis demonstrated that both quantification platforms generated similar patterns; Pearson's correlation coefficients (r) for these eight genes range between 0.63 and 0.96, with an average of 0.79 (Figure S1, Table S4).

In order to overview the genome-wide transcriptomic changes between control and submerged treatment samples irrespective of their genotypic identities, we averaged the expression signal of each gene across the six genotypes in the control and treatment groups, respectively, and then visualized overall transcriptomic change using the MapMan software (Usadel et al., 2005) to identify metabolic pathways and biological processes regulated under submergence. The results of these comparisons show that genes associated with cell wall modification, fermentation, transcription regulation, and hormone biosynthesis were all up-regulated in the submerged group, while those involved in cell wall degradation, mitochondrial electron transport, photosynthesis, and secondary metabolism were mostly down-regulated (Figure 2).

A submerged rice shoot is mainly composed of elongated coleoptile (Figure 1), a specialized organ that does not produce chlorophyll during anaerobic germination. Our results confirm the expression of photosynthetic genes in control green shoots, while photosynthetic activity was reduced in the submerged coleoptiles (the “light reactions” illustrated in Figure 2). Previous results have shown that the coleoptile is the first organ to senesce when rice seedlings are germinated in air (Inada et al., 2000; Kawai and Uchimiya, 2000). Our results further demonstrate that submergence delays the processes of senescence and cell death, as a significant number of genes were actively induced when seedlings were germinated under water. These genes triggered a series of downstream transcriptional regulators that promote adequate metabolic and morphological adjustments to cope with submergence stress.

To further examine whether the gene expression was affected by (1) submergence treatment only, (2) the genotype only, or (3) both genotype and submergence, we performed multi-factor linear model testing on our whole transcriptome dataset using DESeq. Applying a P-value of less than 0.05 as the cutoff, our results show that 3,597 genes were significantly affected by submergence regardless of their genetic backgrounds (i.e., the “treatment effect” test in Table S5), while 5,100 genes were differentially expressed across six genotypes irrespective of submergence (i.e., the “genotype effect” test in Table S5), and the expression of 471 genes was affected by submergence in a genotype-dependent manner (i.e., the “interaction” test in Table S5). These results therefore provide the first confirmation that gene expression under submergence can be regulated differently depending on genetic background; in other words, variation of coleoptile elongation in diverse rice accessions could potentially be determined by genotype-specific gene expression. We performed GO enrichment analysis separately for up-regulated and down-regulated gene sets of the 3,597 submergence-responsive genes from linear modeling by the AgriGO software (Du et al., 2010). The results of this analysis reveal a similar pattern to MapMan visualization (Figure 2)—that genes responsible for oxidative stress responses, transcriptional regulation, post-translational regulation, and cell wall organization are specifically enriched in the up-regulated gene set, while those responsible for carbon fixation, carboxylic acid metabolism, and a number of other macromolecular metabolic processes are all down-regulated (Table S6).

In order to investigate in detail the distribution of DE genes in our six diverse transcriptomic profiles, we next analyzed gene expression in submerged and control samples genotype-by-genotype independently using negative binomial testing in DESeq, and then applied a P-value less than 0.05 as a cutoff to generate a DE gene list for each genotype. Using this approach, we identified 701 and 808 DE genes in the two extremely tolerant landrace genotypes 8391 and 8753, total 692 genes in the moderately tolerant cultivar Nipponbare genotype, and just 193 in the sensitive cultivar genotype IR64 (Figure S2). According to the MA plots (Figure S3), all DE genes showed significant differential expression (log₂ Fold Change > 2 or < −2) and strong expression (mean expression > 100); this suggests that a P-value < 0.05 is a reasonable threshold in this study to identify differential expressed genes. In addition, 83 DE genes (i.e., 68 up-regulated and 15 down-regulated) were commonly present in four rice varieties (Nipponbare, IR64, 8391, 8753), while a number of genotype-specific DE genes can also be observed (Figure S2, Table S7). Our findings suggest that rice coleoptile elongation when submerged at the germination stage has been fine-tuned based on differences in gene expression pattern in different genetic backgrounds.

In the two RILs derived from the Nipponbare and IR64 genotypes, we identified 527 and 929 DE genes in F291 and F274-2a, respectively, while 93 genes were uniquely expressed in recombinants. These data imply that unique expression of these genes contributes to the transgressive phenotype and leads to rapid coleoptile growth (Figure S2, Table S7).

As revealed by MapMan analysis (Figure S4), because transcriptomic activity was significantly reduced in the sensitive IR64 cultivar comparing to other tolerant genotypes, it is possible that poor coleoptile elongation in this case could be attributed to weak transcriptomic changes. Unlike the other five genotypes, IR64 exhibited poor transcriptomic responses under submergence, both in terms of the number of responsive genes and the level of responsiveness (Figure 3 and Figure S5). Thus, based on the results of pairwise DE analyses, we identified sets of genes of interest according to the following criteria: (1) Genes that were responsible for fundamental submergence responses irrespective of genotype; (2) Genes that were differentially expressed in the five tolerant lines but that responded poorly in the sensitive IR64 genotype; (3) Genes that promoted coleoptile elongation in the four rapidly-growing genotypes, and; (4) Genes that exhibited an absence of transcripts in the sensitive IR64 genotype.

Several previous studies have investigated the expression of target genes and the biological pathways involved in coleoptile growth during anaerobic germination (Magneschi et al., 2009; Miro and Ismail, 2013; Lee et al., 2014). In this study, we confirm the regulation of previously reported submergence responsive pathways in tolerant genotypes, although some modulation defects were found in the sensitive cultivar IR64, accounting for its poor seedling growth performance when submerged. Additional genes potentially involved in cell wall biosynthesis or modification pathways were strongly regulated in submerged seedlings of highly tolerant genotypes, promoting faster coleoptile elongation to escape submergence.

In order to comprehensively discuss the molecular basis of coleoptile growth variation across our six diverse genotypes, we selected three major biological processes that involve large numbers of the candidate genes identified in this study and that have been studied intensively in the context of rice anaerobic germination (Figure 5, Figure S6). We compared and discussed the expression patterns of these target genes across our six genotypes.

Rice seedlings take up water in order to germinate; during the first 48 h of this process, metabolic activity is reactivated and starch in the endosperm is degraded, followed by hydrolysis and glycolysis. The energy generated by these biochemical processes is then used for embryonic tissue growth, and the first organ to protrude is the coleoptile. Elongation of this structure requires energy for cell growth and maintenance; when O₂ is limited, carbohydrate catabolism shifts to anaerobic respiration via pyruvate metabolism, correlated with an increase in the expression of genes involved in fermentative pathways (Lasanthi-Kudahettige et al., 2007; Narsai et al., 2009).

Building on these previous observations, we investigated the genes involved in starch and sucrose mobilization across our six genotypes, including those encoding α-amylase, sucrose synthase (SuSy), ADP-Glc PPase, UTP-Glc-1-P uridylyltransferase, phosphoglucomutase (PGM), and hexokinase (Table S9). Two of these enzymes in particular are considered essential for carbohydrate catabolism in rice seedling germination, α-amylase for starch degradation and sucrose synthase for sucrose breakdown. Previous work has shown that RAmy3D expression can be induced in response to an O₂ deficit during early germination (Guglielminetti et al., 1995a,b; Perata et al., 1997; Hwang et al., 1998, 1999), and that expression between 24 and 72 h after germination is positively correlated with shoot growth (Ismail et al., 2009). Our data show that of the ten α-amylase encoding genes we analyzed, only RAmy3D (LOC_Os08g36910) was up-regulated in the five tolerant genotypes and was barely detectable in IR64. At the same time, correlation between the induction level of this gene and coleoptile elongation rate was not linear in our tolerant varieties, an observation that could be explained by our late sample collection time. We observed a similar regulation pattern to that reported in previous work (Magneschi and Perata, 2009) in the case of the Nipponbare genotype; the activity of glycolytic genes, with the exception of those encoding sucrose synthase 1 (SuSy1), hexokinase 7 (HXK7), and PGM, were accelerated in the anoxic embryo and coleoptile tissue (Lasanthi-Kudahettige et al., 2007; Narsai et al., 2009). Indeed, as previously described, the genes involved in glycolytic processes, including those encoding G3PDH (LOC_Os02g38920) and PPDK (LOC_Os03g31750), were all significantly induced, albeit at different levels, in the five tolerant genotypes but not in the sensitive IR64 (Figure 5). This differential regulation could control the rate of glycolysis and influence pyruvate production for fermentative pathways, resulting in variation in coleoptile growth.

Another example of differential modulation of the carbohydrate metabolism pathway is expression of the gene encoding phosphoenolpyruvate carboxylase (PEPC, LOC_Os08g27840). Because strong suppression of this gene was detected in tolerant lines and a lower level of reduction was observed in the IR64 genotype (Figure 5), we hypothesize that the presence of PEPC in the submerged IR64 coleoptile might result in insufficient pyruvate for ethanol production. Thus, because the gene encoding PEPC is significantly expressed in aerobic tissue; when it is unrepressed it will convert PEP into oxaloacetate (Figure 5), bypassing alcohol fermentation.

Previous work has also demonstrated that the activities of enzymes involved in fermentative ethanol production, including pyruvate decarboxylase (PDC) and alcohol dehydrogenase (ADH), are significantly higher in seeds germinating under hypoxia in the tolerant “Khaiyan” rice variety than they are in the sensitive “IR42” variety (Ismail et al., 2009). We detected significant induction of genes encoding both PDC (LOC_Os05g39310) and ADH (LOC_Os11g10480 and LOC_Os11g10510) in the five tolerant genotypes but not in IR64 (Figure 5), suggesting that ethanol production could be impaired in the sensitive genotype via transcriptional regulation, leading to insufficient energy generation for coleoptile growth.

Several previous studies have suggested that because ATP-consuming processes are not favored by hypoxia, PPi should be considered as an alternative energy donor to potentially maintain cell growth (Gibbs and Greenway, 2003; Magneschi and Perata, 2009; Atwell et al., 2015). Thus, to elucidate the contribution of ATP- and PPi-dependent enzymes in anaerobic germination, we investigated the expression patterns of genes encoding ATP-dependent phosphofructokinases (ATP-PFK), PPi-dependent phosphofructokinases (PPi-PFK or PPi-PFP), vacuolar proton pyrophosphatase (V-PPase), and vacuolar proton ATPase (V-ATPase). Our results show that genes encoding two ATP-PFKs (OsPFK4 and OsPFK5), one PPi-PFK (OsPFPA3), and one V-PPase (OVP3) were significantly induced in hypoxia samples of our five tolerant genotypes, while the degree of induction in the IR64 genotype was merely detectable (Table S10). The expression levels of PPi-PFK genes in our hypoxia samples were significantly higher when compared to ATP-PFK, which suggests that elongation of the coleoptile may rely heavily on PPi-PFK to phosphorylate fructose-6-phosphate to fructose-1,6-bisphosphate in glycolysis. In addition, the vacuolar proton pyrophosphatase gene family has been shown to actively pump H⁺ from the cytosol into the vacuole to maintain cytosolic pH homeostasis (Maeshima, 2000). In particular, the OVP3 gene (LOC_Os02g55890) has been previously reported as responsive to anoxia (Liu et al., 2009); this gene was detected as differentially regulated across the six genotypes in this study, suggesting an important role in submergence tolerance.

Cell elongation contributes significantly to coleoptile growth, and the genes encoding several cell wall loosening proteins, such as expansins, have been shown to be uniquely expressed under anoxia (Lasanthi-Kudahettige et al., 2007). Thus, to further investigate the contribution of cell wall-related genes to coleoptile elongation across our diverse genotypes, we examined the responses of over 100 genes involved in cell wall growth and loosening from the literature (Cosgrove, 2005) and from MapMan annotation bins. The results of this analysis show that 20 genes were differentially regulated (P < 0.05) by submergence in at least one genotype (Figure 6, Table S11). Interestingly, however, these results did not demonstrate a consistent expression pattern across diverse genotypes; for example, the gene for expansin B6 (EXPB6, LOC_Os10g40700) was significantly induced in the 8391 and 8753 genotypes but not in IR64 and the two RILs varieties, while the pectinesterase gene (LOC_Os04g51340) was only up-regulated in the two RILs. The results show that in the 8,391 genotype which exhibited the strongest coleoptile elongation significant induction of five cell wall-related genes was observed. At the same time, xyloglucan endotransglucosylase/hydrolase (XET) and pectinesterase genes were only slightly induced in this genotype (Figure 6).

In addition to genes related to the cell wall, those that respond to ethylene also play important roles in tolerance to submergence, including the ethylene response factor (ERF) genes Sub1 and SNORKEL (Fukao et al., 2006; Xu et al., 2006; Fukao and Bailey-Serres, 2008a,b; Hattori et al., 2009). Because previous studies have demonstrated increased ethylene accumulation in young seedlings of the flood-tolerant rice variety “Khao Hlan On” in comparison to the sensitive IR42, this suggests that ethylene biosynthesis genes are involved in submergence tolerance during germination (Ismail et al., 2009). Indeed, of the candidate genes identified in this study, two ERF genes (LOC_Os01g21120 and LOC_Os07g47790) were induced in all genotypes, although higher expression was seen in tolerant lines compared to the sensitive IR64 (Figure S6). Hence, because ethylene biosynthesis relies on the rapid oxidation of l-aminocyclopropane l-carboxylic acid (ACC), and because 22 genes that encode ACC oxidase (ACO) have been annotated in the rice genome, we investigated their expression patterns (Table S12). The results of this investigation show that one ACO gene (LOC_Os01g39860) was strongly induced in all genotypes when submerged, while another (LOC_Os09g27820) was expressed at higher levels in the five tolerant genotypes but was barely detectable in IR64.

In summary, although cell elongation and ethylene signaling pathways play important roles in the submergence response of rice seedlings, with some genes conservatively regulated in all genotypes, still more are regulated in a genotype-specific pattern, contributing to the higher tolerance of these lines. The results of this study illustrate the complexity of transcriptomic fine-tuning in response to submergence in rice seedlings from diverse genotypes.

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.