Spatio-temporal resolution of mRNA profiles during rice grain development
Gene expression was analyzed in five sub-regions of the developing rice grain comprising the CC, NE, OVT, EN and AL tissues that were collected by laser-capture micro dissection at 4, 8 and 16 DAF. RNA sequence reads were quality checked and over 20,000 genes were identified with statistically significant expression (reads higher than 10, p < 0.001). Genes with expression specific to each sub-region and grain development stages were compared using hierarchical clustering and principal component analysis (PCA). The biological replicates from individual sub-regions clustered together and PCA collectively identified 79% variance in our dataset (Fig. 1, Fig. S1). The dataset was further analyzed to separate genes for their tissue- and stage-specific expression. The genes with a preferentially higher expression level (with the highest number of reads, p < 0.001) in a particular sub-region or at a particular developmental time point were assigned to be tissue-specific and stage-specific differentially expressed genes (DEGs), respectively. In total, 10,037 DEGs were identified. In order to further characterize the sub-region and stage-specific genes, gene ontology (GO) analysis was conducted on DEGs using PlantGSEA (Yi et al., 2013). Significantly overrepresented GO terms and metabolic pathways (p < 0.01) were used for further comparison. Selected genes were also analyzed for their expression in these tissues by quantitative PCR (qPCR), generally validating the RNA-Seq results (Fig. S2).
Tissue-specific DEGs. In OVT, 635 and 753 DEGs were identified at 4 and 8 DAF, respectively. The most significantly enriched GO categories for these genes are iron binding, transporter, transcriptional activity, and phytohormone biosynthesis and signalling. Among the 1,123 DEGs identified as OVT-specific at 16 DAF, the over-represented GO categories included protein translation, ribosomal activity, protein folding and macromolecular metabolic process (Fig. 2a, Fig. S3a and Table S1). In the EN-specific set of genes, 143, 186, and 185 DEGs were identified at 4, 8 and 16 DAF, respectively. The enriched GO terms for EN-specific genes were associated with metabolism and providing energy, such as starch and glucose metabolism, glycolysis, transportation, and protein and lipid metabolism. Starch biosynthesis and microtubule-associated genes were also highly expressed in EN at 4 and 8 DAF, suggesting that endosperm cells are still enlarging during grain filling (Fig. 2b and Table S1). In the case of NE, one DEG at 4 DAF, 147 DEGs at 8 DAF, and 16 DEGs at 16 DAF were identified. The enriched GO terms for the 8 DAF NE-specific DEGs included transferase, metabolic activity, transcription, and coenzyme binding, consistent with the role of the NE in supporting EN development (Fig. 2c). Similarly, consistent with the suggested role of CC during grain development, the overrepresented GO terms among the 53, 572, and 99 DEGs identified in the CC tissue at 4, 8 and 16 DAF, respectively, included photosynthesis, transporter activity, metal binding, catalytic activity, and metabolic processes (Fig. 2d and Table S1). Since AL tissue was collected only at 16 DAF, we performed pairwise comparisons of gene expression in the AL tissue with CC, EN, NE, and OVT tissues at 16 DAF. This revealed 1,438, 3,829, 178, and 2,550 DEGs, respectively. Genes associated with macromolecule metabolism, amino acid metabolism, secondary metabolism, and lipid metabolism as well as vitamin B biosynthesis and transport were expressed in AL. These results suggest an important role of the AL in synthesis and transport of metabolites during grain filling in rice. Additionally, some stress responsive genes were also identified among the AL-specific DEGs (Fig. 2e, Fig. S3b and Table S2).
Stage-specific DEGs. Among the DEGs, 1,752 DEGs at 4 DAF, 773 DEGs at 8 DAF, and 2,190 DEGs at 16 DAF were identified as growth stage-specific genes (Fig. 1). The genes in the 4 DAF-specific set were related to chromatin assembly, DNA replication, cell mitosis and division in NE, EN and CC tissues, and to hormone regulation and glycolysis in OVT. The 8 DAF-specific genes were related to major CHO metabolism and transporter activity in OVT and EN. Over-represented GO categories of stage-specific gene sets indicate a transition to translation and protein biosynthesis to establish protein localization and secondary metabolism, especially in CC, OVT and EN tissues at 16 DAF (Fig. 3), which is consistent with previous studies on seed maturation (Gillies et al. 2012).
Diverse and complex regulatory networks controlling grain filling in rice
In order to explore the role of transcription factors (TFs) in rice grain development, we surveyed about 1,500 expressed genes encoding TFs in the rice genome. Of these, 470 TF genes were differentially expressed in our data set (Fig. 4a, and Table S3). The OVT tissue had the maximum number of DEGs encoding TFs, comprising 21.3%, 21.5% and 15.5% of TFs at 4, 8 and 16 DAF, respectively (Fig. 4b). Diverse types of TF-encoding genes were found to be differentially expressed, with the majority of these OVT-specific TFs belonging to bHLH and MYB families of TFs (Fi. 4c). Several of these TF family members, such as ABI3/VP1, ARF, Aux/IAA, B3, ARR and GRAS, have important roles in hormone regulation. Additionally, the TFs known for their roles in vascular development, cell growth regulation and differentiation, such as MYB, SBP, MADS, NAC, bZIP and HB family members, were also differentially expressed in the OVT. Growth stage-specific differences in TF expression profiles were also observed by enrichment analysis and co-expression network analysis. Genes encoding members of MADS and bZIP TF families showed higher expression at 4 and 8 DAF, while those encoding MYB, HB and WRKY TFs had elevated expression at 16 DAF in the OVT (Fig. 4d, Fig. S4a). The CC tissue-specific dataset revealed a higher proportion of TFs as compared to NE, EN and AL tissue-specific sets, especially at 8 DAF (13%) (Fig. 4e). The enriched TF families included WRKY, HB and Zinc finger TFs, most of which are known to be involved in seed development (Fig. 4e) (Zhang et al. 2011; Joseph et al. 2014). Chi-square test and co-expression network analysis showed that zinc-finger and AP2-ERF TF families were overrepresented in CC at 4 and 8 DAF (Table S3, Fig. S4b). In addition, nine TFs were highly expressed in EN at 16 DAF, and most of these are involved in seed maturation or germination (Kim et al. 2008; Lasserre et al. 2008). Two genes encoding AP2-like ethylene-responsive TFs (Os04g0653600 and Os05g0437100) and three heat shock TFs (Os02g0232000, Os06g0553100 and Os10g0419300) were more highly expressed in the AL tissue (Table S3). OsWRKY71 (Os02g0181300) had the highest expression level in the AL tissue, as was previously reported (Zhang et al. 2004). Together, these data suggest that reprogramming events in the developing rice grains comprise a complex coordination of different TFs in each sub-region tissue and at different developmental time points.
Expression profiles of genes related to hormone biosynthesis and signalling
Since plant hormones regulate grain development, we analyzed the expression of genes related to auxin (IAA), gibberellin (GA), brassinosteroid (BR), cytokine (CK), abscisic acid (ABA) and ethylene (ET) biosynthesis and metabolism. In general, BR-, CK- and GA-related genes had higher expression levels at 4 and 8 DAF in all sub-regions. The genes related to IAA biosynthesis and metabolism had relatively stable expression in OVT during 4 and 8 DAF, but exhibited slightly decreased expression at 16 DAF (Fig. 5, Figure S5a and 5b, and Table S4).
The genes encoding IAA biosynthesis enzymes, including OsYUCCA7, OsYUCCA4, OsASA1 and OsASA2 (Mano and Nemoto 2012), were highly expressed in the OVT tissue at 4 and 8 DAF, and showed reduced expression at 16 DAF. Genes encoding IAA transporters and related TFs, such as Aux/IAA, were expressed higher in OVT at 4 and 8 DAF. IAA signalling genes, however, were highly expressed in the AL tissue as compared to the EN at 16 DAF (Fig. 5a). GA20ox is a key enzyme in GA biosynthesis (Fleet and Sun 2005) and among GA20ox encoding genes in rice, Os01g09300 had higher expression level in OVT at 4 DAF. Notably, GA-related genes were expressed at relatively low levels in EN at every time point (Fig. S5a). CK biosynthesis-related genes including OsRR2 and OsRR6 (Hirano et al. 2008) were also preferentially expressed in the OVT tissue at 8 DAF (Fig. 5b). Only few BR biosynthesis-related genes were differentially expressed in our datasets. Among them, OsDWARF (Os03g0602300) encodes the BR enzyme C-6 oxidase, which is involved in seed development and panicle elongation (Mori et al. 2002); this gene was differentially expressed in CC and OVT tissues at 8 and 16 DAF. Several of the BR signalling-related genes were found to be differentially expressed, especially in the OVT tissue. OsCYP genes involved in BR synthesis are known to regulate seed size (Wu et al. 2008); they were expressed highly in the OVT tissue at 4 and 8 DAF (Fig. 5c). Genes encoding enzymes for ABA biosynthesis exhibited higher expression in the AL and the EN tissues, while the ABA signalling-related genes were preferentially expressed in the AL tissue at 16 DAF (Fig. 5d). This result is in a good agreement with previous reports, in which genes involved in ABA biosynthesis and signalling, including NCED, AAO, bZIP and PP2C (HVA), were enriched in the EN in late developmental stages (Xue et al. 2012) (Fig. 5d). The expression pattern of genes encoding ethylene (ET) signal transduction genes was quite stable during 4, 8 and 16 DAF (Fig. S5b). The expression of ethylene responsive TF genes was low in EN at 16 DAF.
Transporters Are Preferentially Expressed In Ovt
In order to gain insights into nutrient transport routes in the rice grain, we extracted the list of genes encoding different transporters from MapMan (https://mapman.gabipd.org) and examined their expression patterns in our set of DEGs. Overall, the expression level of genes encoding nutrient transporters, including transport ATPases and sugar-, nitrate- and sulfate-transporters, was higher in OVT, CC, and NE at 4 and 8 DAF. The expression level of genes encoding amino acid and protein transporters was higher at 16 DAF in OVT, CC and AL. These genes in general had a low expression in EN as compared to the other tissues (Fig. 6, Supplementary Fig, S5c and 5d, and Table S5).
Fifteen genes encoding P- and V-type transport ATPases showed the highest expression at 8 DAF in OVT and CC. Three of these genes expressed higher in CC at 8 DAF and among these, Os03g0183900 is a plasma membrane H+-ATPase that is specifically expressed in seeds (Baxter et al. 2003) (Fig. 6a). Genes encoding protein and sugar transporters expectedly had higher expression levels at 4 and 8 DAF in EN, while their expression decreased at 16 DAF in EN but increased in AL (Fig. 6b and 6c). We also identified several genes for sugar transporters that were preferentially expressed in OVT at 4 and 8 DAF. Six putative sugar transporters showed higher expression at 4 DAF in the OVT (Os02g0229400, Os04g0453400, Os04g0453200, Os03g0197100, Os11g0135300, and Os03g0101300 (Fig. 6c), supporting the view that OVT has a significant role in transporting sugars and other metabolites to the developing EN and embryo. Three genes encoding putative sugar transporters (Os03g0594400, Os12g0132500, and Os04g0511400) had their highest expression levels at 4 DAF in EN. The expression of known genes encoding sucrose transporters, including OsSUT1, OsSUT2, OsSUT3, OsSUT4 and OsSUT5 (Zhu et al. 2002), was similar at 4, 8 and 16 DAF in all sub-regions. Most of the genes encoding amino acid transporters (Zhao et al. 2012), including OsAAP2, OsAAP3, OsAAP5, OsAAP6 and OsAAP7, were preferentially expressed at high levels in OVT but not in EN at 4, 8 and 16 DAF (Fig. S5c). However, OsCAT11, OsATL6 and OsPROT2 showed an increased expression in EN at 16 DAF (Fig. S5c). The genes encoding nitrate, phosphate and sulfate transporters were preferentially expressed in OVT at 4 and 8 DAF but preferentially expressed only later in NE, CC and AL at 16 DAF (Fig. S5d). For example, the high affinity nitrate transporter OsNAR2.1 (Feng et al. 2011) was expressed highly in OVT at 4 DAF (Fig. S5d). OsPHO1;2 and OsPHO1;1, which tightly regulate phosphate loading into xylem in rice and plant fitness (Secco et al. 2010; Jabnoune et al. 2013), were upregulated in OVT at 8 DAF and in CC at 16 DAF, respectively. OsSULT3;3 and OsSULT3;4, which showed highest expression levels in AL at 16 DAF and in OVT at 8 DAF, respectively, are known to be induced by heavy metal stress (Buchner et al. 2004; Kumar et al. 2011) (Fig. S5d). OsSULT3;5 was the only gene among the genes encoding sulfate transporters that showed the highest expression level in EN at 16 DAF (Fig. S5d).
The expression of genes encoding iron- and zinc-related transporters, such as NRAMP, ABC, HMA and ZIP (Lanquar et al. 2005; Kobayashi and Nishizawa 2012), was higher at 4 and 8 DAF, but significantly decreased at 16 DAF (Fig. 6d, Fig. 7 and Tables S5 and S6). As expected, EN had the lowest expression of genes encoding metal transporters, although they were more highly expressed in AL at 16 DAF (Fig. 6d). Expression of NRAMP gene family members including OsNRAMP1, OsNRAMP2, OsNRAMP3 and OsNRAMP5 was higher in OVT and CC but lower in EN at 4 and 8 DAF (Fig. 6d). Similarly, a group of HMA family genes showed significantly higher expression level in OVT at 4 and 8 DAF, and in NE and AL at 16 DAF (Fig. 7b). Among the ZIP family genes, OsZIP4 and OsZIP7 showed higher expression level in OVT at 4 and 8 DAF, while OsZIP3 showed highest expression in EN (Fig. 6d). Genes encoding proteins for iron acquisition and transport such as OsFRDL2 (Yokosho et al. 2016), OsOPT, OsYS and OsCOPT3 were also expressed more highly in OVT and CC at 4 and 8 DAF (Fig. 6b and 6d and Fig. 7c). Together, the expression pattern of genes encoding metal transporters is consistent with the low levels of iron and zinc in the rice endosperm that is consumed as polished rice.
Novel cis-regulating elements associated to OVT and CC specific genes
Since a significant number of genes that were differentially expressed in OVT and CC encode transcription factors as indicated by GO categorization, we analyzed the DEG sequences to identify known or novel TF regulatory DNA sequence motifs in their promoter or coding regions. We searched approximately 1000 bp proximal and distal to the ATG start codon for sequence motifs and identified them using the MEME and JASPAR plant databases. Several common motifs for TF families associated with seed development were identified among the OVT- and CC-specific DEGs (Table 1, upper panel) (Papi et al. 2000; Sreenivasulu et al. 2006; Wang et al. 2010). For example, the motifs for the two MADS family TFs SOC1 and PI were found in OVT and CC DEGs at 4, 8 and 16 DAF. Similarly, the high mobility group box HMG-I/Y for the HMG TF family and the idl motif that is bound by the zinc finger TF family were also enriched in our dataset during all three development stages. Several motifs were enriched only in the OVT specific DEGs. Of note, ABI3 (GCATG), a binding site of ABI/VP3 transcription factor related to auxin signaling, and BES1 (CACGTG) motif, related to BR signalling were enriched only in OVT-specific genes at 8DAF, suggesting an important role of phytohormone related transcription factors in early grain filling. At 16 DAF, AtMYB84 (GGTnGGT) and AtSPL8 (GTAC) were enriched motifs in OVT but no motifs specific to CC were detected as enriched at this development stage (Table 1, upper panel). We also identified a novel motif “GGAGGA” in OVT and CC specific genes at 8 DAF. From the analysis using GO-MEME (Buske et al. 2010) and cross comparison with other plant species, it is likely that this motif is involved in transporter and hormone regulation. Another motif “GCCGCC” was identified as unique, but a high GC content motif is suggested to be commonly found in rice gene promoters (Lenka et al. 2009, 2011). In the distal regions of OVT and CC specific genes, we found motif “GCATGC” as enriched in OVT specific genes at 4 and 8 DAF, and motif ”GAGAGA” enriched in CC specific genes at 8 DAF (Table 1, lower panel). We compared these two motifs with known motifs using TOMTOM (Gupta et al. 2007) software. “GCATGC” was similar to BR signalling-associated regulatory element BES1 and “CCTCC” was similar to MYB regulatory elements (Table 1, lower panel).
Summary of overrepresented cis-regulatory DNA sequence motifs identified from MEME analysis