Introduction
Rice (Oryza sativa L.) is one of the most important staple food crops all over the world, particularly in Asian countries. More than half of the world's population consume rice as a primary source of food and over 3.5 billion people are solely dependent upon rice for at least 20% of their daily required calories (Khush, Reference Khush2013). Under the condition of limited farmland, the main challenge to meet the food demand is increasing rice yield per unit area which leads to the improvement in overall rice yields.
Leaf is the major photosynthetic organ in rice which is directly related to biomass accumulation and grain yield production (Alamin et al., Reference Alamin, Zeng, Qin, Sultana, Jin and Shi2017). It has been proved that appropriate leaf shape improves photosynthesis rate resulting in the higher grain yield in rice (Wu et al., Reference Wu, Fu, Hu, Si, Cheng and Liu2010). As a result, leaf morphology is becoming one of the key concerns in the study of high yield plant breeding, especially high yield rice breeding (Li et al., Reference Li, Xue, Chen, Zhang, Ma, Pan, Zhu, Pan and Zuo2014). The flag leaf or uppermost three leaves were proposed to be long, narrow, moderately rolled (V-shaped), erect and thick for the hybrid rice to be super high yielding (Zhou et al., Reference Zhou, Ma, Liu, Shen and Pan1995; Yuan, Reference Yuan1997; Wu, Reference Wu2009). Particularly, a moderately rolling leaf (RL) of rice can increase the grain yield by improving photosynthetic efficiency and increasing stress resistance through reduced leaf transpiration during drought resulting in the compact tailing and vertical position of the leaf (Lang et al., Reference Lang, Zhang, Gu, Yang and Zhu2004; Zhang et al., Reference Zhang, Xu, Zhu, Qian and Xue2009; Wu et al., Reference Wu, Fu, Hu, Si, Cheng and Liu2010; Zou et al., Reference Zou, Sun, Zhang, Liu, Wu, Tian, Qiu and Lu2011). Thus, identifying moderately RL rice mutants and isolating the responsible genes for controlling leaf rolling will be advantageous for rice breeding with the desired traits (e.g. high yield) and architecture (Zhang et al., Reference Zhang, Xu, Zhu, Qian and Xue2009; Xu et al., Reference Xu, Wang, Long, Huang, Wang, Zhou, Zheng, Sun, Chen and Chen2014; Alamin et al., Reference Alamin, Zeng, Qin, Sultana, Jin and Shi2017).
To date, at least 103 genes/quantitative trait loci (QTLs) for RL trait of rice have been identified by several studies. Most of the existing studies of RL trait of rice are related to the comparison of a RL mutant with a wild type, RL genes/QTLs isolation/identification and investigation of their molecular functions, and cloning/mapping of the identified RL genes/QTLs for developing a rice mutant with desired trait and architecture, which are very time consuming and laborious. So, collective information on RL genes and their comparative analysis are very essential and helpful before starting the experiment and functional studies of a desired RL rice mutant. To our knowledge, there is no collective information on the RL genes reported to date and different comparative bioinformatics analyses of their sequences are still incomplete. Therefore, this bioinformatics study was designed to compare the structures and functions of all RL genes/QTLs, characterized so far throughout several studies, using various bioinformatics analyses including gene structure (exons/introns), conserved domain, phylogenetic, gene ontology (GO), transcription factor (TF), Kyoto Encyclopedia of Genes and Genomes (KEGG) and protein–protein interaction network. Altogether, this study will provide collective information about RL genes and their molecular characteristics till date, and might be helpful for geneticists and rice breeders to develop the rice mutants with desired RL trait and architecture.
Materials and methods
The whole process of this study, from data collection to data analysis, has been illustrated through a schematic diagram in Fig. 1. First, we have searched and collected the RL related genes/QTLs in rice from different sources (published papers, databases and search engines) and then finally selected those genes for which locus identifiers (IDs) was available. After this, we have downloaded the sequences (genomic, CDS and protein sequences) of the finally selected genes and performed various bioinformatics analyses with the sequence data to investigate the structures, functions and diversifications of the selected RL genes. The data collection procedures for RL related genes and the bioinformatics analysis techniques used to analyse the data have been discussed in the next subsections in detail.
Schematic diagram of the study.
Data collection procedures for RL related genes
Identification of RL related genes was done in several ways in this study (Fig. 1). First, we searched and collected all the published papers which were related to RL of different varieties (japonica, indica and hybrid) of rice (Oryza sativa L.). From these collected papers we extracted all the RL related genes. We also used some most popular and widely used publicly available databases such as Oryzabase (http://shigen.nig.ac.jp/rice/oryzabase), GRAMENE (http://archive.gramene.org/), Rice Genome Browser (http://rice.plantbiology.msu.edu/cgi-bin/gbrowse/rice/), Rice Annotation Project (RAP) database (RAP-DB: https://rapdb.dna.affrc.go.jp/), PubMed (https://www.ncbi.nlm.nih.gov/pubmed) and some most popular search engines (e.g. Google Scholar, Google, Yahoo and Bing) to find out the RL genes of rice along with locus ID. In most of the papers that reported locus identifiers (IDs) of RL genes, the Michigan State University (MSU) locus IDs were reported, while some papers (Fujino et al., Reference Fujino, Matsuda, Ozawa, Nishimura, Koshiba, Fraaije and Sekiguchi2008; Segami et al., Reference Segami, Kono, Ando, Yano, Kitano, Miura and Iwasaki2012; Yang et al., Reference Yang, Li, Miao, Gan, Qiao, Chang, Shi and Chen2016) mentioned the RAP locus IDs instead of MSU ID for RL related genes. Also, in some papers, only the names of the RL genes were reported instead of locus ID. MSU/RAP locus IDs of the RL genes, for which locus IDs were not reported, were found using various web databases of rice including Oryzabase (http://shigen.nig.ac.jp/rice/oryzabase), GRAMENE (http://archive.gramene.org/), Rice Genome Browser (http://rice.plantbiology.msu.edu/cgi-bin/gbrowse/rice/) and RAP database (RAP-DB: https://rapdb.dna.affrc.go.jp/). RAP IDs of RL genes were converted to MSU IDs using the RAP database (http://rapdb.dna.affrc.go.jp).
A total of 103 RL related genes/QTLs reported in different kinds of literature and databases were identified in this study. The detailed information of the identified 103 RL genes/QTLs are available in online Supplementary Table S1. Among the total 103 RL genes/QTLs, 9, 11, 16, 9, 7, 5, 14, 3, 12, 9, 2 and 8 genes/QTLs were identified on chromosome 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12, respectively. Note that the genes RL1, RL3 and RL11(t) have been identified on two different chromosomes (1, 12), (3, 12) and (7, 4), respectively. Among the 103 RL genes, ND1 (Li et al., Reference Li, Xiong, Li, Cui, Tang, Zhang, Pauly, Cheng and Zhou2009), DNL3 (Shi et al., Reference Shi, Wei, Adedze, Sheng, Tang, Hu and Wang2016) and NRL1 (Hu et al., Reference Hu, Zhu, Zeng, Gao, Guo, Fang, Zhang, Dong, Yan and Liu2010; Wu et al., Reference Wu, Fu, Hu, Si, Cheng and Liu2010) are identical; RL9 (Yan et al., Reference Yan, Yan, Zeng, Yang, Fang, Tian, Sun, Cheng and Gu2008) and SLL1 (Zhang et al., Reference Zhang, Xu, Zhu, Qian and Xue2009) are identical; and COW1 (Woo et al., Reference Woo, Park, Su'udi, Yang, Park, Back, Park and An2007) and NAL7 (Fujino et al., Reference Fujino, Matsuda, Ozawa, Nishimura, Koshiba, Fraaije and Sekiguchi2008) are identical (online Supplementary Table S1). So, we got a total of 96 distinct RL genes/QTLs after considering one gene from each group of similar genes. Among 96 distinct genes/QTLs, all the QTLs (19 QTLs) were excluded from the analysis since locus ID was not available for any of the QTLs. Then among the remaining 77 genes, we excluded 35 genes for which locus IDs were not found. Finally, we included only 42 RL genes in our analysis for which locus IDs were available. As a result, we confined all of our analyses with only those 42 RL genes: OsACL1 (Li et al., Reference Li, Shi, Li, Shen, Wang, An and Zhang2010), OsADL1 (Hibara et al., Reference Hibara, Obara, Hayashida, Abe, Ishimaru, Satoh, Itoh and Nagato2009), OsAGO1a (Li et al., Reference Li, Xue, Zuo, Chen, Zhang, Li, Zhu, Ma, Pan and Pan2013a), OsAGO7 (Shi et al., Reference Shi, Wang, Wan, Shen, Wang and Zhang2007), OsARF18 (Huang et al., Reference Huang, Li and Zhao2016), OsARVL4 (Wang et al., Reference Wang, Xu, Nian, Shen, Lai, Hu, Zeng, Ge, Fang, Zhu, Qian and Zhang2016), OsAS2 (Ma et al., Reference Ma, Wang, Guo and Zhang2009), OsCFL1 (Wu et al., Reference Wu, Li, He, Waßmann, Yu, Qin, Schreiber, Qu and Gu2011), OsDCL1 (Liu et al., Reference Liu, Li, Li, Liu, Cao, Chu and Cao2005), OsHB4 (Zhang et al., Reference Zhang, Zhang, Srivastava, Pan, Bai, Fang, Shi and Zhu2018), OsI_14279 (Wang et al., Reference Wang, Wang, Chen, Liang, Huang and Yi2017), OsLBD3-7 (Li et al., Reference Li, Zou, Zhang, Shao, Liu, Liu, Li and Zhao2016), OsLC2 (Zhao et al., Reference Zhao, Hu, Guo, Qian and Xue2010), OsLRRK1 (Zhou et al., Reference Zhou, Wang, Wu, Yang, Liu, Yan, Tang, Zhao, Zhu and Lin2018), OsMYB103L (Yang et al., Reference Yang, Li, Liu, Ji, Hao, Zhao, Li, Chen, Cheng and Zhu2014), OsNAL1 (Qi et al., Reference Qi, Qian, Bu, Li, Chen, Sun, Liang, Zhou, Chu and Li2008), OsNAL11 (Wu et al., Reference Wu, Luo, Chen, Tao, Huang, Wang, Chen and Xiao2016; Zhao et al., Reference Zhao, Luo, Jiang, Yang and Zha2017), OsNAL2 (Cho et al., Reference Cho, Yoo, Zhang, Pandeya, Koh, Hwang, Kim and Paek2013), OsNAL3 (Cho et al., Reference Cho, Yoo, Zhang, Pandeya, Koh, Hwang, Kim and Paek2013), OsNAL7 (Fujino et al., Reference Fujino, Matsuda, Ozawa, Nishimura, Koshiba, Fraaije and Sekiguchi2008), OsNAL9 (Li et al., Reference Li, Wu, Hu, Xing, Qian, Si, Sun, Wang, Fu and Liu2013b), OsNRL1 (Hu et al., Reference Hu, Zhu, Zeng, Gao, Guo, Fang, Zhang, Dong, Yan and Liu2010; Wu et al., Reference Wu, Fu, Hu, Si, Cheng and Liu2010), OsNRL4 (Liang et al., Reference Liang, Qin, Zeng, Zheng, Jin and Shi2016), OsREL1 (Chen et al., Reference Chen, Xie, Gao, Wang, Sun, Liu, Zhu, Peng, Zhao and Liu2015), OsREL2 (Yang et al., Reference Yang, Li, Miao, Gan, Qiao, Chang, Shi and Chen2016), OsRFS (Cho et al., Reference Cho, Lee, Gi, Yim, Koh, Kang and Paek2018; Yoshimura et al., Reference Yoshimura, Ideta and Iwata1997), OsRL14 (Fang et al., Reference Fang, Zhao, Cong, Sang, Du, Wang, Li, Ling, Yang and He2012), OsRL15 (Lee et al., Reference Lee, Woo, Lee, Lee, Kim and Koh2016), OsRL16 (Liu et al., Reference Liu, Kong, You, Zhong, Jiang, Zhao and Wan2015), OsRL9 (Yan et al., Reference Yan, Yan, Zeng, Yang, Fang, Tian, Sun, Cheng and Gu2008), OsRRK1 (Ma et al., Reference Ma, Zhao, Shangguan, Shi, Zeng, Wu, Chen, You, Zhu, Du and He2017), OsRoc5 (Zou et al., Reference Zou, Sun, Zhang, Liu, Wu, Tian, Qiu and Lu2011), OsSCL1 (Shen, Reference Shen2011), OsSFL1 (Alamin et al., Reference Alamin, Zeng, Qin, Sultana, Jin and Shi2017), OsSLL2 (Zhang et al., Reference Zhang, Wu, Jiang, Wang, Zhang, Guo, Wu and Wan2015), OsSND2 (Ye et al., Reference Ye, Wu, Chen, Liu, Wu, Liu and Fu2018), OsSRL1 (Xiang et al., Reference Xiang, Zhang, Qian and Xue2012), OsSRL2 (Liu et al., Reference Liu, Li, Liu, Tang, Sun, Li, Shen, Du and Cheng2016), OsSRS5 (Segami et al., Reference Segami, Kono, Ando, Yano, Kitano, Miura and Iwasaki2012), OsYABBY1 (Dai et al., Reference Dai, Zhao, Ma, Hu, Hedden, Zhang and Zhou2007), OsYABBY6 (Xia et al., Reference Xia, Tang, Yang, Li, Wang, Lu, Liu and Lin2017) and OsZHD1 (Xu et al., Reference Xu, Wang, Long, Huang, Wang, Zhou, Zheng, Sun, Chen and Chen2014). Then we downloaded the genomic sequence, coding sequence (CDS) and protein sequence from the Rice Genome Annotation Project database (http://rice.plantbiology.msu.edu/index.shtml) for each of the 42 RL genes included in our study. For convenience, we defined these 42 finally selected RL genes as ‘rolling leaf genes of interest’ in our study. The collected RL dataset that is used in this research is deposited at http://www.bbcba.org/softwares/Rolling_Leaf_Dataset.zip.
Bioinformatics analysis of RL genes
We analysed the genomic sequences, CDS sequences and protein sequences of the RL genes using various software packages and web-based tools. The gene structure analysis was done based on genomic sequences and CDS using the web-based bioinformatics tool Gene Structure Display Server (GSDS2.0: http://gsds.cbi.pku.edu.cn) (Hu et al., Reference Hu, Jin, Guo, Zhang, Luo and Gao2015). In this analysis, we investigated the exon-intron structure along with the number of exons and introns of 42 RL genes. The NCBI ‘Batch CD-Search’ web-based tool (Marchler-Bauer and Bryant, Reference Marchler-Bauer and Bryant2004; Marchler-Bauer et al., Reference Marchler-Bauer, Lu, Anderson, Chitsaz, Derbyshire, DeWeese-Scott, Fong, Geer, Geer, Gonzales, Gwadz, Hurwitz, Jackson, Ke, Lanczycki, Lu, Marchler, Mullokandov, Omelchenko, Robertson, Song, Thanki, Yamashita, Zhang, Zhang, Zheng and Bryant2011; Marchler-Bauer et al., Reference Marchler-Bauer, Derbyshire, Gonzales, Lu, Chitsaz, Geer, Geer, He, Gwadz, Hurwitz, Lanczycki, Lu, Marchler, Song, Thanki, Wang, Yamashita, Zhang, Zheng and Bryant2015) in Conserved Domain Database (CDD: http://www.ncbi.nlm.nih.gov/cdd) was used to perform conserved domain analysis based on the protein sequences of the identified RL genes. In our study, MEGA V6 (Tamura et al., Reference Tamura, Stecher, Peterson, Filipski and Kumar2013) software was used for the multiple sequence alignment and phylogenetic analysis. Multiple alignments of protein sequences of RL genes was conducted using the ‘ClustalW’ method (Thompson et al., Reference Thompson, Higgins and Gibson1994) and was viewed in a compatible printable format using GeneDoc Version 2.6.002 software. Similar to Rahman et al. (Reference Rahman, Yang, Xu, Munyampundu and Cai2016), the maximum likelihood method was used to construct the phylogenetic tree based on the multiple aligned protein sequences.
GO, TF and KEGG pathway analysis
The web-based tool AgriGO v2.0 (http://systemsbiology.cau.edu.cn/agriGOv2/) was used for GO analysis to investigate the functional enrichments of the RL genes of interest. Oryza sativa japonica was selected for species and Singular Enrichment Analysis was selected as an analysis tool with different settings of reference backgrounds. A GO term was considered significantly enriched among a set of genes if P-value (P) is <0.05. We constructed the heatmap of GO term versus Gene along with the bar plot of GO term versus −log10(P) for all significant GO terms. TFs were identified by using the PlantTFDB 4.0 database (http://planttfdb.cbi.pku.edu.cn/) (Jin et al., Reference Jin, Tian, Yang, Meng, Kong, Luo and Gao2017). KEGG Automated Annotation Server (KAAS: http://www.genome.jp/tools/kaas/) (Yoshizawa et al., Reference Yoshizawa, Itoh, Okuda, Moriya and Kanehisa2007) was used for the KEGG pathway enrichment analysis to get the summary of the gene pathway network.
Protein network association construction
Functional protein association networks analysis was done using the web-based resource STRING V10.5 (http://string-db.org/) (Szklarczyk et al., Reference Szklarczyk, Morris, Cook, Kuhn, Wyder, Simonovic, Santos, Doncheva, Roth, Bork, Jensen and Mering2017) to examine the protein–protein interactions among the RL genes. The prediction of interactions was done with the default settings of the STRING database.
Results
Summary of genomic information of identified RL genes
The comprehensive information of 42 RL genes of rice including locus ID, gene location, genomic sequence length, CDS length, number of introns and exons, predicted protein length/size, molecular weight (Mol. Wt., kDa) and isoelectric point (pI), is shown in Table 1. Among the 42 RL genes, 3, 6, 9, 3, 1, 3, 5, 1, 3, 3, 2 and 3 genes were located on chromosomes 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12, respectively. The length of RL genes varied from 315 to 16638 base pairs (or nucleotides). The RL gene OsREL1 has the minimum genomic length (315 nucleotides) and OsADL1 has the maximum genomic length (16638 nucleotides). The RL gene OsSLL2 has the lowest CDS length (225 nucleotides) and protein length (74 amino acids) while the gene OsRFS has the highest CDS length (6579 nucleotides) and protein length (2192 amino acids). Gene OsRFS has the highest molecular weight (243660 kDa) and OsSLL2 has the lowest molecular weight (8408.14 kDa). The gene OsREL1 corresponds to the minimum value of pI (4.2637) whereas the gene OsNRL4 exhibits the maximum value of pI (11.8634).
Gene name, Gene location, Gene length (number of nucleotides), CDS length (number of nucleotides), No. of introns and exons, protein length (number of amino acids), Molecular Weight (Mol.Wt.) in kDa, isoelectric point (pI) of rolling leaf genes
Gene structure analysis of RL genes
The gene structure and exon-intron pattern of 42 RL genes of interest were comparatively analysed and presented in online Supplementary Fig. S1. The number of introns in the RL genes of interest ranges from 0 to 29. The gene OsADL1 has the maximum number of exons (30) whereas the genes OsACL1, OsI_14279, OsNAL3, OsNRL4, OsREL1, OsSCL1, OsSLL2 and OsZHD1 have the lowest number of exon(s) (only one exon). OsAGO1a and OsSRL2 have the second and third highest number of exons (23 and 19, respectively). Each of the RL genes OsDCL1, OsHB4 and OsRL16 gene has 18 exons whereas OsRFS has 11 exons; and the OsLBD3-7 gene has 10 exons. Furthermore, every gene in the group (OsNAL9, OsRL15 and OsRoc5); the gene OsRL9; and each gene in the gene set (OsRRK1, OsYABBY1 and OsYABBY6) have 9, 7 and 6 exons, respectively, in our study. Moreover, each candidate of OsNAL1 and OsSRL1 genes; every single gene in the set (OsARVL4, OsLC2, OsLRRK1, OsNAL7, OsREL2 and OsSRS5); each of the genes in the group (OsAGO7, OsARF18, OsMYB103L and OsSND2); and each candidate of the gene cluster (OsAS2, OsCFL1, OsNAL11, OsNAL2, OsNRL1, OsRL14 and OsSFL1) have 5, 4, 3, 2 exons, respectively. Although each of the genes OsAS2, OsCFL1, OsNAL11, OsNAL2, OsNRL1, OsRL14 and OsSFL1 has two exons and one intron, the length of the exons and introns are different from each other indicating diversity in genes structure. The genes OsNAL1, OsSFL1, OsRL9, OsSRL1, OsYABBY1 and OsYABBY6 contain an intron with an unusually large size (online Supplementary Fig. S1). Genes OsADL1, OsNAL1, OsNAL7, OsSLL2 and OsSRL2 have long upstream; OsAS2, OsI_14279, OsDCL1, and OsSRL1 longer downstream; OsCFL1 and OsDCL1 do not have upstream; and OsNAL2, OsNAL3, OsNRL4, OsREL1, OsRFS and OsRL9 do not have upstream and downstream. The gene structure analysis shows that the genes OsADL1, OsDCL1, OsHB4, OsRRK1 and OsSRL2 are conserved in the sense that they have short length introns (online Supplementary Fig. S1).
Domain analysis of rice RL proteins
Online Supplementary Table S2 and Fig. S2(b) display the results of the domain analysis of the 42 RL proteins identified in this study using the CDD web-tool of NCBI. Conserved domain analysis shows that the proteins OsACL1, OsNRL4, OsREL1, OsSLL2, OsSRL1, OsSRL2 and OsNAL1 do not contain any database-recognized domain in spite of having an important role in leaf rolling of rice. The RL proteins OsARVL4, OsAS2, OsCFL1, OsLRRK1, OsNAL11, OsNAL2, OsNAL3/OsWOX3, OsNAL9, OsRL14, OsRL15, OsRL16, OsRL9/OsSLL1, OsRRK1, OsSCL1, OsSFL1, OsSND2, OsSRS5, OsYABBY1 and OsYABBY6 contain one conserved domain; proteins OsAGO1a, OsARF18, OsI_14279, OsLC2, OsREL2 and OsZHD1 contain two conserved domains; OsADL1, OsLBD3-7, OsNAL7/OsCOW1, OsNRL1 and OsRoc5 proteins contain three conserved domains; OsHB4 and OsMYB103L contain four conserved domains; protein OsAGO7 contains five conserved domains; OsDCL1 contains six conserved domains; and OsRFS contains seven conserved domains. So, the RL protein OsRFS contains the maximum number of domains (seven domains), and the protein OsARVL4, OsAS2, OsCFL1, OsLRRK1, OsNAL11, OsNAL2, OsNAL3/OsWOX3, OsNAL9, OsRL14, OsRL15, OsRL16, OsRL9/OsSLL1, OsRRK1, OsSCL1, OsSFL1, OsSND2, OsSRS5, OsYABBY1 and OsYABBY6 contain the minimum number of domains (one domain). The OsAGO1a and OsAGO7 proteins encode a common conserved domain PIWI. The RL proteins OsRL9/OsSLL1 and OsMYB103L contain one similar domain ‘Myb-like DNA-binding domain’. The OsAGO7 and OsDCL1 RL proteins contain a PAZ domain. A single similar Homeobox domain was encoded by each of the proteins OsHB4, OsNAL2, OsNAL3/OsCOW1, OsRoc5 and OsZHD1. The RL protein OsMYB103L contains two SANT domains and two Myb_DNA_binding domains. OsRFS contains two CHROMO domains. Proteins OsRoc5 and OsHB4 encodes one Basic leucine zipper (bZIP) domain. One common domain YABBY belongs to both the proteins OsYABBY1 and OsYABBY6. The RL proteins OsADL1, OsARF18, OsARVL4, OsAS2, OsCFL1, OsI_14279, OsLBD3-7, OsLC2, OsLRRK1, OsNAL11, OsNAL7, OsNAL9, OsNRL1, OsREL2, OsRL14, OsRL15, OsRL16, OsRRK1, OsSCL1, OsSFL1, OsSND2 and OsSRS5 do not have any common domain.
Phylogenetic analysis of RL related genes
To investigate the evolution and diversification of RL genes, we constructed a maximum likelihood based phylogenetic tree on the alignment of full length RL proteins (Fig. 2). The phylogenetic tree showed that the 42 RL proteins clustered into five major groups (I–V) with the groups III and IV further divided into two subgroups while group V separated into four subgroups with robust bootstrap support. The RL proteins OsAS2, OsNAL1, OsSND2, OsRL9, OsCFL1 and OsNRL4 gathered into the group I, while the RL proteins OsRRK1, OsLRRK1, OsSRS5 and OsACL1 formed the group II. RL proteins OsNRL1, OsLBD3-7, OsRL14, OsNAL11 and OsRL16 constructed the subgroup IIIa while the proteins OsRoc5 and OsHB4 congregate in subgroup IIIb. Proteins OsRFS, OsRL15, OsNAL7, OsZHD1, OsNAL9, OsADL1 and OsSRL2 clustered into subgroup IVa whereas OsSFL1 and OsSRL1 grouped into subgroup IVb. The RL proteins OsREL2, OsARVL4, OsI_14279, OsDCL1 and OsSLL2; OsMYB103L and OsSCL1; OsAGO1a, OsAGO7, OsYABBY1 and OsYABBY6; and OsNAL2, OsNAL3, OsREL1, OsLC2 and OsARF18 appeared in subgroup Va, Vb, Vc and Vd, respectively. These results can be helpful for further evolutionary and functional studies of RL related genes.
Phylogenetic tree of 42 rolling leaf genes of interest identified in this study. The tree was constructed based on multiple aligned sequences by maximum likelihood method with bootstrap of 1000 in MEGA6. Multiple sequence alignment was performed using ClustalW program in MEGA6. The coloured shapes indicate different clusters of RL proteins. The roman numerals I-Vd indicate groups and subgroups of RL proteins.
GO, TF and KEGG analysis of RL genes
In order to better understand the functional involvement and to characterize the selected RL genes, different types of gene enrichment analyses were conducted. The GO analysis was performed to explain the biological importance of RL genes. The ‘GO analysis toolkit’ of the agriGO v2.0 database was used to perform the GO analysis. Fig. 3(a) represents the heatmap of ‘GO term versus gene’ along with the bar plot of –log10(P) against all significant GO terms (P < 0.05). All of the significant GO terms enriched among the RL genes and their description, the number in the input list and P-value are shown in online Supplementary Table S3. From the analysis result, it was observed that 38 GO terms were significantly enriched for the biological process whereas very few significant overrepresentations were found for cellular component and molecular function (Fig. 3(a) and online Supplementary Table S3). We found that the 42 RL genes clustered into three groups (Group I, Group II and Group III) based on the enriched GO terms. The genes in three groups were enriched in various biological processes. The genes in group I were mostly involved in different types of regulations including regulation of transcription (GO:0045449), regulation of nitrogen compound metabolic process (GO:0051171), regulation of cellular biosynthetic process (GO:0031326), regulation of biosynthetic process (GO:0009889), regulation of primary metabolic process (GO: 0080090) and so on. The genes of this group were also involved in response to the abiotic stimulus (GO:0009628) and various biosynthetic processes (GO:0044249 – cellular biosynthetic process, GO:0034645 – cellular macromolecule biosynthetic process, GO:0009059 – macromolecule biosynthetic process). This observation indicated that the genes in group I are more effective and biologically more functional. The genes in group II were mainly involved in the biological process (GO:0008150) whereas genes in group III were mostly involved in multicellular organismal development (GO:0007275), cell differentiation (GO:0030154), primary metabolic process (GO:0044238) and macromolecule metabolic process (GO:0043170). The biological process ‘multicellular organismal development’ (GO:0007275; P < 4.92 × 10−10) was the most representative process that enriched for a large number of genes belonging to group I and group III.
The enriched GO terms for all RL genes and KEGG analysis results. (a) The GO terms indicated by BP are involved in the biological process, the GO terms indicated by MF are involved in molecular function and the GO terms indicated by CC are involved in the cellular component. (b) Bar chart for all significant KEGG pathways of all RL genes where X-axis represents significant KEGG pathway terms and Y-axis represents the number of genes. (c) Pie chart for all allocated KEGG pathways of all RL genes where the significant pathway terms are grouped and indicated by different colour.
TF analysis showed that a total of 13 genes were related to 10 different TF families (Table 2). Gene sets (OsRoc5, OsHB4), (OsNAL3, OsNAL2) and (OsYABBY1, OsYABBY6) belonged to the TF families HD-ZIP, WOX and YABBY, respectively. Genes OsARF18, OsRL9, OsSCL1, OsAS2, OsMYB103L, OsSND2 and OsZHD1 were the members of the TF families ARF, G2-like, GRAS, LBD, MYB, NAC and ZF-HD, respectively. The KEGG pathway analysis was carried out using the KEGG database to know the biological process in the RL genes of interest (Fig. 3(b) and (c)). Results showed that only 14 out of 42 RL genes were involved in the KEGG pathways (Fig. 3(b)). Of those 14 genes, 50% genes were involved in the metabolism pathway, followed by cellular processes (36%) and 7% genes were involved in both genetic information processing and organismal system in this study (Fig. 3(c)). KEGG pathway enrichment analysis (Fig. 3(b)) also showed that three RL genes functioned in the metabolic pathways; two genes functioned in the cell growth and death, cellular community-eukaryotes pathways; and one gene functioned in the pathways: biosynthesis of secondary metabolites, carbohydrate metabolism, lipid metabolism, amino acid metabolism, translation, transport and catabolism, and ageing. Online Supplementary Table S4 represented the KEGG orthologous of the RL genes of interest along with their description.
Summary of transcription factors (TFs) identified in this study
Network analysis of RL related proteins
Network analysis (Fig. 4) showed that there is no protein–protein interaction for the 23 candidate proteins OsREL2, OsNAL2, OsYABBY1, OsLRRK1, OsACL1, OsARVL4, OsLC2, OsREL1, OsSLL2, OsARF18, OsRL16, OsNRL4, OsNAL3, OsRL15, OsZHD1, OsSRL1, OsNAL1, OsI_14279, OsRFS, OsRRK1, OsRL14, OsNAL11 and OsSFL1. The RL protein OsRL9 is associated with the proteins OsRoc5, OsADL1, OsNAL7, OsAS2, OsNRL1, OsYABBY6, OsHB4 and OsAGO7. The RL proteins within each of the pairs (OsMYB103L and OsSND2), (OsSRS5 and OsNAL9), (OsRoc5 and OsCFL1), (OsRoc5 and OsADL1), (OsADL1 and OsAGO7), (OsLBD3-7 and OsAGO1a), (OsSCL1 and OsDCL1) and (OsSRL2 and OsYABBY6) are associated with each other. Proteins OsAGO1a, OsAGO7 and OsDCL1 are associated with each other. Network analysis predicted the association of OsRL9 with OsRoc5, OsADL1, OsNAL7, OsAS2, OsNRL1, OsYABBY6, OsHB4 and OsAGO7; and association of OsRoc5 with OsADL1 by text mining. The interaction between OsRoc5 and OsCFL1 was experimentally determined as well as predicted by text mining. The interaction between the RL proteins within each of the pairs (OsMYB103L and OsSND2) and (OsSRS5 and OsNAL9) was predicted by text mining and co-expression analysis. Also, the association between OsLBD3-7 and OsAGO1a was experimentally determined as well as predicted by co-expression analysis. The RL protein OsAGO1a and OsAGO7 exhibited association which was known from curated databases as well as predicted by gene text mining and protein homology. The association of OsAGO1a with OsDCL1 and OsAGO7 with OsDCL1 was experimentally determined and known from curated databases, and also predicted by text mining and co-expression analysis.
Protein–protein interaction network of 42 rolling leaf genes of interest. The different methods of prediction of interactions have been indicated with the different coloured connective lines. Protein–protein interaction network was built using the web-based tool STRING V10.5 (http://string-db.org/).
Discussion
The ultimate goal of rice breeding is the development of rice mutants with super-high yield and stress tolerance traits. In rice, appropriate leaf rolling is considered to be an important agronomic element, particularly, a moderate leaf rolling is regarded as a crucial phenotypic trait of the ideal rice plant as it is important for increasing grain yield (Yuan, Reference Yuan1997; Lang et al., Reference Lang, Zhang, Gu, Yang and Zhu2004; Wu, Reference Wu2009; Zhang et al., Reference Zhang, Xu, Zhu, Qian and Xue2009). However, to the best of our knowledge, there is no study where all the RL genes/QTLs and their genomic information have been put together along with their genome-wide comparative analysis from the different bioinformatics point of view. In this study, we listed up 103 RL genes/QTLs along with their genomic information reported in various studies till date and performed various comparative analyses from different bioinformatics standpoints with the RL genes of interest.
We selected 42 genes as the genes of interest in our analysis among the total 103 identified RL genes/QTLs due to the availability of their locus IDs in this study (Table 1). The maximum number of RL genes was found on chromosome 3 and minimum genes were found on chromosomes 5 and 8 among the 42 RL genes. The RL gene OsADL1 has the longest genomic length and OsREL1 has the shortest genomic length. OsSLL2 has the lowest value of CDS and protein length while the gene OsRFS has the highest value of CDS and protein length. OsRFS has the highest molecular weight and OsSLL2 has the lowest molecular weight. The gene OsREL1 corresponds to the minimum value of pI whereas the gene OsNRL4 exhibits the maximum value of pI. The comparative gene structure analysis revealed that some RL genes contained only one exon, some genes did not have upstream and some genes did not contain any upstream and downstream (online Supplementary Fig. S1). Most of the genes (60%) contained 0 to 3 introns. Although some genes contained an equal number of exons and introns the length of exons and introns were different from each other indicating structural diversity. From the gene structure analysis (online Supplementary Fig. S1) it can be inferred that RL genes are diverse in their structures from each other in terms of conserveness and structure of exons and introns.
Domain analysis revealed that very few of the RL proteins contain the domains of the same family. Both the RL genes proteins OsNAL2 and OsNAL3 contained only one and identical domain named ‘Homeobox domain’ (online Supplementary Table S2 and Fig. S2(b)). According to Cho et al. (Reference Cho, Yoo, Zhang, Pandeya, Koh, Hwang, Kim and Paek2013), the proteins NAL2 and NAL3 were paralogs and encoded an identical homeobox 3A (OsWOX3A) protein which supported our findings. Both the proteins OsYABBY1 and OsYABBY6 contained only one domain (YABBY domain) which is identical in both proteins (online Supplementary Table S2 and Fig. S2(b)). As reported by Toriba et al. (Reference Toriba, Harada, Takamura, Nakamura, Ichikawa, Suzaki and Hirano2007), both the proteins OsYABBY1 and OsYABBY6 belonged to the YABBY family which is in favour of our findings. Most of the RL proteins contained the domains of various families indicating the diversity of RL proteins. This study showed that RL proteins were diverse in terms of domains they contained and their domain structure. Based on the phylogenetic analysis of 42 RL protein sequences, we grouped the RL proteins into five major groups. Again groups III, IV and V were divided into two, two and four subgroups, respectively. We found that both the RL proteins OsNAL2 and OsNAL3 belonged to same subgroup Vd, and the proteins OsYABBY1 and OsYABBY6 belonged to the same subgroup Vc. These findings are supported by Cho et al. (Reference Cho, Yoo, Zhang, Pandeya, Koh, Hwang, Kim and Paek2013) and Toriba et al. (Reference Toriba, Harada, Takamura, Nakamura, Ichikawa, Suzaki and Hirano2007).
To explore functional similarity and characterize a set of genes various enrichment analyses, such as GO and KEGG analyses, were performed. The biological process ‘multicellular organismal development’ (GO:0007275; P < 4.92E-10) was found as the most representative process that enriched for a large number of genes that belonged to groups I and III. The GO term GO:0007275 (multicellular organismal development) was found as the most significant. From the heatmap of the GO terms against gene, three distinct clusters of genes were observed and the genes that belonged to different clusters were involved in different biological processes. We found that OsNAL2 and OsNAL3 are involved in the regulation of the biosynthetic process, which is similar to the findings by Cho et al. (Reference Cho, Kang, Lee, Lee and Paek2016). TFs have significant functions for controlling plant propagation, maturation and react to unfavourable situations including drought, chill, salinity and high temperature (Zhang et al., Reference Zhang, Rerksiri, Liu, Zhou, Xiong, Xiang, Chen and Xiong2013). For example, myb-protein encoding genes have been reported to function as regulators of cell differentiation (Kirik and Bäumlein, Reference Kirik and Bäumlein1996) and NAC TFs involved in various mechanisms include developmental process, flower creation, seed growth, hormone signalling, responses to abiotic stress and ageing in plants (Nuruzzaman et al., Reference Nuruzzaman, Manimekalai, Sharoni, Satoh, Kondoh, Ooka and Kikuchi2010). TF analysis showed that 10 different TFs families (such as HD-ZIP, WOX, ARF, G2-like, GRAS, LBD, MYB, NAC and ZF-HD) involved in leaf rolling, where all of these TFs were responsible for the different aspects of plant developments (Table 2). Metabolomics plays a significant role in fundamental plant biology and applied biotechnology. Several current studies have found that the analysis of the metabolite composition of mutants can assist the assignment of functions to genes (Schauer and Fernie, Reference Schauer and Fernie2006). Results of KEGG pathway analysis showed that half of the RL genes were involved in the metabolism pathway including biosynthesis of secondary metabolites, carbohydrate metabolism, lipid metabolism, amino acid metabolism and metabolic pathways (Fig. 3(b) and (c)). According to Li et al. (Reference Li, Zhang, Gan, Qiao, Yang, Miao, Wang, Zhang, Liu, Li, Shi and Chen2017), the proteins responsible for leaf rolling in rice are involved in different pathways including biosynthesis of phenylpropanoid, phenylalanine metabolism, and sucrose and starch metabolism. As reported by Chen et al. (Reference Chen, Wan, Shen, Zhou, Huang, Chu and Guan2018), RL proteins in Brassica napus are involved in different pathways including ‘metabolic pathways’, ‘biosynthesis of secondary metabolites’ and ‘starch and sucrose metabolism’. So, the results of previous studies of RL support our finding.
The network analysis was used for investigating the association among the 42 RL proteins. Results showed that the RL protein OsRL9 was associated with a maximum number of proteins (9 proteins). We found that the proteins OsAGO1a, OsAGO7 and OsDCL1 were associated with each other. Protein homology and text mining predicted the association between the proteins OsAGO1a and OsAGO7 which was also known from the curated databases. All the four methods (experimentally determination, determination from curated databases, prediction using text mining and prediction by co-expression analysis) identified the association between the genes of each of the pairs of RL proteins (OsAGO1a and OsDCL1) and (OsAGO7 and OsDCL1). The protein OsRoc5 was associated with OsCFL1 and OsADL1; OsMYB103L was associated with OsSND2; and OsLBD3-7 was associated with OsAGO1a. Chen et al. (Reference Chen, Sheng, Cai, Li, Wei, Xie, Jiao, Shao, Tang, Wang and Hu2019) reported direct interaction between the proteins OsNAL9 and OsNRL3 whereas due to the unavailability locus ID for the gene OsNRL3 investigation of the interaction between the proteins OsNAL9 and OsNRL3 was not possible. However, in this study, we found the interaction of OsNAL9 with OsSRS5. We found 23 RL proteins (OsREL2, OsNAL2, OsYABBY1, OsLRRK1, OsACL1, OsARVL4, OsLC2, OsREL1, OsSLL2, OsARF18, OsRL16, OsNRL4, OsNAL3, OsRL15, OsZHD1, OsSRL1, OsNAL1, OsI_14279, OsRFS, OsRRK1, OsRL14, OsNAL11 and OsSFL1) which were not associated with any other RL proteins. This indicates that more than 50% of the RL proteins of interest are not associated with each other.
Altogether, this is the first study where we have studied and put the genetic information of all identified RL genes/QTLs till now and their genomic information along with their genome-wide comparative analysis from the different bioinformatics point of view. This study will enable the rice breeders and researchers to get collective information on RL genes/QTLs along with their comparative analysis results. As a result, the rice breeders will be able to develop a super-high-yield rice mutant with moderate leaf rolling and desired architecture.
Conclusion
Our present study identified 103 RL genes/QTLs in rice genome characterized by several studies to date. Among 103 RL genes/QTLs, 42 genes, for which locus IDs were available, were finally selected in our analysis. This study provides a comparative genomic analysis of the selected 42 RL genes to explore their structures, functions and diversifications using different bioinformatics techniques. Gene structure analysis shows that the selected RL genes are diverse in structures in terms of intron-exon. Domain analysis reveals that the RL proteins contain different types of domains except for some proteins indicating their functional diversification. Phylogenetic analysis has clustered the RL proteins into five major groups with groups III, IV and V divided into some subgroups. GO analysis indicates that a total of 42 significant GO terms enriched, TF analysis shows 10 different TFs families and KEGG analysis shows that 14 genes are involved in the KEGG pathways from all RL genes. Protein–protein interaction network shows that 23 RL proteins are not associated with any other RL proteins and the protein OsRL9 is associated with maximum proteins. Therefore, we may conclude that the RL genes have different types of functions for controlling leaf rolling in rice. These results provide insight into the structures, predicted functions, evolutionary information and diversifications of RL related genes in rice which will be a basis for further functional analysis of selected candidate genes. The application of machine learning methods has been booming in the past several years (Khatun et al., Reference Khatun, Hasan and Kurata2019a, Reference Khatun, Hasan and Kuratab; Hasan et al., Reference Hasan, Basith, Khatun, Lee, Manavalan and Kurata2020a, Reference Hasan, Manavalan, Khatun and Kuratab, Reference Hasan, Manavalan, Shoombuatong, Khatun and Kuratac, Reference Hasan, Manavalan, Shoombuatong, Khatun and Kuratad, Reference Hasan, Schaduangrat, Basith, Lee, Shoombuatong and Manavalane; Mosharaf et al., Reference Mosharaf, Hassan, Ahmed, Khatun, Moni and Mollah2020), prompting researchers to design machine learning architectures to predict diverse plant species. As for future studies on structures, functions and diversifications of plant genes, machine learning methods will be employed to propel the RL related genes prediction to the next level.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S1479262120000404.
Acknowledgements
We would like to thank the learned reviewers whose valuable comments and suggestions helped to strengthen this paper. Also we would like to acknowledge HEQEP sub-project (CP-3603,W2, R3), and Bioinformatics Laboratory, Department of Statistics, University of Rajshahi, for providing access to use all the required resources during this research work.
