Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Genome-wide analysis of spatiotemporal expression patterns during rice leaf development

BMC Genomics volume 22, Article number: 169 (2021) Cite this article

Abstract

Background

Rice leaves consist of three distinct regions along a proximal-distal axis, namely the leaf blade, sheath, and blade-sheath boundary region. Each region has a unique morphology and function, but the genetic programs underlying the development of each region are poorly understood. To fully elucidate rice leaf development and discover genes with unique functions in rice and grasses, it is crucial to explore genome-wide transcriptional profiles during the development of the three regions.

Results

In this study, we performed microarray analysis to profile the spatial and temporal patterns of gene expression in the rice leaf using dissected parts of leaves sampled in broad developmental stages. The dynamics in each region revealed that the transcriptomes changed dramatically throughout the progress of tissue differentiation, and those of the leaf blade and sheath differed greatly at the mature stage. Cluster analysis of expression patterns among leaf parts revealed groups of genes that may be involved in specific biological processes related to rice leaf development. Moreover, we found novel genes potentially involved in rice leaf development using a combination of transcriptome data and in situ hybridization, and analyzed their spatial expression patterns at high resolution. We successfully identified multiple genes that exhibit localized expression in tissues characteristic of rice or grass leaves.

Conclusions

Although the genetic mechanisms of leaf development have been elucidated in several eudicots, direct application of that information to rice and grasses is not appropriate due to the morphological and developmental differences between them. Our analysis provides not only insights into the development of rice leaves but also expression profiles that serve as a valuable resource for gene discovery. The genes and gene clusters identified in this study may facilitate future research on the unique developmental mechanisms of rice leaves.

Background

Leaves, which are the main site of photosynthesis in higher plants, are usually polarized along three axes: proximal-distal, adaxial-abaxial, and medial-lateral. Tissues arranged along these axes have characteristic morphologies and functions. As leaves are derived from immature cell populations protruding from the shoot apical meristems (SAM), their morphology and functions must be acquired during the course of development. Leaf development is a tightly orchestrated process incorporating multiple events crucial to organogenesis: axis determination, pattern formation, and identity establishment. Additionally, the growth of leaf primordia, which relies on cell proliferation and differentiation, is precisely regulated both temporally and spatially to produce typically shaped leaves.

The morphology of leaves varies greatly among species and developmental phases and environments, and this variation is driven by differences in leaf genetic programs among species. Hence, the mechanisms regulating leaf morphogenesis should be studied in a wide variety of species. Most information currently available has been obtained from analyses of the model eudicot plant Arabidopsis. A number of genes regulating leaf development have been identified in Arabidopsis [1], and the molecular mechanisms of leaf morphogenesis in various species have been elucidated based on information obtained from Arabidopsis.

Grasses belong to the monocot clade, and their leaf morphology is distinct from that of Arabidopsis. Although grass leaves are polarized along the same three axes as those of other plants, they are unique in that distinct regions with differing morphology and function are located along the proximal-distal axis (Fig. 1a). The leaf blade is the distal part of the leaf; it has a flat structure and is rich in mesophyll cells, in which photosynthesis occurs. The leaf sheath is located at the basal part of the leaf and has a thick structure that protects inner leaves and provides structural support to the blade. The boundary region between the blade and sheath comprises the lamina joint, ligule, and auricle. The lamina joint acts as a hinge that allows the blade to bend abaxially, thereby optimizing light capture by the blade. Each of the three regions undergoes different developmental processes. In addition, spatiotemporal coordination of tissue differentiation during development contributes to the final leaf morphology. As with Arabidopsis, tissue differentiation in grass leaves proceeds in the basipetal direction, suggesting that these processes are under precise spatial and temporal control by genetic mechanisms.

Fig. 1
figure 1

Overview of the samples used for microarray analysis. a Morphology of mature rice leaves. A, auricle; B, leaf blade; J, lamina joint; L, ligule; S, leaf sheath. b Samples used for microarray analysis. A seedling and the separated leaves of the seedling are shown. Yellow boxes indicate the 12 parts sampled for microarray analysis. Arrows show the boundaries between the leaf blade and sheath. See Table 1 for details. c Magnified view of the leaf primordium at stage P3. d Magnified view of the shoot apex containing the shoot apical meristem and P1 and P2 leaf primordia. Scale bars: 1 cm in (a), (b); 20 μm in (c), (d)

Full size image

To date, several genes that are important for grass leaf morphology have been identified. Related to the organs and tissues that are differentiated along the proximal-distal axis, LIGULELESS1, a member of the SQUAMOSA PROMOTER BINDING-LIKE (SPL) gene family, is essential for the differentiation of organs in the blade-sheath boundary region [2,3,4]. Homologs of Arabidopsis BLADE-ON-PETIOLE genes in rice are important for sheath development [5]. Class I KNOX genes important for the maintenance of the SAM are believed to provide proximal cues for leaf primordia [6, 7]. Meanwhile, cell proliferation patterns have been reported to affect grass leaf morphology. Cell proliferation in immature tissues of leaf primordia is controlled by protein complexes encoded by two gene families, GROWTH-REGULATING FACTORS (GRFs) and GRF-INTERACTING FACTORS (GIFs) [8,9,10,11]. Changes in their protein complex composition reportedly serve as a switch for the transition from cell division to cell expansion [12]. However, few of the genes that play important roles in grass leaf development have been identified.

For elaboration of leaf development and morphology, gene expression should be precisely regulated both temporally and spatially. Thus, revealing the expression patterns of genes during development would contribute significantly to understanding the genetic mechanisms behind the process of leaf establishment. Transcriptome analysis is a powerful method for exploring gene expression dynamics at both the genome-wide and single-gene levels. To date, transcriptome analysis of leaf development has been performed in species including Arabidopsis [13, 14], maize [15,16,17,18,19,20,21,22], and rice [23,24,25]. In particular, transcriptome changes accompanied by tissue differentiation have been intensively studied in the developing leaf blade in maize. However, no study to date has reported the temporal transcriptomic changes occurring from leaf initiation to leaf maturation in rice. Furthermore, most transcriptome studies of grass leaves have been performed on tissues from only one of the regions along the longitudinal axis, making it difficult to draw direct comparisons among regions. Therefore, to fully elucidate gene expression profiles and characterize gene function, it is necessary to investigate spatiotemporal changes in leaf transcriptomes from each region along the longitudinal axis of the leaf.

In this study, we performed transcriptome analysis of rice leaf development using the Agilent rice 44 K microarray, which is compatible with the rice expression database RiceXPro [26, 27]. Our experimental design included a broad range of developmental stages and several distinct regions along the leaf longitudinal axis, which allowed us to capture overall transcriptome dynamics throughout leaf development. Our data analysis uncovered trends in the expression patterns of certain gene clusters during leaf development and revealed relationships between developmental events and those gene clusters. In addition, we performed in situ hybridization with 49 selected genes based on the data from our transcriptome analysis. As a result, we identified multiple genes with localized expression in tissues characteristic of grass leaves. The present work provides a foundation for future analyses of genes with novel functions in rice leaf development.

Results

Experimental design for microarray dataset

Rice leaf ontogeny, i.e., the developmental process from initiation to maturation, is described in Itoh et al. (2005) [28]. Briefly, according to the staging system based on plastochron numbers (Pn), the P1 leaf primordium protrudes from the SAM and then grows to surround the SAM at stage P2. During the P1 and P2 stages, the leaf primordium consists of undifferentiated cells with no morphological characters. During the P3 stage, the boundary between the blade and sheath is established, and the future blade and sheath parts can be distinguished. In addition, the ligule primordium is formed in the boundary region at this stage. Although most of the P3 leaf primordium is comprised of undifferentiated cells, the outermost cells on the distal side of the primordium begin to differentiate into epidermal cells. During stage P4, the leaf blade elongates rapidly, and the difference between the blade and sheath becomes more pronounced. The P4 leaf primordium exhibits a clear gradient of cell states along its longitudinal axis; cells in the proximal region remain undifferentiated, whereas those in the distal region are differentiated. During stage P5, the leaf sheath elongates rapidly, and the growth and maturation of the leaf are completed by the P6 stage, whereas bending of the lamina joint occurs between stages P5 and P6.

To obtain a comprehensive transcriptome of leaf development in rice, we sampled 12 leaf parts representing various stages and components along the longitudinal axis (Fig. 1b–d; Table 1). Rice seedlings at the four-leaf stage were dissected into 12 parts: shoot apex containing the SAM and P1 and P2 leaf primordia (Fig. 1d); entire P3 leaf primordium (Fig. 1c); apical, middle, and proximal parts of P4 leaf blade; P4 leaf sheath; and the leaf blade, sheath, and boundary region of P5 and P6 leaves. Three biological replicates were prepared for each part, and their RNA was hybridized to a 44 K rice microarray (Agilent Technologies, Santa Clara, CA) [26, 29, 30]. Out of 43,144 probes corresponding to 29,864 genes on the rice 44 K microarray platform, 31,996 probes corresponding to 24,022 genes were expressed in at least one sample. Normalized expression levels of these 31,996 probes corresponding to 24,022 genes were used in this study. Distribution of the normalized expression levels of those probes for each sample roughly exhibited the normal distribution centered at zero (Supplemental Figure 1). Pearson correlation analysis showed strong correlations among the three replicates, indicating that our dataset was highly reproducible (Supplemental Figure 2). Furthermore, we verified expression profiles of six genes in P4 stage by real-time RT-PCR analysis (Supplemental Figure 3). The expression patterns of all the genes were consistent with the patterns from microarray analysis. This suggests that our expression data of microarray is considerably accurate and reliable.

