Rice (Oryza sativa) is an important food grain crop that is consumed worldwide (Kubo and Purevdorj, 2004). The human population is estimated to reach approximately 10.7 billion by 2050 (Speidel, 2000) and accordingly, the demand for consumption of rice will also increase. On the other hand, the productivity of rice is not increasing at the same pace due to various reasons such as poor water management, soil depletion, abiotic stresses (e.g., drought, flooding, and salinity), biotic stresses (e.g., damage caused by bacteria, fungi, insects), etc. (Kumar et al., 2017). Therefore, it becomes necessary to understand the underlying mechanisms in crops to sustain production and fulfill the demand of the growing population (Ma and Michailides, 2007). Plants must be able to react rapidly with various abiotic and biotic stress signals and develop efficient defense responses to cope with adverse conditions arising in the field (Kumar et al., 2015). Understanding the crosstalk mechanisms in a stress response could help in the characterization of synergistic and antagonistic mechanisms (Kim et al., 2021). Plant hormones are the key components of these defense and adaptation mechanisms. To facilitate the adaptations to the stresses, various hormonal pathways are upregulated or downregulated. Change in hormonal influence usually affects the degree of resistance or susceptibility to different stresses (O’Brien and Benkova, 2013). Cytokinins (a class of phytohormones) are central regulators of plant growth and development. Cytokinins regulate numerous developmental and physiological processes such as cytokinesis of shoots and roots, reproductive behavior, leaf senescence, and responses to environmental cues, particularly to light and nutrients (Haberer and Kieber, 2002; Werner et al., 2003). Trans-Zeatin (tz) is an active form of cytokinin involved in managing environmental stress. Cytokinin pathway has been widely studied and a huge amount of gene expression data are available in public repositories (Edgar et al., 2002; Leinonen et al., 2011). These data can be better utilized for constructing gene regulatory networks and identifying key genes which will further help in developing improved rice varieties having the ability to produce high yield and resistance to such abiotic stress and adverse conditions (Stuart et al., 2003; Lelandais et al., 2011).

Key genes regulating plant growth and cytokinin responsive genes involved in development process will help in developing better stress tolerant varieties (Zhang et al., 2011, 2012; Li et al., 2019). Co-expression analysis is one of the important ways to construct such a network and identify the most relevant module (Heyer et al., 1999; Hudson et al., 2012; Gaiteri et al., 2014). A statistical approach, Weighted Gene Co-expression Network Analysis (WGCNA) is an effective way to perform such analysis (Tang et al., 2016; Che et al., 2018; Wu et al., 2018) and it has been successfully used previously in identifying important modules and key hub genes related to rice in different conditions (Tan et al., 2017; Zhang et al., 2018). This approach provides an analytical framework for the analysis of microarray and transcriptomic data and helps in relating the gene expression to external traits. Based on various co-expression analysis studies, we found WGCNA as a well-proven approach for the identification of functional modules and co-expressed genes from large biological datasets (Kost et al., 2017; Zhang et al., 2018; Mishra et al., 2021). WGCNA is available as an R package and works on the principle of “guilt by association” (Langfelder and Horvath, 2008), that is, a group of genes are more likely to be associated with each other when they have similar expressions (Gillis and Pavlidis, 2012; van Dam et al., 2018). It uses an unsupervised learning method where large-scale data is clustered by building a tree from bottom to top by connecting the two nearest branches or genes. Modules are formed using the hierarchical clustering method and comprised of genes with similar functions. The modules can be further utilized to identify the important key genes. These key or hub genes are centrally connected genes in different modules which are expected to have an important function and play a critical role in maintaining overall harmony within the cell and development of tissue (Sircar and Parekh, 2015).

In this paper, we have performed co-expression analysis on publicly available microarray data retrieved from the NCBI GEO database (Edgar et al., 2002; Barrett et al., 2013) for cytokinin-responsive genes. We have identified novel key genes in each module using sound statistical approaches of co-expression analysis. Furthermore, we have done the biological validation of the novel key genes using well-established and experimentally validated QTLs (quantitative trait locus) of rice.

The productivity of rice is severely affected due to various biotic and abiotic stresses. Therefore, to develop a variety that is resistant to these stresses, there is an urgent need to identify important hub genes governing the whole production process. Rice genome consists of around 58,000 genes (Cao et al., 2012) and conventional approaches can identify several hundred genes related to these stresses. However, by using conventional approaches, it becomes difficult to identify the few centrally important genes that play important roles in cellular functions to cope with these stresses. This problem can be solved by applying systems biology approaches (Arora et al., 2019). In this study, we have performed a comprehensive analysis on existing data retrieved from NCBI to understand the potential genes and mechanisms involved in such processes by first performing PCA to validate the distribution and uniformity of data evenly (Figure 10) and subsequently obtained significant modules associated with the biological functions regulating the growth and development of the plant. Moreover, hub genes in these modules were also identified that play an important role during cytokinin signaling and are crucial in plant growth and development.