Table 1 Description of samples used for microarray analysis
Full size table

Transcriptome dynamics during rice leaf development

To elucidate the transcriptome dynamics that occur during leaf development, principal component analysis (PCA) was performed using all samples. The first, second, and third principal components (PC1, PC2, and PC3) explained 60.9, 13.2, and 8.1% of the total variance among samples, respectively (Supplemental Figure 4; Additional file 3). Plotting the samples within the three-dimensional space defined by PC1, PC2, and PC3 allowed the relationships among the samples to be visualized, reflecting the two properties of tissue differentiation state and tissue identity (Supplemental Figure 4a; Additional file 3). However, these properties were poorly represented by each principal component (Supplemental Figure 4b and c) due to the fact that PC2 was excessively attracted toward P4Bm, which was located distant from the other samples. Generally, PCA is not robust to outliers, and principal components tend to be attracted toward them, interfering with the detection of the overall dataset structure. To avoid excessive attraction of PC2 toward the outlier P4Bm, we modified PC1, PC2, and PC3 while retaining positional relationships between the samples, as follows. PC1, PC2, and PC3 scores were treated as variables, and PCA was again applied using all samples except P4Bm, resulting in modified principal components (mPC1, mPC2, and mPC3) that were not excessively affected by P4Bm. Altogether, mPC1, mPC2, and mPC3 explained 60.9, 11.9, and 9.5% of the total variance among samples, respectively, and successfully captured the characteristic patterns of the dataset (Fig. 2; Additional file 4).

Fig. 2
figure 2

Principal component analysis (PCA) score plot of samples based on modified principal components. a The space defined by mPC1, mPC2, and mPC3. A three-dimensional model allowing interactive rotation is available in additional file 4. See Supplemental Figure 4 for the PCA score plot based on the original principal components. b The space defined by mPC1 and mPC2. The gray arrow represents the regression curve for the distribution of samples in this space. c The space defined by mPC1 and mPC3. The proportions of the total variance explained by mPC1, mPC2, and mPC3 are shown in parentheses. Samples collected at the same stage are shown in the same color. Samples with different tissue identities are indicated by different symbols: shoot apex, square; P3 leaf, circle; blade, triangle; blade-sheath boundary, diamond; sheath, inverted triangle

Full size image

Two groups were separated by mPC1, namely, immature tissues represented by SA, P3, P4S, and P4Bb, mature tissues represented by P4Ba, and samples derived from P5 and P6 stage leaves. This result suggests that mPC1 represented the differences between immature and mature tissues (Fig. 2b). Conversely, mPC2 characterized samples with intermediate tissue differentiation, most notably P4Bm, suggesting that mPC2 represented the transient state of the transcriptome during tissue differentiation (Fig. 2b). Thus, an arrow fitting the distribution of the samples in the space defined by mPC1 and mPC2 would represent the change in transcriptome dynamics associated with tissue differentiation from the immature state through the transient state to the mature state. Collectively, mPC1 and mPC2 explained 72.8% of the total variance among samples, suggesting that tissue differentiation state has profound effects on the leaf transcriptome. Moreover, samples derived from the P4 leaf exhibited large transcriptomic variations, whereas all P4-stage leaf samples including sheath samples were aligned along the arrow. This result indicates that the shift in the transcriptome associated with leaf maturation is found throughout the leaf during P4, coinciding with intensive basipetal tissue differentiation at stage P4 [31].

In addition, mPC3 separated samples of the leaf blade—P4Bm, P4Ba, P5B, and P6B—from those of the leaf sheath and blade-sheath boundary region—P5S, P6S, P5BS, and P6BS, indicating that mPC3 represented differences between the leaf blade and sheath (Fig. 2c). On the other hand, only slight differences were observed among immature leaf samples such as P3, P4S, and P4Bb, suggesting that the transcriptomic difference between the leaf blade and sheath becomes more pronounced during maturation.

Overall, our results suggest that the transcriptome of each part of the leaf changes with the progression of tissue differentiation and the acquisition of tissue identity.

Gene expression patterns during leaf development and their associations with gene function and transcriptional regulation

To uncover the major gene expression patterns during rice leaf development, we conducted cluster analysis of genes based on their expression patterns. Prior to cluster analysis, analysis of variance (ANOVA) was applied to detect differentially expressed genes among different parts of the leaf. Of 31,996 probes corresponding to 24,022 genes, 31,043 probes corresponding to 23,350 genes were extracted (p-value = 0.001 when adjusted for the false discovery rate [FDR]). K-means clustering was performed on these probes, and 28 clusters with distinct expression patterns were obtained (Supplemental Figure 5; Supplemental Table 1). For some of the clusters, there were large differences in expression levels among samples and characteristic expression patterns, suggesting that some groups of genes undergo similar changes in gene expression, and that such changes are associated with events during leaf development.

To evaluate how the gene expression patterns and dynamics of these gene clusters are related to the functions of the genes, we conducted Gene Ontology (GO) enrichment analysis on each cluster (Supplemental Figure 6). In addition, given the importance of transcriptional regulation to development, transcription factors and transcriptional regulators were extracted from each cluster (Supplemental Figure 7). These analyses identified the characteristic functions of genes within each cluster, including several genes that may be involved in specific processes during rice leaf development (Fig. 3; Table 2).

Fig. 3
figure 3

Six clusters showing expression patterns supposedly associated with events during leaf development. Grey lines indicate the expression profiles of probes in each cluster. Red lines indicate the mean of all probes in each cluster. Three biological replicates are summarized by median. The number of genes in each cluster is indicated in the upper right of each panel. See Supplemental Figure 5 for all 28 clusters obtained through K-means clustering

Full size image
Table 2 Enriched GO terms, transcription factors, and transcriptional regulators in the six clusters in Fig. 3. See Supplemental Figures 6 and 7 for enriched GO terms and TF/TRs in all 28 clusters obtained through K-means clustering, respectively
Full size table

Cluster 1, a group of genes that was specifically expressed in the shoot apex, was enriched in genes involved in transcriptional regulation (regulation of transcription, p = 4.8e-07). Within this cluster were class I KNOX genes, which are important for SAM maintenance [32], and OsNAM/OsCUC3 genes, which may be involved in organ boundary formation [33]. Thus, Cluster 1 was predicted to include genes related to SAM function and leaf initiation. This cluster also included genes in the BBM clade of the PLETHORA family [34], which are expressed in crown root primordia [35], and OsTB1, which is expressed in axillary buds [36]. Because the shoot apex tissue used in this study contained stem tissue as well as the SAM and leaf primordia, the presence of root- and axillary bud-related genes in this cluster was not surprising.

Cluster 2 contains genes that were highly expressed in tissues undergoing active cell proliferation. In this cluster, GO terms associated with cell division and cytokinesis (microtubule-based movement, p = 8.6e-05) were detected. Moreover, it contained ANT clade genes of the PLETHORA family [34], GRF family genes [12], and OsGIF1/MKB3 [11], which have been described as promoters of cell proliferation in leaf primordia. Thus, Cluster 2 was expected to contain important genes related to cell proliferation in leaf primordia.

Cluster 9 genes were highly expressed in the middle parts of the P4 leaf blade and P5 leaf sheath. GO analysis revealed that this cluster contains class III peroxidases (response to oxidative stress, p = 4.2e-07). Some class III peroxidases regulate reactive oxygen species homeostasis in the apoplast, thereby affecting cell-wall stiffness [37, 38]. GO terms for cell-wall remodeling enzymes including XTHs (carbohydrate metabolism, p = 1.3e-3) [39] were also enriched in Cluster 9. Thus, Cluster 9 appears to be enriched in genes involved in the control of cell-wall extensibility and cell elongation.

Cluster 12 genes were mostly expressed at high levels in the mature leaf blade. GO terms for genes involved in photosynthesis (photosynthesis, p = 7.9e-21; carbon utilization by fixation of carbon dioxide, p = 8.5e-10; and electron transport, p = 3.3e-07) were enriched. This cluster included C2C2-CO-like family genes [40] and phytochrome-interacting bHLH factors (PIFs) [41]. Thus, Cluster 12 was expected to contain a high concentration of genes associated with photosynthesis and light-mediated signal transduction.

Cluster 7 genes were highly expressed in immature samples and samples from mature tissues from the sheath and blade-sheath boundary region. This cluster was enriched in genes involved in carbohydrate metabolism (main pathways of carbohydrate metabolism, p = 3.7e-05). The leaf sheath is believed to act as sink tissue for carbohydrates prior to heading [42]. In addition, this cluster includes OsBOP genes that are important for sheath development [5]. Thus, Cluster 7 was enriched in genes involved in the carbohydrate sink function and sheath development.

Cluster 15 consisted of genes that were preferentially expressed in the mature sheath and blade-sheath boundary region. This cluster was enriched in genes related to GO terms for protein kinases (protein amino acid phosphorylation, p = 6.2e-05). The blade-sheath boundary contains the lamina joint, which bends between stages P5 and P6. Various phytohormones and environmental stressors affect the bending process [43, 44], and many protein kinases exhibit temporal changes in expression during lamina-joint bending [25]. Thus, we assumed that Cluster 15 contains genes involved in lamina-joint bending.

Taken together, these analyses revealed genome-wide gene expression patterns during rice leaf development, and relationships between expression pattern and function were found in certain gene clusters.

Identification of genes with localized expression during leaf development

To identify novel genes that play important roles in early development and subsequent morphogenesis and tissue formation in the rice leaf, we selected genes from Clusters 1 to 8 using K-means clustering analysis. These clusters contain genes that tended to be highly expressed during the early stages and weakly expressed during the later stages, and thus were expected to include candidate genes. Forty-nine genes, most of them involved in transcriptional regulation, were selected from the gene clusters, and their spatial expression patterns around the shoot apex were examined through in situ hybridization. Examples of these gene expression patterns are described below.

Cluster 2 included the PLATZ family transcription factor Os02g0172800, which is a co-ortholog of ORESARA15 [45]. This gene was expressed in the basal part of immature leaves, and its expression was strongest in the abaxial side of the leaf primordium (Fig. 4a, b).

Fig. 4
figure 4

Spatial expression patterns of several genes identified through K-means clustering analysis. The expression patterns of (a, b) Os02g0172800, c OsbHLH080, d, e OsbHLH166, f Os05g0363500, and g, h OsGH3–4 around the shoot apex. a, c, d, f, and g are longitudinal sections of the shoot apices, while (b), e, and h are transverse sections of the shoot apices. The cluster number to which each gene belongs is indicated in the upper left. L, ligule; A, aerenchyma; D, diaphragm; V, vascular. Scale bars: 100 μm

Full size image

Tissue-specific expression patterns of three genes in Cluster 3 were detected. Expression of OsbHLH080, a member of the bHLH gene family, was observed in the abaxial base of leaf primordia and developing ligules (Fig. 4c). Another bHLH gene, OsbHLH166, was expressed mainly in the presumptive lysigenous aerenchyma areas of leaf primordia (Fig. 4d, e). Os05g0363500, which encodes a WD40 repeat-containing protein, was expressed in the developing diaphragms of leaf sheaths (Fig. 4f).

OsGH3–4, a member of the GH3 family involved in auxin conjugation, was placed in Cluster 7 and exhibited elevated expression levels mainly in the leaf sheath and blade-sheath boundary region. Expression of OsGH3–4 was detected in both the central domain of the SAM and tissues adaxially adjacent to vascular bundles in the leaf sheath and midrib, where bundle sheath extension cells differentiate (Fig. 4g, h).

Additionally, Cluster 7 contained four OsARF paralogs belonging to the ARF6/8 subfamily [46]. These genes, OsARF6/12/17/25, were expressed at the basal part of the leaf primordium around stage P3, with especially high expression levels at the margin (Fig. 5a–d, arrowheads). Moreover, these four OsARFs were also expressed in the developing ligule and marginal parts of the blade-sheath boundary region that were expected to differentiate into auricles (Fig. 5a–d). In addition to these common expression patterns among the four OsARFs, the characteristic expression of each gene was also observed. OsARF6 was expressed in the epidermis of the basal region of P4 leaf blades (Fig. 5a), whereas OsARF12 was expressed throughout the sheath and basal region of the blades (Fig. 5b). OsARF17 expression was detected in the adaxial epidermis at the blade-sheath boundary at stage P3 and in the lamina joint of P4 leaves (Fig. 5c), whereas OsARF25 was strongly expressed in the lamina joint of P4 leaves and at the base of sheaths in stages P4 and P5 (Fig. 5d).

Fig. 5
figure 5

Spatial expression patterns of four OsARFs and the accumulation pattern of mature miR167. The expression patterns of a OsARF6, b OsARF12, c OsARF17, d OsARF25, and e mature miR167 around the shoot apex in longitudinal sections. Arrowheads in (a)–(d) indicate expression in the basal part of the leaf margin. L, ligule; M, marginal part of the blade-sheath boundary. Scale bars: 100 μm

Full size image

In addition, ARF6/8 orthologs including the four OsARFs listed above are known targets of miR167, which is an evolutionarily conserved microRNA in seed plants [46]. To clarify the relationships between the expression patterns of the four OsARFs and miR167, the accumulation patterns of mature miR167 in the shoot apex were examined through in situ hybridization using a probe containing BNANC, which is a bridged nucleic acid derivative of a locked nucleic acid (Fig. 5e). Signals specific to mature miR167 were detected at higher levels in the distal part and lower levels in the basal part, indicating that the accumulation of miR167 and the expression of the four identified OsARFs were largely exclusive.

In addition to the genes described above, we explored a number of genes with localized expression during leaf development based on the list obtained through K-means clustering (Supplemental Figure 8; Supplemental Table 2). Therefore, our strategy of using K-means clustering facilitated the selection of genes that were differentially expressed among stages or tissues. In addition, the genes identified here should be further analyzed to uncover their functions in rice leaf development.

Discussion

Leaf primordia, which originate from the flank of the SAM, undergo various developmental processes from differentiation to maturation. Leaf primordia of grasses develop distinct regions along the longitudinal axis, and each region differentiates into a morphologically unique structure. Tissue differentiation in grass leaves during development proceeds in a basipetal direction, which has been extensively studied using the maize developing leaf blade as a model [16, 47,48,49,50,51,52,53,54]. These studies revealed that basipetal tissue differentiation is accompanied by dynamic changes in mRNA, proteins, and various metabolites. However, the resolution of these studies is limited spatially and temporally, as only parts of a single developing leaf blade at a specific stage were used, whereas other stages and tissues were not examined. Thus, our study is unique in terms of providing genome-wide expression profiles both across developmental stages and distinct tissues along the longitudinal axis.

In this study, the PCA results suggest that the leaf transcriptome dynamically changes among stages of tissue differentiation, as shown in the space defined by mPC1 and mPC2 (Fig. 2b). On the other hand, the difference between the leaf blade and sheath, which was represented by mPC3, accounted for less transcriptomic variation among samples than did mPC1 or mPC2 (Fig. 2c), suggesting that organ identity has a smaller effect on the leaf transcriptome than tissue differentiation. In particular, marked transcriptomic changes were observed in the P4 leaf along its longitudinal axis (Fig. 2b). The transcriptomes of the proximal parts, P4S and P4Bb, resembled that of the leaf primordium in early stages, such as P3. On the other hand, the distal part (P4Ba) exhibited transcriptomic similarities to leaf parts at later stages, such as P5 and P6. Thus, our results suggest that P4 is the stage at which tissue differentiation proceeds in the basipetal direction along the whole leaf and dynamic transcriptome changes occur.

Along the longitudinal axis of developing grass leaves, it has been suggested that cell proliferation and elongation occur actively in the basal and middle regions [55] and that photosynthetic activity and related gene expression is high in the apical region [16, 53]. Our cluster and GO analyses revealed that the genes in Clusters 2, 9, and 12 may be involved in cell proliferation, cell elongation, and photosynthesis, respectively (Fig. 3; Table 2). The expression patterns of these clusters in various parts of the P4 leaf are consistent with the previous findings; namely, different developmental events occur in a single leaf primordium along its longitudinal axis [16, 53, 55]. Thus, P4 is the stage wherein drastic developmental reprograming occurs in rice. Moreover, acquisition of photosynthetic competence is initiated between the P3 and P4 stages [24], which is in accordance with our finding that genes in Cluster 12 were weakly expressed at stage P3 and strongly expressed at stage P4 (Fig. 3; Table 2). In addition to these clusters, we identified Cluster 7 genes as being highly expressed in the sheath and blade-sheath boundary regions (Fig. 3; Table 2). OsBOP genes, which are regulators of leaf sheath identity, were included in this cluster and highly expressed at stage P3, suggesting that a number of upstream or downstream genes exhibit polarized expression patterns along the longitudinal axis from an early stage of leaf development.

The rice leaf has distinct morphological features that are not present in Arabidopsis. The development of these features is likely controlled in part by genetic mechanisms unique to grasses. Previous studies have identified several genetic factors in the morphogenesis of characteristic organs in grass leaves. LIGULELESS1 (LG1) is required for the development of structures at the boundary between the blade and sheath [2,3,4]. DROOPING LEAF is important for the development of the midrib [56, 57]. However, most of the genetic factors and their networks that underpin grass leaf morphogenesis remain unknown.