Amongst the identified modules, the hub genes identified in these modules were observed to be involved in various processes. For example, in the turquoise module, the top hub gene (Os04g0295400) is located in chr4 and its important function is to encode Jasmonate-induced protein. Though little is known about this function in rice, Jasmonate-induced proteins are already reported for immunity and development in other plants (Wasternack and Hause, 2013; Campos et al., 2014). In-depth characterization of this gene is further required as this family of genes are reported to play important roles such as in the defense systems against lethal disease and bacterial blight (Yamada et al., 2012; Taniguchi et al., 2014) and it may also be involved in stress management as both these stresses (abiotic and biotic stress) are interlinked with each other (Cao et al., 2017). Similarly, in the blue module, we identified the Os07g0424400 hub gene in chr7 that played an important role in cellulose synthase A catalytic subunit 3 [UDP-forming] that governs a major mechanism of cell wall formation (Taylor et al., 2000) and ultimately helps in supporting the plant growth against abiotic stress. In the brown module, the top potential hub gene, Os01g0743600, is located on chr1 with the reported function of ATP-dependent protease La domain containing protein (Koodathingal et al., 2009) and it is one of the key components in providing protection against the harmful effects of unfolded proteins. It is activated by stress conditions in the endoplasmic reticulum (ER) and it supports plant defense as well as response to abiotic stresses (Bao and Howell, 2017).

The top potential key gene in the green module, Os07g0171300, is still not fully characterized, but annotation results suggest it has a key influence in the protein kinase-like domain superfamily which is believed to be a conserved protein domain mainly involved in most signaling and regulatory processes in the eukaryotic cell, and as a switch in controlling biological processes such as metabolism, transcription, cell moment, apoptosis, etc.

In the yellow module, the identified key gene Os12g0566000 codes for boron transporter 1, and is located at LOC4352546 of Chr12 and mapped with QTL identification number AK100510. Boron is essential for maintaining the integrity of plants cell walls. It exhibits an important structural role in shaping the cell by providing mechanical strength via cross-linking of cell wall rhamnogalacturonan II (RG-II) to form a stable three-dimensional pectic network which contributes to the mechanical properties of cell wall structure (Funakawa and Miwa, 2015). It is also reported that boron expression deficiency inhibits plant photosynthetic capacity (Zhao and Oosterhuis, 2002; Kastori et al., 2008) and directly impacts the total yield of the crop.

In the red module, the identified key gene Os09g0442400 codes for protein Gamete expressed 1, and is located at LOC4347181 of Chr9 and mapped with QTL identification number AK106970. Gamete expressed 1 protein is mainly responsible for fertilization. It has a dual function during gametophyte development and early embryogenesis and it is required for correct pollen maturation (Engel et al., 2005). In the pink module, the identified key gene Os07g0187700 codes for protein SEC12-like protein1, and is located at LOC4342601 of Chr7 and mapped with QTL identification number AK111777. Phosphate Transporter 1 (PHT1) is a plant-specific SEC12 gene that encodes phosphate transporter involved in phosphate uptake by facilitating the trafficking of PHT1-1/PHT1; 1 from the ER to the plasma membrane that enables the ER exit of a high-affinity phosphate transporter (González et al., 2005). The top key hub gene identified in the grey module, Os07g0131600, codes for HEX6 protein which is one of the hexose carrier proteins, and is located at LOC4342334 of Chr7 mapped with QTL identification number AK068296. HEX6 protein has an active uptake of hexose with an important role in glucose/hydrogen symport (Boles and Hollenberg, 1997). These different hub genes directly or indirectly govern the main function of the positive build-up of overall cell growth but there are leftover top key genes Os07g0171300 in the green module and Os11g0491400 in the black modules that need to be further characterized to understand their role at different stages of crop development.

Further, the identified key hub genes were visualized using the STRING database for protein-protein interaction (Szklarczyk et al., 2017). The STRING database helped in the identification of direct (physical) interactions and indirect (functional) interactions as long as the interactions were specific and biologically meaningful. Out of the 9 obtained modules, seven genes were found to be associated with neighboring proteins (Figure 11).

Similarly, we identified other important hub genes in each module (Table 1 to Table 9) which are not yet fully explored with respect to cytokinin signaling to maintain the harmony of cell and rice growth mechanisms. Although we identified 36 key genes, we were interested in understanding the role of genes in various processes such as BP, CC, and MF to further delineate the role of these genes in playing different roles in the development of cytokinin-related responses. The Gene Ontology terms in these processes include transcription (GO:0006351), auxin-activated signaling pathway (GO:0009734), MAPK cascade (GO:0000165), and regulation of transcription in BP (GO:0006355) (Figure 12A), whereas the GO terms mostly constituted in CC include nucleus (GO:0005634), cytoplasm (GO:0005737), chloroplast (GO:0009507), and microtubule (GO:0005874) (Figure 12B) and composed of molecular functions metal ion binding (GO:0046872), protein serine/threonine phosphatase activity (GO:0004722), ATP-dependent RNA helicase activity (GO:0004004), and RNA polymerase II regulatory region sequence-specific DNA binding (GO:0000977) are some (Figure 12C). Likewise, the Kegg pathway also revealed the enrichment of GO terms such as plant hormone signal transduction (osa04075), starch and sucrose metabolism (osa00500), and diterpenoid biosynthesis (osa00904) (Figure 12D). In the absence of wet lab experiments, these identified hub genes were validated with the help of well-known QTLs and pathways using an in-silico approach. The analysis indicates the involvement of identified hub genes in these stress conditions as these genes are found to be associated with biotic and abiotic stress-related QTLs and pathways. These results indicate the role of identified hub genes in the regulation of plant growth and development. However, these hub genes need further attention at the molecular level through wet lab experiments to improve the traits which will be useful in enhancing the productivity of the crop.