In general, to explore novel genes that play roles in the developmental process of interest, an efficient screening method is required. Expression profiling is one practical method of gene discovery. Transcriptome analysis and in situ hybridization are representative methods of expression analysis with different advantages. Transcriptome analysis can provide global expression profiles at the genome-wide scale, whereas in situ hybridization can reveal the spatial expression pattern of a gene at the tissue and cellular levels. By combining these methods, many studies have attempted to identify novel genes that function in a developmental process of interest [15, 22, 58,59,60,61,62,63,64,65].

In this study, we identified a number of genes with localized expression patterns that may be associated with morphological features unique to rice and grasses. For example, Os02g0172800 was expressed in immature leaf primordia at higher levels in the abaxial side of the basal part of the leaf primordia (Fig. 4a, b). This expression pattern is similar to that of MKB3, which is a positive regulator of cell proliferation in leaf primordia [11]. Furthermore, Os02g0172800 is a co-ortholog of ORESARA15, which reportedly promotes cell proliferation through a genetic pathway mediated by AN3, an ortholog of MKB3 [45]. It has been suggested that the formation of the leaf sheath requires a decreasing gradient of cell proliferation from the abaxial side of the leaf primordium to the adaxial side, which may be associated with the expression pattern of MKB3 [11]. Thus, Os02g0172800 might be also involved in cell proliferation in leaf primordia, and its expression pattern might reflect the unique cell proliferation pattern during the development of the leaf sheath.

Expression of OsbHLH80 was detected in the developing ligule (Fig. 4c). Ligule development was found to be disrupted by dysfunction of brassinosteroid (BR) signaling [66]. OsBC1, a paralog of OsbHLH80, is considered to regulate lamina-joint bending in response to BRs [67]. Thus, characterizing the roles of OsbHLH80 in ligule development in association with BR signaling is a worthwhile topic for future research.

OsGH3–4 was expressed in tissues adaxially adjacent to vascular bundles in the leaf sheath and midrib, where bundle sheath extension cell differentiation occurs (Fig. 4g, h). In addition, tissue-specific expression of OsbHLH166 and Os05g0363500 were observed in the presumptive region of lysigenous aerenchyma and developing diaphragms, respectively (Fig. 4d–f). During the development of the leaf sheath and midrib, a group of parenchyma cells collapse to form lysigenous aerenchyma, whereas parenchyma cells adaxially adjacent to the vascular bundles remain intact and differentiate into bundle sheath extension cells. Additionally, some parenchyma cells transform into stellate parenchyma cells in diaphragms that vertically separate the lysigenous aerenchyma [68]. Although these structures are well-developed in the leaves of rice and some wetland plants, the genetic networks regulating the development of these structures remain unknown. It has been demonstrated that the formation of constitutive aerenchyma in rice roots is regulated by auxin signaling [69]. Thus, OsGH3–4 may play an important role in the patterning of bundle sheath extension cells through the modulation of the spatial patterns of auxin accumulation in leaf primordia. Functional analysis of OsbHLH166 and Os05g0363500 would provide insights into the developmental programs underlying the formation of aerenchyma and elaboration of the unique cell shape of stellate parenchyma cells forming diaphragms.

OsARF6/12/17/25 exhibited polarized expression along the longitudinal axis and was expressed in developing ligules as well as the marginal parts of the blade-sheath boundary, where auricle differentiation occurs (Fig. 5a–d). Moreover, OsARF17 exhibited a unique expression pattern in the adaxial epidermis of the blade-sheath boundary in leaf primordia (Fig. 5c). Similarly, OsLG1 exhibited localized expression in marginal parts of the blade-sheath boundary (Supplemental Figure 9). Recent research suggested that a member of the ARF6/8 subfamily acts downstream of the OsLG1 ortholog in wheat [4]. Thus, these genes may be important for the development of organs in the blade-sheath boundary region. Furthermore, OsARF25 exhibited strong expression in the leaf-sheath pulvinus, which is a gravisensitive tissue in the leaf base that is involved in shoot bending in grasses [70] (Fig. 5d). This expression pattern suggests that OsARF25 may play an important role in the bending capability of shoots at the leaf-sheath pulvinus in response to gravity. In addition to the leaf-sheath pulvinus, OsARF25 was strongly expressed in the lamina joint (Fig. 5d). Several genes controlling shoot bending were previously reported to be expressed in both the lamina joint and leaf-sheath pulvinus [71,72,73]. These gene expression patterns suggest that the lamina joint and leaf-sheath pulvinus use similar genetic programs to achieve bending.

Post-transcriptional regulation by miR167 has been reported to be required for attaining the correct spatial expression patterns of ARF6/8 during reproductive development in Arabidopsis [74]. Our analyses revealed mutually exclusive expression of miR167 and OsARF6/12/17/25 along the longitudinal axis of leaf primordia. This result indicates that miR167 post-transcriptionally downregulates OsARF6/12/17/25 expression in the distal part of the leaf primordium (Fig. 5). A similar accumulation pattern for mature miR167 was reported in Arabidopsis cotyledons [75]. Thus, post-transcriptional regulation of ARF6/8 orthologs by miR167 in the distal part of the leaf primordium may be evolutionally conserved between rice and Arabidopsis.

Our microarray analysis had some limitations in resolution despite the sampling of extensive leaf stages and regions. For example, the sample SA contained the SAM, P1 and P2 leaf primordia, and immature stem tissues below the shoot apex, which makes it impossible to distinguish differences in expression profiles among these tissues (Fig. 1d). To improve our transcriptomic resolution, the separation of small tissues through laser microdissection, as performed in maize transcriptome studies on different domains of the SAM [15, 22], may facilitate more detailed gene expression profiling at the early stages of rice leaf development. In addition, it is possible that gene expression would change as plants grow mature. As a preliminary test for the possibility, we analyzed temporal expression profiles of six genes in shoot samples from the second to the forth leaf stage by real-time RT-PCR analysis (Supplemental Figure 3). Four genes are equally expressed among the three developmental stages (Supplemental Figure 3a to d), whereas OsPsbR3 showed a decreasing trend (Supplemental Figure 3e). OsLOX2;2 showed extremely low expression in all the developmental stages, presumably because the shoot tissue used in this experiment did not contain the distal part of leaf blade where this gene is highly expressed (Supplemental Figure 3f). The result indicated that expression of some genes could change even within the short growth period. To clarify expression profiles of genes throughout the developmental stages of plant life cycle, a large-scale transcriptome analysis is required in future studies. Moreover, the 44 K microarray platform used in this study does not cover all transcripts in the rice genome. RNA sequencing analysis is required to obtain the expression profiles of all transcripts and identify a greater number of novel genes related to leaf development. However, our expression data are indispensable for elucidating the global transcriptome and true nature of individual gene expression levels in rice due to their compatibility with various datasets in the expression profile database RiceXPro (http:// ricexpro.dna.affrc.go.jp) [26, 27].

Conclusions

We identified a large number of genes that exhibit localized and unique expression patterns during rice leaf development (Figs. 4 and 5; Supplemental Figure 8; Supplemental Table 2). Due to the recent development of CRISPR/Cas9 technology, any genes in the rice genome can be easily knocked out [76]. Thus, reverse genetics strategies to reveal gene function can be applied to the genes identified in this study. Our findings will provide the foundation for future research on the development of grass leaves and contribute to the elucidation of genetic programs unique to grasses.

Methods

Plant materials and growth conditions

Rice (Oryza sativa L. ssp. japonica cv. Nipponbare) seeds were obtained from the National Agriculture and Food Research Organization. They were sown in germination boxes, and the seedlings were grown in a growth chamber (14-h light period at 30 °C and 10-h dark period at 25 °C). At 14–15 days after germination, seedlings in which the tip of the fifth leaf had just emerged from the fourth leaf were collected. The seedlings were dissected under a dissecting microscope to separate tissues at different developmental stages and locations along the longitudinal axis of the leaves. The dissected tissues were used for RNA extraction of microarray analysis.

RNA extraction and microarray analysis

Total RNA was extracted from the collected samples using a RNeasy Mini Kit (Qiagen, Hilden, Germany), and labeling was performed using Quick Amp Labeling Kit, One-Color (Agilent Technologies) in the presence of cyanine-3 (Cy3)-CTP according to the manufacturer’s protocol. The resulting Cy3-labeled cRNA was purified using a RNeasy Mini Kit (Qiagen, Hilden, Germany). A total of 1000 ng Cy3-labeled cRNA was fragmented and hybridized onto a slide of the rice 4 × 44 K microarray RAP-DB (G2519F#15241; Agilent Technologies). After washing, the slides were scanned on an Agilent G2505B DNA microarray scanner, and background correction of the raw Cy3 signals was performed using Feature Extraction 10.5.1.1 software (Agilent Technologies).

Statistical analysis

The processed raw signal intensities of all probes were transformed to log2 scale, and normalization procedures were performed, including 75th-percentile normalization for inter-array comparison and baseline correlation of each probe to the median value of all the samples, using GeneSpring GX12 software (Agilent Technologies). In this study, we used 31,996 probes corresponding to 24,022 loci, which had raw signal intensities > 50 in at least one of the 36 microarray datasets. A density plot of probes for each sample was drawn using ggridges package in R [77]. Pearson correlation analysis was performed on all samples using the ggcorrplot package in R [78]. PCA based on the variance-covariance matrix was first performed on all samples using the prcomp function in R. Subsequently, PCA was applied again to all samples except for P4Bm, using the PC1, PC2, and PC3 scores as variables to compute modified principal components (mPC1, mPC2, and mPC3). The distribution of the samples in the space defined by mPC1 and mPC2 was approximated with a quadratic formula using the nls function in R. One-way ANOVA with Benjamini and Hochberg FDR correction was applied to detect the probes that were differentially expressed among samples (FDR-adjusted p-value = 0.001) using GeneSpring GX12. The extracted probes were classified into 28 clusters through K-means cluster analysis based on Euclidean distance using the cclust package in R [79]. GO enrichment analysis was performed on the genes in each cluster using a tool based on a rice gene coexpression database, RiceFREND [80]. GO terms corresponding to fewer than five genes were discarded. As GO terms are organized in a hierarchical structure, only the child terms were retained. Transcription factor and transcriptional regulator genes were extracted from each cluster using the list in the Plant Transcription Factor Database [81]. In addition, the enrichment of each gene family in each cluster was statistically tested by comparing the number of family members in each cluster with that obtained from the rice 44 K microarray platform, using one-sided Fisher’s exact tests with Benjamini and Hochberg FDR correction with the fisher.test and p.adjust functions in R. The expression profiles shown in Supplemental Figure 3, 5 and Fig. 3, and the heatmaps in Supplemental Figsure 6 and 7 were drawn using the ggplot2 package in R [82].

Gene annotation used for analyses

The gene annotation used in this study was obtained from RAP-DB [29, 30]. The list of genes associated with the biosynthesis, catabolism, and signaling of six phytohormones was derived from Hirano et al. (2008) [83]. For the extraction of transcription factor and transcriptional regulator genes from the rice genome, the list from the Plant Transcription Factor Database [81] was used.

In situ hybridization

Rice (O. sativa ssp. japonica cv. Nipponbare) seeds were obtained from the National Agriculture and Food Research Organization. They were sown in soil, and the seedlings were grown in a growth chamber (14-h light period at 30 °C and 10-h dark period at 25 °C). At 14–15 days after germination, the shoot apex of each four-leaf stage seedling was dissected and fixed in 4% paraformaldehyde in 0.1 M sodium phosphate buffer for 24 h at 4 °C, and then dehydrated in a graded ethanol series. The ethanol was replaced with Histo-Clear (National Diagnostics, Atlanta, GA), and the samples were embedded in Paraplast Plus (Leica, Wetzlar, Germany). Paraffin sections (thickness, 8 μm) were placed on microscope slides and coated with 3-aminopropyl triethoxysilane (Matsunami Glass, Osaka, Japan). To generate probes for 30 genes (see Supplemental Table 3), the corresponding full-length cDNA clones were obtained from NIAS Genebank [84] and used as templates. For 20 genes, cDNA fragments were amplified using PCR and cloned into the pCR Blunt II TOPO vector (Invitrogen, Carlsbad, CA) using primers specific to each gene (Supplemental Table 3). The cDNA was amplified via PCR or digested with restriction enzymes and then transcribed using digoxigenin-labeled antisense riboprobes using a MAXIscript In Vitro Transcription kit (Life Technologies, Carlsbad, CA) with digoxigenin-11-UTP (Roche, Basel, Switzerland), or using T7, SP6 (Takara, Shiga, Japan), or T3 (Roche) RNA polymerase with DIG-RNA labeling mix (Roche). For detecting mature miR167, a 3′-end digoxigenin-labeled probe with the sequence 5′-AgAtCaTgCtGgCaGcTtCa-3′ was used, where uppercase and lowercase letters represent BNANC[NMe] and DNA, respectively. In situ hybridization and immunological detection of the hybridization signals were performed as described by Kouchi and Hata (1993) [85].

Real-time RT-PCR analysis

For verification of expression profiles, the same RNA samples from the four parts of P4 stage: P4S, P4Bb, P4Bm and P4Ba (Table 1) as those used for microarray were examined. For examination of temporal expression profiles, RNA samples from 1 cm of the most basal part of rice shoots at the second, third and fourth leaf stage were used. RNA was treated with TURBO DNA-free™ Kit (Invitrogen, Carlsbad, CA), and cDNA was synthesized using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA). Quantitative RT- PCR was performed with the StepOne™ Real-Time PCR System (Life Technologies, Carlsbad, CA) using the TaqMan Fast Universal PCR Master Mix and FAM-labeled TaqMan probes for each gene. OsRAD6 (Os03g0791800) was used as an internal standard [31]. Three technical replicates were performed for each sample. The primers and TaqMan probes for each gene are listed in Supplemental Table S4.

Availability of data and materials

The microarray datasets generated and analyzed during the current study are available in the Gene Expression Omnibus (GEO) repository through accession number GSE159047 [https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/geo/query/acc.cgi?acc=GSE159047]. The list of transcription factor and transcriptional regulator genes are available in the Plant Transcription Factor Database [http://plntfdb.bio.uni-potsdam.de/v3.0/].

Abbreviations

SPL :

SQUAMOSA PROMOTER BINDING-LIKE

SAM:

Shoot apical meristem

GRF :

GROWTH-REGULATING FACTORS

GIF :

GRF-INTERACTING FACTORS

Pn:

Plastchron number

PCA:

Principal component analysis

ANOVA:

Analysis of variance

FDR:

False discovery rate

GO:

Gene ontology

PIF :

Phytochrome-interacting bHLH factors

BNANC :

Bridged nucleic acid

LG1 :

LIGULELESS1

BR:

Brassinosteroid

RT-PCR:

Reverse transcription PCR

GEO:

Gene Expression Omnibus

References

  1. Tsukaya H. Leaf Development The Arabidopsis Book, vol. 11; 2013. p. e0163.

    Google Scholar 

  2. Moreno MA, Harper LC, Krueger RW, Dellaporta SL, Freeling M. liguleless1 encodes a nuclear-localized protein required for induction of ligules and auricles during maize leaf organogenesis. Genes Dev. 1997;11(5):616–28.

    Article  CAS  PubMed  Google Scholar 

  3. Lee J, Park JJ, Kim SL, Yim J, An G. Mutations in the rice liguleless gene result in a complete loss of the auricle, ligule, and laminar joint. Plant Mol Biol. 2007;65(4):487–99.

    Article  CAS  PubMed  Google Scholar 

  4. Liu K, Cao J, Yu K, Liu X, Gao Y, Chen Q, Zhang W, Peng H, Du J, Xin M, Hu Z, Guo W, Rossi V, Ni Z, Sun Q, Yao Y. Wheat TaSPL8 modulates leaf angle through Auxin and Brassinosteroid signaling. Plant Physiol. 2019;181(1):179–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Toriba T, Tokunaga H, Shiga T, Nie F, Naramoto S, Honda E, Tanaka K, Taji T, Itoh JI, Kyozuka J. BLADE-ON-PETIOLE genes temporally and developmentally regulate the sheath to blade ratio of rice leaves. Nat Commun. 2019;10(1):619.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Jackson D. Double labeling of KNOTTED1 mRNA and protein reveals multiple potential sites of protein trafficking in the shoot apex. Plant Physiol. 2002;129(4):1423–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Ramirez J, Bolduc N, Lisch D, Hake S. Distal expression of knotted1 in maize leaves leads to reestablishment of proximal/distal patterning and leaf dissection. Plant Physiol. 2009;151(4):1878–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Liu H, Guo S, Xu Y, Li C, Zhang Z, Zhang D, Xu S, Zhang C, Chong K. OsmiR396d-regulated OsGRFs function in floral organogenesis in rice through binding to their targets OsJMJ706 and OsCR4. Plant Physiol. 2014;165(1):160–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. He Z, Zeng J, Ren Y, Chen D, Li W, Gao F, Cao Y, Luo T, Yuan G, Wu X, Liang Y, Deng Q, Wang S, Zheng A, Zhu J, Liu H, Wang L, Li P, Li S. OsGIF1 positively regulates the sizes of stems, leaves, and grains in Rice. Front Plant Sci. 2017;8:1730.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Zhang D, Sun W, Singh R, Zheng Y, Cao Z, Li M, Lunde C, Hake S, Zhang Z. GRF-interacting factor1 regulates shoot architecture and meristem determinacy in maize. Plant Cell. 2018;30(2):360–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Shimano S, Hibara KI, Furuya T, Arimura SI, Tsukaya H, Itoh JI. Conserved functional control, but distinct regulation, of cell proliferation in rice and Arabidopsis leaves revealed by comparative analysis of GRF-INTERACTING FACTOR 1 orthologs. Development. 2018;145(7):dev159624.

    Article  PubMed  CAS  Google Scholar 

  12. Nelissen H, Eeckhout D, Demuynck K, Persiau G, Walton A, van Bel M, Vervoort M, Candaele J, De Block J, Aesaert S, Van Lijsebettens M, Goormachtig S, Vandepoele K, Van Leene J, Muszynski M, Gevaert K, Inzé D, De Jaeger G. Dynamic changes in ANGUSTIFOLIA3 complex composition reveal a growth regulatory mechanism in the maize leaf. Plant Cell. 2015;27(6):1605–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Schmid M, Davison TS, Henz SR, Pape UJ, Demar M, Vingron M, Schölkopf B, Weigel D, Lohmann JU. A gene expression map of Arabidopsis thaliana development. Nat Genet. 2005;37(5):501–6.

    Article  CAS  PubMed  Google Scholar 

  14. Klepikova AV, Kasianov AS, Gerasimov ES, Logacheva MD, Penin AA. A high resolution map of the Arabidopsis thaliana developmental transcriptome based on RNA-seq profiling. Plant J. 2016;88(6):1058–70.

    Article  CAS  PubMed  Google Scholar 

  15. Brooks L 3rd, Strable J, Zhang X, Ohtsu K, Zhou R, Sarkar A, Hargreaves S, Elshire RJ, Eudy D, Pawlowska T, Ware D, Janick-Buckner D, Buckner B, Timmermans MC, Schnable PS, Nettleton D, Scanlon MJ. Microdissection of shoot meristem functional domains. PLoS Genet. 2009;5(5):e1000476.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Li P, Ponnala L, Gandotra N, Wang L, Si Y, Tausta SL, Kebrom TH, Provart N, Patel R, Myers CR, Reidel EJ, Turgeon R, Liu P, Sun Q, Nelson T, Brutnell TP. The developmental dynamics of the maize leaf transcriptome. Nat Genet. 2010;42(12):1060–7.

    Article  CAS  PubMed  Google Scholar 

  17. Sekhon RS, Lin H, Childs KL, Hansey CN, Buell CR, de Leon N, Kaeppler SM. Genome-wide atlas of transcription during maize development. Plant J. 2011;66(4):553–63.

    Article  CAS  PubMed  Google Scholar 

  18. Wang P, Kelly S, Fouracre JP, Langdale JA. Genome-wide transcript analysis of early maize leaf development reveals gene cohorts associated with the differentiation of C4 Kranz anatomy. Plant J. 2013;75(4):656–70.

    Article  CAS  PubMed  Google Scholar 

  19. Liu WY, Chang YM, Chen SC, Lu CH, Wu YH, Lu MY, Chen DR, Shih AC, Sheue CR, Huang HC, Yu CP, Lin HH, Shiu SH, Ku MS, Li WH. Anatomical and transcriptional dynamics of maize embryonic leaves during seed germination. Proc Natl Acad Sci U S A. 2013;110(10):3979–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Yu CP, Chen SC, Chang YM, Liu WY, Lin HH, Lin JJ, Chen HJ, Lu YJ, Wu YH, Lu MY, Lu CH, Shih AC, Ku MS, Shiu SH, Wu SH, Li WH. Transcriptome dynamics of developing maize leaves and genomewide prediction of cis elements and their cognate transcription factors. Proc Natl Acad Sci U S A. 2015;112(19):E2477–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Dong L, Qin L, Dai X, Ding Z, Bi R, Liu P, Chen Y, Brutnell TP, Wang X, Li P. Transcriptomic analysis of leaf sheath maturation in maize. Int J Mol Sci. 2019;20(10):E2472.

    Article  PubMed  CAS  Google Scholar 

  22. Knauer S, Javelle M, Li L, Li X, Ma X, Wimalanathan K, Kumari S, Johnston R, Leiboff S, Meeley R, Schnable PS, Ware D, Lawrence-Dill C, Yu J, Muehlbauer GJ, Scanlon MJ, Timmermans MCP. A high-resolution gene expression atlas links dedicated meristem genes to key architectural traits. Genome Res. 2019;29(12):1962–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Jiao Y, Tausta SL, Gandotra N, Sun N, Liu T, Clay NK, Ceserani T, Chen M, Ma L, Holford M, Zhang HY, Zhao H, Deng XW, Nelson T. A transcriptome atlas of rice cell types uncovers cellular, functional and developmental hierarchies. Nat Genet. 2009;41(2):258–63.

    Article  CAS  PubMed  Google Scholar 

  24. van Campen JC, Yaapar MN, Narawatthana S, Lehmeier C, Wanchana S, Thakur V, Chater C, Kelly S, Rolfe SA, Quick WP, Fleming AJ. Combined chlorophyll fluorescence and Transcriptomic analysis identifies the P3/P4 transition as a key stage in Rice leaf photosynthetic development. Plant Physiol. 2016;170(3):1655–74.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Zhou LJ, Xiao LT, Xue HW. Dynamic cytology and transcriptional regulation of Rice Lamina joint development. Plant Physiol. 2017;174(3):1728–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Sato Y, Antonio BA, Namiki N, Takehisa H, Minami H, Kamatsuki K, Sugimoto K, Shimizu Y, Hirochika H, Nagamura Y. RiceXPro: a platform for monitoring gene expression in japonica rice grown under natural field conditions. Nucleic Acids Res. 2011;39:D1141–8.

    Article  CAS  PubMed  Google Scholar 

  27. Sato Y, Takehisa H, Kamatsuki K, Minami H, Namiki N, Ikawa H, Ohyanagi H, Sugimoto K, Antonio BA, Nagamura Y. RiceXPro version 3.0: expanding the informatics resource for rice transcriptome. Nucleic Acids Res. 2013;41:D1206–13.

    Article  CAS  PubMed  Google Scholar 

  28. Itoh J, Nonomura K, Ikeda K, Yamaki S, Inukai Y, Yamagishi H, Kitano H, Nagato Y. Rice plant development: from zygote to spikelet. Plant Cell Physiol. 2005;46(1):23–47.

    Article  CAS  PubMed  Google Scholar 

  29. Rice Annotation Project. The Rice annotation project database (RAP-DB): 2008 update. Nucleic Acids Res. 2008;36:D1028–33.

    Article  CAS  Google Scholar 

  30. Sakai H, Lee SS, Tanaka T, Numa H, Kim J, Kawahara Y, Wakimoto H, Yang CC, Iwamoto M, Abe T, Yamada Y, Muto A, Inokuchi H, Ikemura T, Matsumoto T, Sasaki T, Itoh T. Rice annotation project database (RAP-DB): an integrative and interactive database for rice genomics. Plant Cell Physiol. 2013;54(2):e6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Dechkrong P, Yoshikawa T, Itoh JI. Morphological and molecular dissection of leaf development in wild-type and various morphogenetic mutants in Rice. Am J Plant Sci. 2015;6:1215–32.

    Article  Google Scholar 

  32. Tsuda K, Ito Y, Sato Y, Kurata N. Positive autoregulation of a KNOX gene is essential for shoot apical meristem maintenance in rice. Plant Cell. 2011;23(12):4368–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Hibara K, Nagato Y. OsNAM and OsCUC3 are expressed specifically in organ boundaries. Rice Genet Newslett. 2006;23:96–7.

    Google Scholar 

  34. Horstman A, Willemsen V, Boutilier K, Heidstra R. AINTEGUMENTA-LIKE proteins: hubs in a plethora of networks. Trends Plant Sci. 2014;19(3):146–57.

    Article  CAS  PubMed  Google Scholar 

  35. Li P, Xue H. Structural characterization and expression pattern analysis of the rice PLT gene family. Acta Biochim Biophys Sin Shanghai. 2011;43(9):688–97.

    Article  CAS  PubMed  Google Scholar 

  36. Takeda T, Suwa Y, Suzuki M, Kitano H, Ueguchi-Tanaka M, Ashikari M, Matsuoka M, Ueguchi C. The OsTB1 gene negatively regulates lateral branching in rice. Plant J. 2003;33(3):513–20.

    Article  CAS  PubMed  Google Scholar 

  37. Lu D, Wang T, Persson S, Mueller-Roeber B, Schippers JH. Transcriptional control of ROS homeostasis by KUODA1 regulates cell expansion during leaf development. Nat Commun. 2014;7(5):3767.

    Article  CAS  Google Scholar 

  38. Raggi S, Ferrarini A, Delledonne M, Dunand C, Ranocha P, De Lorenzo G, Cervone F, Ferrari S. The Arabidopsis class III peroxidase AtPRX71 negatively regulates growth under physiological conditions and in response to Cell Wall damage. Plant Physiol. 2015;169(4):2513–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Franková L, Fry SC. Biochemistry and physiological roles of enzymes that 'cut and paste' plant cell-wall polysaccharides. J Exp Bot. 2013;64(12):3519–50.

    Article  PubMed  CAS  Google Scholar 

  40. Gangappa SN, Botto JF. The BBX family of plant transcription factors. Trends Plant Sci. 2014;19(7):460–70.

    Article  CAS  PubMed  Google Scholar 

  41. Nakamura Y, Kato T, Yamashino T, Murakami M, Mizuno T. Characterization of a set of phytochrome-interacting factor-like bHLH proteins in Oryza sativa. Biosci Biotechnol Biochem. 2007;71(5):1183–91.

    Article  CAS  PubMed  Google Scholar 

  42. Hirose T, Endler A, Ohsugi R. Gene expression of enzymes for starch and sucrose metabolism and transport in leaf sheaths of Rice (Oryza sativa L.) during the heading period in relation to the sink to source transition. Plant Production Science. 1999;2(3):178–83.

    Article  Google Scholar 

  43. Luo X, Zheng J, Huang R, Huang Y, Wang H, Jiang L, Fang X. Phytohormones signaling and crosstalk regulating leaf angle in rice. Plant Cell Rep. 2016;35(12):2423–33.

    Article  PubMed  CAS  Google Scholar 

  44. Ruan W, Guo M, Xu L, Wang X, Zhao H, Wang J, Yi K. An SPX-RLI1 module regulates leaf inclination in response to phosphate availability in Rice. Plant Cell. 2018;30(4):853–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Kim JH, Kim J, Jun SE, Park S, Timilsina R, Kwon DS, Kim Y, Park SJ, Hwang JY, Nam HG, Kim GT, Woo HR. ORESARA15, a PLATZ transcription factor, mediates leaf growth and senescence in Arabidopsis. New Phytol. 2018;220(2):609–23.

    Article  CAS  PubMed  Google Scholar 

  46. Finet C, Berne-Dedieu A, Scutt CP, Marlétaz F. Evolution of the ARF gene family in land plants: old domains, new tricks. Mol Biol Evol. 2013;30(1):45–56.

    Article  CAS  PubMed  Google Scholar 

  47. Majeran W, Friso G, Ponnala L, Connolly B, Huang M, Reidel E, Zhang C, Asakura Y, Bhuiyan NH, Sun Q, Turgeon R, van Wijk KJ. Structural and metabolic transitions of C4 leaf development and differentiation defined by microscopy and quantitative proteomics in maize. Plant Cell. 2010;22(11):3509–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Pick TR, Bräutigam A, Schlüter U, Denton AK, Colmsee C, Scholz U, Fahnenstich H, Pieruschka R, Rascher U, Sonnewald U, Weber AP. Systems analysis of a maize leaf developmental gradient redefines the current C4 model and provides candidates for regulation. Plant Cell. 2011;23(12):4208–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Nelissen H, Rymen B, Jikumaru Y, Demuynck K, Van Lijsebettens M, Kamiya Y, Inzé D, Beemster GT. A local maximum in gibberellin levels regulates maize leaf growth by spatial control of cell division. Curr Biol. 2012;22(13):1183–7.

    Article  CAS  PubMed  Google Scholar 

  50. Nelissen H, Sun XH, Rymen B, Jikumaru Y, Kojima M, Takebayashi Y, Abbeloos R, Demuynck K, Storme V, Vuylsteke M, De Block J, Herman D, Coppens F, Maere S, Kamiya Y, Sakakibara H, Beemster GTS, Inzé D. The reduction in maize leaf growth under mild drought affects the transition between cell division and cell expansion and cannot be restored by elevated gibberellic acid levels. Plant Biotechnol J. 2018;16(2):615–27.

    Article  CAS  Google Scholar 

  51. Facette MR, Shen Z, Björnsdóttir FR, Briggs SP, Smith LG. Parallel proteomic and phosphoproteomic analyses of successive stages of maize leaf development. Plant Cell. 2013;25(8):2798–812.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Tausta SL, Li P, Si Y, Gandotra N, Liu P, Sun Q, Brutnell TP, Nelson T. Developmental dynamics of Kranz cell transcriptional specificity in maize leaf reveals early onset of C4-related processes. J Exp Bot. 2014;65(13):3543–55.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Wang L, Czedik-Eysenberg A, Mertz RA, Si Y, Tohge T, Nunes-Nesi A, Arrivault S, Dedow LK, Bryant DW, Zhou W, Xu J, Weissmann S, Studer A, Li P, Zhang C, LaRue T, Shao Y, Ding Z, Sun Q, Patel RV, Turgeon R, Zhu X, Provart NJ, Mockler TC, Fernie AR, Stitt M, Liu P, Brutnell TP. Comparative analyses of C4 and C3 photosynthesis in developing leaves of maize and rice. Nat Biotechnol. 2014;32(11):1158–65.

    Article  PubMed  CAS  Google Scholar 

  54. Czedik-Eysenberg A, Arrivault S, Lohse MA, Feil R, Krohn N, Encke B, Nunes-Nesi A, Fernie AR, Lunn JE, Sulpice R, Stitt M. The interplay between carbon availability and growth in different zones of the growing maize leaf. Plant Physiol. 2016;172(2):943–67.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. Nelissen H, Gonzalez N, Inzé D. Leaf growth in dicots and monocots: so different yet so alike. Curr Opin Plant Biol. 2016;33:72–6.

    Article  PubMed  Google Scholar 

  56. Yamaguchi T, Nagasawa N, Kawasaki S, Matsuoka M, Nagato Y, Hirano HY. The YABBY gene DROOPING LEAF regulates carpel specification and midrib development in Oryza sativa. Plant Cell. 2004;16(2):500–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Strable J, Wallace JG, Unger-Wallace E, Briggs S, Bradbury PJ, Buckler ES, Vollbrecht E. Maize YABBY genes drooping leaf1 and drooping leaf2 regulate plant architecture. Plant Cell. 2017;29(7):1622–41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Furutani I, Sukegawa S, Kyozuka J. Genome-wide analysis of spatial and temporal gene expression in rice panicle development. Plant J. 2006;46(3):503–11.

    Article  CAS  PubMed  Google Scholar 

  59. Yadav RK, Girke T, Pasala S, Xie M, Reddy GV. Gene expression map of the Arabidopsis shoot apical meristem stem cell niche. Proc Natl Acad Sci U S A. 2009;106(12):4941–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Takacs EM, Li J, Du C, Ponnala L, Janick-Buckner D, Yu J, Muehlbauer GJ, Schnable PS, Timmermans MC, Sun Q, Nettleton D, Scanlon MJ. Ontogeny of the maize shoot apical meristem. Plant Cell. 2012;24(8):3219–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Li G, Wang D, Yang R, Logan K, Chen H, Zhang S, Skaggs MI, Lloyd A, Burnett WJ, Laurie JD, Hunter BG, Dannenhoffer JM, Larkins BA, Drews GN, Wang X, Yadegari R. Temporal patterns of gene expression in developing maize endosperm identified through transcriptome sequencing. Proc Natl Acad Sci U S A. 2014;111(21):7582–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Slane D, Kong J, Berendzen KW, Kilian J, Henschen A, Kolb M, Schmid M, Harter K, Mayer U, De Smet I, Bayer M, Jürgens G. Cell type-specific transcriptome analysis in the early Arabidopsis thaliana embryo. Development. 2014;141(24):4831–40.

    Article  CAS  PubMed  Google Scholar 

  63. Johnston R, Wang M, Sun Q, Sylvester AW, Hake S, Scanlon MJ. Transcriptomic analyses indicate that maize ligule development recapitulates gene expression patterns that occur during lateral organ initiation. Plant Cell. 2014;26(12):4718–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Itoh J, Sato Y, Sato Y, Hibara K, Shimizu-Sato S, Kobayashi H, Takehisa H, Sanguinet KA, Namiki N, Nagamura Y. Genome-wide analysis of spatiotemporal gene expression patterns during early embryogenesis in rice. Development. 2016;143(7):1217–27.

    Article  CAS  PubMed  Google Scholar 

  65. Feng N, Song G, Guan J, Chen K, Jia M, Huang D, Wu J, Zhang L, Kong X, Geng S, Liu J, Li A, Mao L. Transcriptome profiling of wheat inflorescence development from spikelet initiation to floral patterning identified stage-specific regulatory genes. Plant Physiol. 2017;174(3):1779–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Kir G, Ye H, Nelissen H, Neelakandan AK, Kusnandar AS, Luo A, Inzé D, Sylvester AW, Yin Y, Becraft PW. RNA interference knockdown of BRASSINOSTEROID INSENSITIVE1 in maize reveals novel functions for Brassinosteroid signaling in controlling plant architecture. Plant Physiol. 2015;169(1):826–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Jang S, An G, Li HY. Rice leaf angle and grain size are affected by the OsBUL1 transcriptional activator complex. Plant Physiol. 2017;173(1):688–702.

    Article  CAS  PubMed  Google Scholar 

  68. Matsukura C, Kawai M, Toyofuku K, Barrero RA, Uchimiya H, Yamaguchi J. Transverse vein differentiation associated with gas space formation - fate of the middle cell layer in leaf sheath development of Rice. Ann Bot. 2000;85(1):19–27.

    Article  Google Scholar 

  69. Yamauchi T, Tanaka A, Inahashi H, Nishizawa NK, Tsutsumi N, Inukai Y, Nakazono M. Fine control of aerenchyma and lateral root development through AUX/IAA- and ARF-dependent auxin signaling. Proc Natl Acad Sci U S A. 2019;116(41):20770–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Kaufman PB, Brock TG, Song I, Rho YB, Ghosheh NS. How cereal grass shoots perceive and respond to gravity. Am J Bot. 1987;74(9):1446–57.

    Article  CAS  PubMed  Google Scholar 

  71. Yoshihara T, Iino M. Identification of the gravitropism-related rice gene LAZY1 and elucidation of LAZY1-dependent and -independent gravity signaling pathways. Plant Cell Physiol. 2007;48(5):678–88.

    Article  CAS  PubMed  Google Scholar 

  72. Tan L, Li X, Liu F, Sun X, Li C, Zhu Z, Fu Y, Cai H, Wang X, Xie D, Sun C. Control of a key transition from prostrate to erect growth in rice domestication. Nature genet. 2008;40(11):1360–4.

    Article  CAS  PubMed  Google Scholar 

  73. Wu X, Tang D, Li M, Wang K, Cheng Z. Loose plant Architecture1, an INDETERMINATE DOMAIN protein involved in shoot gravitropism, regulates plant architecture in rice. Plant Physiol. 2013;161(1):317–29.

    Article  CAS  PubMed  Google Scholar 

  74. Wu MF, Tian Q, Reed JW. Arabidopsis microRNA167 controls patterns of ARF6 and ARF8 expression, and regulates both female and male reproduction. Development. 2006;133(21):4211–8.

    Article  CAS  PubMed  Google Scholar 

  75. Ágyi Á, Havelda Z. Analysis of gradient-like expression of miR167 in Arabidopsis thaliana embryonic tissue. J Plant Biol. 2013;56:336–44.

    Article  CAS  Google Scholar 

  76. Endo M, Mikami M, Toki S. Multigene knockout utilizing off-target mutations of the CRISPR/Cas9 system in rice. Plant Cell Physiol. 2015;56(1):41–7.

    Article  CAS  PubMed  Google Scholar 

  77. Claus O. Wilke (2020). ggridges: Ridgeline Plots in 'ggplot2'. R package version 0.5.2. https://CRAN.R-project.org/package=ggridges. Accessed 2 Jan 2021.

  78. Kassambara A. ggcorrplot: Visualization of a Correlation Matrix using 'ggplot2'. R package version 0.1.3. https://CRAN.R-project.org/package=ggcorrplot. Accessed 1 Oct 2020.

  79. Dimitriadou E. cclust: Convex Clustering Methods and Clustering Indexes. R package version 0.6–21. https://CRAN.R-project.org/package=cclust. Accessed 1 Oct 2020.

  80. Sato Y, Namiki N, Takehisa H, Kamatsuki K, Minami H, Ikawa H, Ohyanagi H, Sugimoto K, Itoh J, Antonio BA, Nagamura Y. RiceFREND: a platform for retrieving coexpressed gene networks in rice. Nucleic Acids Res. 2013;41:D1214–21 http://ricefrend.dna.affrc.go.jp/.

    Article  CAS  PubMed  Google Scholar 

  81. Pérez-Rodríguez P, Riaño-Pachón DM, Corrêa LG, Rensing SA, Kersten B, Mueller-Roeber B. PlnTFDB: updated content and new features of the plant transcription factor database. Nucleic Acids Res. 2010;38:D822–7 http://plntfdb.bio.uni-potsdam.de/v3.0/.

    Article  PubMed  CAS  Google Scholar 

  82. Wickham H. ggplot2: elegant graphics for data analysis. Springer-Verlag New York ISBN 978-3-319-24277-4. https://ggplot2.tidyverse.org. Accessed 1 Oct 2020.

  83. Hirano K, Aya K, Hobo T, Sakakibara H, Kojima M, Shim RA, Hasegawa Y, Ueguchi-Tanaka M, Matsuoka M. Comprehensive transcriptome analysis of phytohormone biosynthesis and signaling genes in microspore/pollen and tapetum of rice. Plant Cell Physiol. 2008;49(10):1429–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. NIAS Genebank. http://www.dna.affrc.go.jp/distribution/ Accessed 1 Oct 2020.

  85. Kouchi H, Hata S. Isolation and characterization of novel nodulin cDNAs representing genes expressed at early stages of soybean nodule development. Mol Gen Genet. 1993;238:106–19.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank Dr. Yoshiaki Nagamura and Ms. Ritsuko Motoyama (NARO) for helping the design and support for microarray analysis, respectively. We also thank Dr. Takeshi Izawa (The University of Tokyo) and Dr. Kenichiro Hibara (Kibi International University) for helpful comments, and Ryuichi Soga and Hidetoshi Teshima (Institute for Sustainable Agro-ecosystem Services, The University of Tokyo) for their assistance in cultivating rice plants.

Funding

This work was supported partly by a project of the Ministry of Agriculture, Forestry and Fisheries of Japan (Genomics for Agricultural Innovation, RTR0002, RTR0006), and partly by a grant from Japan Society for the Promotion of Science, KAKENHI [20H02955 and 20H05406 to JI].

Author information

Authors and Affiliations

  1. Graduate School of Agricultural and Life Sciences, University of Tokyo, Tokyo, 113-8657, Japan

    Masayuki Miya & Jun-Ichi Itoh

  2. Graduate School of Agriculture, Kyoto University, Kyoto, 606-8502, Japan

    Takanori Yoshikawa

  3. Institute of Crop Science, National Agriculture and Food Research Organization (NARO), Tsukuba, 305-8518, Japan

    Yutaka Sato

Authors
  1. Masayuki Miya
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Takanori Yoshikawa
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yutaka Sato
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jun-Ichi Itoh
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

MM and JI designed experiments; MM, TY, YS and JI performed experiments; MM and JI wrote the article. All authors have read and approved the manuscript.

Corresponding author

Correspondence to Jun-Ichi Itoh.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

No conflicts of interest declared.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplemental Figure S1.

Distribution of the normalized intensity values of probes for each sample. Supplemental Figure S2. Heatmap showing Pearson correlation coefficient (PCC) values representing the relationships between samples. Supplemental Figure S3. Verification of expression profiles of six genes by real-time RT-PCR analysis. Supplemental Figure S4. Principal component analysis score plot of samples based on the original principal components. Supplemental Figure S5. Expression patterns of 28 clusters obtained through K-means analysis. Supplemental Figure S6. Heatmap of Gene Ontology (GO) terms overrepresented in each cluster. Supplemental Figure S7. Numbers of transcription factors and transcriptional regulators in each cluster.

Additional file 2: Supplemental Figure S8.

In situ hybridization of genes showing localized expression during leaf development. Supplemental Table 2. Annotation of genes in Fig. 4, 5. Supplemental Figure S9. Spatial expression pattern of LIGULELESS1.

Additional file 3:.

Three-dimensional model of Supplemental Figure 4. A three-dimensional model of PCA score plot of samples based on the original principal components. The proportions of the total variance explained by PC1, PC2, and PC3 are shown in parentheses. Samples collected at the same stage are shown in the same color. Samples with different tissue identities are indicated by different symbols: shoot apex, square; P3 leaf, circle; blade, triangle; blade-sheath boundary, diamond; sheath, inverted triangle. Red arrows represent the directions of the modified principal components (mPC1, mPC2, and mPC3) shown in Fig. 2.

Additional file 4:.

Three-dimensional model of Fig. 2. A three-dimensional model of PCA score plot of samples based on the modified principal components. The proportions of the total variance explained by PC1, PC2, and PC3 are shown in parentheses. Samples collected at the same stage are shown in the same color. Samples with different tissue identities are indicated by different symbols: shoot apex, square; P3 leaf, circle; blade, triangle; blade-sheath boundary, diamond; sheath, inverted triangle. Red arrows represent the directions of the original principal components (PC1, PC2, and PC3) shown in Supplemental Figure 4.

Additional file 5: Supplemental Table S1

. The list of genes categorized into 28 clusters by K-means analysis.

Additional file 6: Supplemental Table S3

. The list of full-length cDNA clones and cloning primers used for preparation of probes for in situ hybridization.

Additional file 7: Supplemental Table S4

. The list of primers and TaqMan probes used for real-time RT-PCR analysis.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Miya, M., Yoshikawa, T., Sato, Y. et al. Genome-wide analysis of spatiotemporal expression patterns during rice leaf development. BMC Genomics 22, 169 (2021). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-021-07494-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-021-07494-5

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us