Skip to main content

ORIGINAL RESEARCH article

Front. Genet., 07 February 2022
Sec. Genomics of Plants and the Phytoecosystem
This article is part of the Research Topic Harnessing Cytokinin Biology in Crop Biofortification and Enhanced Food Security View all 16 articles

Identification of Potential Cytokinin Responsive Key Genes in Rice Treated With Trans-Zeatin Through Systems Biology Approach

  • 1Centre for Agricultural Bioinformatics, ICAR-Indian Agricultural Statistics Research Institute, New Delhi, India
  • 2National Institute of Animal Science, Rural Development Administration, Jeonju, South Korea
  • 3ICAR-National Institute for Plant Biotechnology, New Delhi, India
  • 4Agricultural Education Division, Indian Council of Agricultural Research, New Delhi, India

Rice is an important staple food grain consumed by most of the population around the world. With climate and environmental changes, rice has undergone a tremendous stress state which has impacted crop production and productivity. Plant growth hormones are essential component that controls the overall outcome of the growth and development of the plant. Cytokinin is a hormone that plays an important role in plant immunity and defense systems. Trans-zeatin is an active form of cytokinin that can affect plant growth which is mediated by a multi-step two-component phosphorelay system that has different roles in various developmental stages. Systems biology is an approach for pathway analysis to trans-zeatin treated rice that could provide a deep understanding of different molecules associated with them. In this study, we have used a weighted gene co-expression network analysis method to identify the functional modules and hub genes involved in the cytokinin pathway. We have identified nine functional modules comprising of different hub genes which contribute to the cytokinin signaling route. The biological significance of these identified hub genes has been tested by applying well-proven statistical techniques to establish the association with the experimentally validated QTLs and annotated by the DAVID server. The establishment of key genes in different pathways has been confirmed. These results will be useful to design new stress-resistant cultivars which can provide sustainable yield in stress-specific conditions.

Introduction

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.

Materials and Methods

Data Summary

Microarray data related to cytokinin effect on root and leaves of rice were retrieved from the NCBI GEO database with reference series GSE6719 (Hirose et al., 2007). In this experiment, roots and leaves of rice seedlings were treated with trans-zeatin (an active form of type-A cytokinin hormone) for 30 and 120 min. Three independent replicate treatments were performed for roots and leaves at each time point. To identify cytokinin responsive genes, dimethyl sulfoxide (DMSO) treated roots and leaves for each time point (30 and 120 min) were used as a control. Three independent replicates of DMSO treatment were performed for each organ per time point. Then, the microarray data were generated by using the Affymetrix GeneChip® rice genome array which contains probes to query approximately 48,564 japonica and 1,269 indica transcripts. There are 24 samples having accession series (GSE6719). The number of transcripts (features) in the expression data was 57,381, of which annotation data was available for only 23,850 transcripts. We have used only the annotated transcripts for WGCNA-based co-expression analysis.

Weighted Gene Co-Expression Analysis

To perform the analysis, we first studied the distribution pattern of grown samples using transcript count data and performed principal component analysis. Following, uniform distribution pattern we have used the R package “WGCNA” for constructing the network (Langfelder and Horvath, 2008). The elements in the co-expression matrix are defined as the weighted value of correlation coefficients. Gene co-expression of a pair of genes i and j were calculated using Eq. 1 of an unsigned network:

Aij=|cor(Ei,Ej)|β(1)

where Ei and Ej consist of expression vector profiles of genes i and j across multiple samples and Aij is the adjacency of the unsigned network. Pearson’s correlation coefficient was used to identify the similarity of genes. The absolute value of correlation is raised to a power β to create the adjacency matrix. WGCNA emphasizes high correlation by raising the absolute value of correlation to a power β ≥ 1. The correlation network among the genes is built based on the adjacency matrix. Each gene is considered as a node and an edge between two genes is built if it passes a set threshold of co-expression strength. A soft threshold value of 6 was used to determine the significant edges for connecting pairs of genes. Furthermore, modules within the networks were identified by calling the R function “cutreeDynamic” available in the WGCNA package, which helps in identifying the minimum number of large modules with a strong association of genes. The function was applied after identifying the soft threshold of the adjacency matrix, and cut them respectively to get modules related to different functions. These identified modules were further used to identify the key or hub genes responsible for governing specific biological functions in cell development.

In order to identify the modules that were highly related, module similarity was quantified by eigengenes correlation (Shi et al., 2015). The eigengenes of a module are defined as the eigenvector associated with the first principal component of the expression matrix (Langfelder and Horvath, 2007). Highly related modules were identified using module significance. Module significance is defined as the average gene significance of all the genes in the modules to access the association of a module to the phenotype.

Visualization of Network and Hub Genes Identification

Various methods have been used for the identification of hub genes from a large dataset (Bader and Hogue, 2003; Chin et al., 2014). These methods mainly focused on hub gene identification, based only on gene connection degrees in the gene co-expression network. A heat map was generated to compare the expression pattern across the samples and assigned different colors to genes with similar values (Zhao et al., 2014). Genes are interconnected in each module and possess specific functions. After identifying the most significant module, identification of key or hub genes was carried out using an R package, “dhga” (Differential Hub Gene Analysis) (Das et al., 2017).

Gene Annotation and Gene Ontology Analysis

The annotation of expressed genes was performed using the Institute for Genomic Research (TIGR) database which includes information regarding biological processes (BP), molecular function (MF), and cellular components (CC) (Ouyang et al., 2007). Furthermore, the identified genes from identified modules were submitted to the OryzaExpress database (Hamada et al., 2011) and converted the ids into probeID (Affymatrix). These ids were subsequently converted into geneid using DAVID (the database for annotation, visualization, and integrated discovery) (Huang et al., 2007) for KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis web server to decipher the role of genes varies from BP, CC, MF, and KEGG pathway analysis respectively against oryza sativa database and observed in REVIGO database (Supek et al., 2011). We found the involvement of these hub genes by observing the changes in the expression level of treated and control (normal) conditions [S1_file].

Validation of Identified Hub Genes

We validated the identified hub genes using in-silico approach. In the in-silico approach, we identified the chromosome-wise distribution of the hub genes and mapped these genes on chromosomes using web server oryzabase (https://shigen.nig.ac.jp/rice/oryzabase/) (Kurata and Yamazaki, 2006). Further, we used a statistical framework approach to test the association of these hub genes with well-known experimentally validated and genetically rich QTLs data reported for biotic and abiotic stress conditions (Jaiswal et al., 2002). We have used an R package “GSAQ” (https://cran.r-project.org/web/packages/GSAQ) (Das et al., 2018) for mapping the identified hub genes with QTLs on respective chromosomes. GSAQ provides a platform to associate the selected hub genes to the corresponding overlapped QTL-IDs with their genomic positions. Furthermore, the identified hub genes were also validated through pathway analysis.

Results

Weighted Co-expression Network Construction and Module Identification

Clustering of the samples (Figure 1) suggests that there is no outlier present in the data. Power β was obtained through two types of graphs given in Figure 2: 1) soft threshold values of β (x-axis) vs scale-free topology model fit scaled R2 (y-axis) and 2) soft threshold values of β (x-axis) vs mean connectivity scores (y-axis). The optimal value of β obtained using these graphs is 6 with R2value 0.8. This value of β was further used to produce hierarchical clustering (Figure 3).

FIGURE 1
www.frontiersin.org

FIGURE 1. Clustering dendrogram of samples based on their Euclidean distance.

FIGURE 2
www.frontiersin.org

FIGURE 2. Analysis of network topology for various soft-thresholding powers. The left panel shows the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The right panel displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis).

FIGURE 3
www.frontiersin.org

FIGURE 3. Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned merged module colors and the original module colors.

A dynamic hierarchical tree algorithm was used to divide the clustering tree constructed from the differentially expressed genes, resulting in 24 different co-expression modules in the data which were named as blue (3,972 genes), brown (3,064 genes), green (1801 genes), yellow (2,387 genes), red (1,272 genes), and turquoise (5,670 genes), Black (1,115 genes), Cyan (275 genes), dark green (85 genes), dark red (113), dark turquoise (83 genes), green yellow (338 genes), grey (35 genes), grey60 (162 genes), light cyan (246 genes), light green (144 genes), light yellow (141 genes), magenta (644 genes), midnight blue (252 genes), pink (835 genes), purple (442 genes), royal blue (122 genes), salmon (296 genes), and tan (326 genes) (Figure 4). After clustering, the genes were grouped into modules (subnetworks) depicted in different colors for easy identification (Figure 5).

FIGURE 4
www.frontiersin.org

FIGURE 4. Visualization of the eigengenes network representing the relationships among the modules. The top panel shows a hierarchical clustering dendrogram of the modules based on the dissimilarity score of Ei and Ej given by 1−cor (Ei, Ej). The heatmap in the bottom panel shows the eigengenes adjacency Aij= (1 + cor (Ei, Ej))/2.

FIGURE 5
www.frontiersin.org

FIGURE 5. Visualizing the gene network using a heatmap plot. The heatmap depicts the Topological Overlap Matrix (TOM). The light red color represents low overlap and progressively darker shades of red represent higher overlap. Blocks of dark colors along the diagonal depict modules. The gene dendrogram and module assignment are also shown along the left side and the top of the Figure.

Moreover, in order to explore the relationship between identified modules and the experimental samples (traits), we calculated and tested the correlation coefficients. These correlation coefficients along with their p-values for module-trait relationship were depicted (Figure 6). In Figure 6, the red color shows a strong positive correlation and the blue color displays a strong negative correlation. Furthermore, the centralized hub genes were identified from these modules through statistical analysis with the help of R-package “dhga” using a weighted gene score. A total of 36 hub genes were identified in 9 modules and their detailed description (function/annotation, location, and accession number) are mentioned in Table 1; Table 2; Table 3; Table 4; Table 5; Table 6: Table 7; Table 8; Table 9.

FIGURE 6
www.frontiersin.org

FIGURE 6. Module-trait associations where each row corresponds to a module or eigengene and a column to a trait. Each cell contains the corresponding correlation and p-value. The table is color-coded by correlation according to the color legend mentioned in the Figure.

TABLE 1
www.frontiersin.org

TABLE 1. Identified hub genes with respect to their level of connectivity in the turquoise module.

TABLE 2
www.frontiersin.org

TABLE 2. Identified hub genes with respect to their level of connectivity in the blue module.

TABLE 3
www.frontiersin.org

TABLE 3. Identified hub genes with respect to their level of connectivity in the brown module.

TABLE 4
www.frontiersin.org

TABLE 4. Identified hub genes with respect to their level of connectivity in the green module.

TABLE 5
www.frontiersin.org

TABLE 5. Identified hub genes with respect to their level of connectivity in the yellow module.

TABLE 6
www.frontiersin.org

TABLE 6. Identified hub genes with respect to their level of connectivity in the red module.

TABLE 7
www.frontiersin.org

TABLE 7. Identified hub genes with respect to their level of connectivity in the black module.

TABLE 8
www.frontiersin.org

TABLE 8. Identified hub genes with respect to their level of connectivity in the pink module.

TABLE 9
www.frontiersin.org

TABLE 9. Identified hub genes with respect to their level of connectivity in the grey module.

Validation of Identified Hub Genes

These 36 identified hub genes were mapped on rice chromosomes and 32 genes were found to be located at various rice chromosomes (Figure 7). In the process of finding the association of these genes with the well-known QTLs related to salt, cold, drought, bacterial stress, we found that 17 out of the 36 identified hub genes were associated with these QTLs mapped on various chromosomes (Figure 8). We also performed pathway analysis and found that the expression of these genes either increased or decreased during the period of treatment in five different pathways. Heat map analysis was conducted through DAVID to produce a matrix of enriched GO terms with the identified genes. The green and black color shading on the heat map matrix indicate a positive and negative correlation between the enriched GO term and identified DMSO and tz-treated hub genes, respectively (Figure 9).

FIGURE 7
www.frontiersin.org

FIGURE 7. Mapped hub genes on different chromosomes.

FIGURE 8
www.frontiersin.org

FIGURE 8. Hub gene distribution on chromosome to the corresponding QTL ids.

FIGURE 9
www.frontiersin.org

FIGURE 9. The classification of identified hub genes represented by heat map using DAVID software. It is noted that eight genes (y-axis) are classified into four classes (x-axis).

Discussion

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.

FIGURE 10
www.frontiersin.org

FIGURE 10. Distribution of samples based on different stress at different time intervals and obtained uniform distribution of data in all four coordinates (Szklarczyk et al., 2017).

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).

FIGURE 11
www.frontiersin.org

FIGURE 11. The identified top hub genes are denoted with circles and are known as the first shell of interactors. Each color symbol signifies a specific interaction either known or predicted, as mentioned in previous studies.

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.

FIGURE 12
www.frontiersin.org

FIGURE 12. Gene ontology analysis viz biological process, cellular component, and molecular function terms associated with different modules performed using REVIGO.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

DM, DA, NB, AS, SM, AK, and SS conducted the experiments. DA, NB, and SS analyzed the data, DM and DA drafted the manuscript. MF, SL, AR, PP, and KC proofread the manuscript. AR and KC revised the manuscript. All authors agreed to the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Authors duly acknowledge the support for fund and manpower received from ICAR sponsored CABin scheme network project entitled “Agricultural Bioinformatics and Computational Biology” and also the computational support received from Advanced Supercomputing Hub for OMICS Knowledge in Agriculture (ASHOKA).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.780599/full#supplementary-material

Abbreviations

QTL, Quantitative trait locus; WGCNA, Weighted correlation network analysis; GEO, Gene expression omnibus; tz, trans-zeatin; DMSO, Dimethyl sulfoxide; dhga: Differential hub gene analysis; GSAQ, Gene Set Analysis with QTL; ER, Endoplasmic reticulum; KEGG, Kyoto Encyclopedia of Genes and Genomes; DAVID, Database for annotation, visualization, and integrated discovery; TOM, Topological overlap matrix; GO, Gene ontology; PHT1, Phosphate transporter 1.

References

Arora, D., Chaudhary, R., and Singh, A. (2019). System Biology Approach to Identify Potential Receptor for Targeting Cancer and Biomolecular Interaction Studies of Indole[2,1-a]Isoquinoline Derivative as Anticancerous Drug Candidate against it. Interdiscip. Sci. Comput. Life Sci. 11 (1), 125–134. doi:10.1007/s12539-017-0249-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Bader, G. D., and Hogue, C. W. (2003). An Automated Method for Finding Molecular Complexes in Large Protein Interaction Networks. BMC Bioinformatics 4 (4), 2–27. doi:10.1186/1471-2105-4-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bao, Y., and Howell, S. H. (2017). The Unfolded Protein Response Supports Plant Development and Defense as Well as Responses to Abiotic Stress. Front. Plant Sci. 8, 344. doi:10.3389/FPLS.2017.00344

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2013). NCBI GEO: Archive for Functional Genomics Data Sets-Update. Nucleic Acids Res. 41, D991–D995. doi:10.1093/NAR/GKS1193

PubMed Abstract | CrossRef Full Text | Google Scholar

Boles, E., and Hollenberg, C. P. (1997). The Molecular Genetics of Hexose Transport in Yeasts. FEMS Microbiol. Rev. 21, 85–111. doi:10.1111/J.1574-6976.1997.TB00346.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Campos, M. L., Kang, J.-H., and Howe, G. A. (2014). Jasmonate-triggered Plant Immunity. J. Chem. Ecol. 40, 657–675. doi:10.1007/S10886-014-0468-3

CrossRef Full Text | Google Scholar

Cao, P., Jung, K.-H., Choi, D., Hwang, D., Zhu, J., and Ronald, P. C. (2012). The Rice Oligonucleotide Array Database: an Atlas of rice Gene Expression. Rice 5, 17. doi:10.1186/1939-8433-5-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, X.-Q., Jiang, Z.-H., Yi, Y.-Y., Yang, Y., Ke, L.-P., Pei, Z.-M., et al. (2017). Biotic and Abiotic Stresses Activate Different Ca2+ Permeable Channels in Arabidopsis. Front. Plant Sci. 8, 83. doi:10.3389/FPLS.2017.00083

PubMed Abstract | CrossRef Full Text | Google Scholar

Che, T., Li, D., Jin, L., Fu, Y., Liu, Y., Liu, P., et al. (2018). Long Non-coding RNAs and mRNAs Profiling during Spleen Development in Pig. PLoS One 13, e0193552. doi:10.1371/JOURNAL.PONE.0193552

PubMed Abstract | CrossRef Full Text | Google Scholar

Chin, C.-H., Chen, S.-H., Wu, H.-H., Ho, C.-W., Ko, M.-T., and Lin, C.-Y. (2014). cytoHubba: Identifying Hub Objects and Sub-networks from Complex Interactome. BMC Syst. Biol. 8, 1–7. doi:10.1186/1752-0509-8-S4-S11

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, S., Meher, P. K., Rai, A., Bhar, L. M., and Mandal, B. N. (2017). Statistical Approaches for Gene Selection, Hub Gene Identification and Module Interaction in Gene Co-expression Network Analysis: An Application to Aluminum Stress in Soybean (Glycine max L.). PLoS One 12, e0169605. doi:10.1371/JOURNAL.PONE.0169605

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, S., Rai, A., Mishra, D. C., and Rai, S. N. (2018). Statistical Approach for Gene Set Analysis with Trait Specific Quantitative Trait Loci. Sci. Rep. 8 (8), 1–12. doi:10.1038/s41598-018-19736-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R., Domrachev, M., and Lash, A. E. (2002). Gene Expression Omnibus: NCBI Gene Expression and Hybridization Array Data Repository. Nucleic Acids Res. 30, 207–210. doi:10.1093/NAR/30.1.207

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, M. L., Holmes-Davis, R., and McCormick, S. (2005). Green Sperm. Identification of Male Gamete Promoters in Arabidopsis. Plant Physiol. 138, 2124–2133. doi:10.1104/PP.104.054213

PubMed Abstract | CrossRef Full Text | Google Scholar

Funakawa, H., and Miwa, K. (2015). Synthesis of Borate Cross-Linked Rhamnogalacturonan II. Front. Plant Sci. 6, 223. doi:10.3389/FPLS.2015.00223

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaiteri, C., Ding, Y., French, B., Tseng, G. C., and Sibille, E. (2014). Beyond Modules and Hubs: the Potential of Gene Coexpression Networks for Investigating Molecular Mechanisms of Complex Brain Disorders. Genes, Brain Behav. 13, 13–24. doi:10.1111/GBB.12106

CrossRef Full Text | Google Scholar

Gillis, J., and Pavlidis, P. (2012). “Guilt by Association” Is the Exception rather Than the Rule in Gene Networks. PLOS Comput. Biol. 8, e1002444. doi:10.1371/JOURNAL.PCBI.1002444

PubMed Abstract | CrossRef Full Text | Google Scholar

González, E., Solano, R., Rubio, V., Leyva, A., and Paz-Ares, J. (2005). PHOSPHATE TRANSPORTER TRAFFIC FACILITATOR1 Is a Plant-specific SEC12-Related Protein that Enables the Endoplasmic Reticulum Exit of a High-Affinity Phosphate Transporter in Arabidopsis [W]. Plant Cell 17, 3500–3512. doi:10.1105/TPC.105.036640

PubMed Abstract | CrossRef Full Text | Google Scholar

Haberer, G., and Kieber, J. J. (2002). Cytokinins. New Insights into a Classic Phytohormone. Plant Physiol. 128, 354–362. doi:10.1104/PP.010773

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamada, K., Hongo, K., Suwabe, K., Shimizu, A., Nagayama, T., Abe, R., et al. (2011). OryzaExpress: an Integrated Database of Gene Expression Networks and Omics Annotations in rice. Plant Cel. Physiol. 52 (2), 220–229. doi:10.1093/pcp/pcq195

PubMed Abstract | CrossRef Full Text | Google Scholar

Heyer, L. J., Kruglyak, S., and Yooseph, S. (1999). Exploring Expression Data: Identification and Analysis of Coexpressed Genes. Genome Res. 9, 1106–1115. doi:10.1101/GR.9.11.1106

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirose, N., Makita, N., Kojima, M., Kamada-Nobusada, T., and Sakakibara, H. (2007). Overexpression of a Type-A Response Regulator Alters rice Morphology and Cytokinin Metabolism. Plant Cel Physiol. 48, 523–539. doi:10.1093/PCP/PCM022

CrossRef Full Text | Google Scholar

Huang, D., Sherman, B. T., Tan, Q., Collins, J. R., Alvord, W. G., Roayaei, J., et al. (2007). The DAVID Gene Functional Classification Tool: a Novel Biological Module-Centric Algorithm to Functionally Analyze Large Gene Lists. Genome Biol. 8 (8), R183. doi:10.1186/GB-2007-8-9-R183

PubMed Abstract | CrossRef Full Text | Google Scholar

Hudson, N. J., Dalrymple, B. P., and Reverter, A. (2012). Beyond Differential Expression: the Quest for Causal Mutations and Effector Molecules. BMC Genomics 13 (13), 356. doi:10.1186/1471-2164-13-356

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaiswal, P., Ware, D., Ni, J., Chang, K., Zhao, W., Schmidt, S., et al. (2002). Gramene: Development and Integration of Trait and Gene Ontologies for rice. Comp. Funct. Genomics 3, 132–136. doi:10.1002/CFG.156

PubMed Abstract | CrossRef Full Text | Google Scholar

Kastori, R., Plesničar, M., Panković, D., and Sakač, Z. (2008). Photosynthesis, Chlorophyll Fluorescence and Soluble Carbohydrates in sunflower Leaves as Affected by boron Deficiency. J. Plant Nutr. 18, 1751. doi:10.1080/01904169509365021

CrossRef Full Text | Google Scholar

Kim, H., Kim, H., Seong, P., Arora, D., Shin, D., Park, W., et al. (2021). Transcriptomic Response under Heat Stress in Chickens Revealed the Regulation of Genes and Alteration of Metabolism to Maintain Homeostasis. Animals 11 (8), 2241. doi:10.3390/ani11082241

PubMed Abstract | CrossRef Full Text | Google Scholar

Koodathingal, P., Jaffe, N. E., Kraut, D. A., Prakash, S., Fishbain, S., Herman, C., et al. (2009). ATP-dependent Proteases Differ Substantially in Their Ability to Unfold Globular Proteins. J. Biol. Chem. 284, 18674–18684. doi:10.1074/JBC.M900783200

CrossRef Full Text | Google Scholar

Kost, M. A., Perales, H. R., Wijeratne, S., Wijeratne, A. J., Stockinger, E., and Mercer, K. L. (2017). Differentiated Transcriptional Signatures in the maize Landraces of Chiapas, Mexico. BMC Genomics 18, 1–14. doi:10.1186/S12864-017-4005-Y

PubMed Abstract | CrossRef Full Text | Google Scholar

Kubo, M., and Purevdorj, M. (2004). The Future of rice Production and Consumption. J. Food Distrib. Res. 35, 128–142. doi:10.22004/AG.ECON.27145

CrossRef Full Text | Google Scholar

Kumar, A., Kumar, S., Kumar, A., Sharma, N., Sharma, M., Singh, K. P., et al. (2017). Homology Modeling, Molecular Docking and Molecular Dynamics Based Functional Insights into rice Urease Bound to Urea. Proc. Natl. Acad. Sci. India, Sect. B Biol. Sci. 88, 1539–1548. doi:10.1007/S40011-017-0898-0

CrossRef Full Text | Google Scholar

Kumar, R. R., Goswami, S., Sharma, S. K., Kala, Y. K., Rai, G. K., Mishra, D. C., et al. (2015). Harnessing Next Generation Sequencing in Climate Change: RNA-Seq Analysis of Heat Stress-Responsive Genes in Wheat (Triticum aestivum L.). OMICS: A J. Integr. Biol. 19 (10), 632–647. doi:10.1089/omi.2015.0097

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurata, N., and Yamazaki, Y. (2006). Oryzabase. An Integrated Biological and Genome Information Database for rice. Plant Physiol. 140, 12–17. doi:10.1104/PP.105.063008

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2007). Eigengene Networks for Studying the Relationships between Co-expression Modules. BMC Syst. Biol. 1 (1), 1–17. doi:10.1186/1752-0509-1-54

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinformatics 9 (9), 1. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Leinonen, R., Sugawara, H., and Shumway, M.International Nucleotide Sequence Database Collaboration (2011). The Sequence Read Archive. Nucleic Acids Res. 39, D19–D21. doi:10.1093/NAR/GKQ1019

PubMed Abstract | CrossRef Full Text | Google Scholar

Lelandais, G., Goudot, C., and Devaux, F. (2011). The Evolution of Gene Expression Regulatory Networks in Yeasts. Comptes Rendus Biologies 334, 655–661. doi:10.1016/J.CRVI.2011.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Liu, Z., Mi, M., Zhang, C., Xiao, Y., Liu, X., et al. (2019). Identification of Hub Genes and Key Pathways Associated with Angioimmunoblastic T-Cell Lymphoma Using Weighted Gene Co-expression Network Analysis. Cmar 11, 5209–5220. doi:10.2147/CMAR.S185030

CrossRef Full Text | Google Scholar

Ma, Z., and Michailides, T. J. (2005). Genetic Structure of Botrytis Cinerea Populations from Different Host Plants in California. Plant Dis. 89, 1083–1089. doi:10.1094/PD-89-1083

PubMed Abstract | CrossRef Full Text | Google Scholar

Mishra, D. C., Arora, D., Kumar, R. R., Goswami, S., Varshney, S., Budhlakoti, N., et al. (2021). Weighted Gene Co-expression Analysis for Identification of Key Genes Regulating Heat Stress in Wheat. Cereal Res. Commun. 49 (1), 73–81. doi:10.1007/s42976-020-00072-7

CrossRef Full Text | Google Scholar

O’Brien, J. A., and Benková, E. (2013). Cytokinin Cross-Talking during Biotic and Abiotic Stress Responses. Front. Plant Sci. 4, 451. doi:10.3389/FPLS.2013.00451

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, S., Zhu, W., Hamilton, J., Lin, H., Campbell, M., Childs, K., et al. (2007). The TIGR Rice Genome Annotation Resource: Improvements and New Features. Nucleic Acids Res. 35, D883–D887. doi:10.1093/NAR/GKL976

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, K., Bing, Z. T., Cao, G. Q., Guo, L., Cao, Y. N., Jiang, H. O., et al. (2015). Identify the Signature Genes for Diagnose of Uveal Melanoma by Weight Gene Co-expression Network Analysis. Int. J. Ophthalmol. 8, 269–274. doi:10.3980/J.ISSN.2222-3959.2015.02.10

CrossRef Full Text | Google Scholar

Sircar, S., and Parekh, N. (2015). Functional Characterization of Drought-Responsive Modules and Genes in Oryza Sativa: a Network-Based Approach. Front. Genet. 6, 256. doi:10.3389/FGENE.2015.00256

PubMed Abstract | CrossRef Full Text | Google Scholar

Speidel, J. J. (2000). Environment and Health: 1. Population, Consumption and Human Health. CMAJ 163(5): 551–556. Available at: /labs/pmc/articles/PMC80465/(Accessed September 3, 2021).

PubMed Abstract | Google Scholar

Stuart, J. M., Segal, E., Koller, D., and Kim, S. K. (2003). A Gene-Coexpression Network for Global Discovery of Conserved Genetic Modules. Science 302, 249–255. doi:10.1126/SCIENCE.1087447

PubMed Abstract | CrossRef Full Text | Google Scholar

Supek, F., Bošnjak, M., Škunca, N., and Šmuc, T. (2011). REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLoS One 6 (7), e21800. doi:10.1371/journal.pone.0021800

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING Database in 2017: Quality-Controlled Protein-Protein Association Networks, Made Broadly Accessible. Nucleic Acids Res. 45, D362–D368. doi:10.1093/NAR/GKW937

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, M., Cheng, D., Yang, Y., Zhang, G., Qin, M., Chen, J., et al. (2017). Co-expression Network Analysis of the Transcriptomes of rice Roots Exposed to Various Cadmium Stresses Reveals Universal Cadmium-Responsive Genes. BMC Plant Biol. 17, 1–18. doi:10.1186/S12870-017-1143-Y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, R., Zhu, W., Song, X., Lin, X., Cai, J., Wang, M., et al. (2016). Genome-Wide Identification and Function Analyses of Heat Shock Transcription Factors in Potato. Front. Plant Sci. 7, 490. doi:10.3389/FPLS.2016.00490

PubMed Abstract | CrossRef Full Text | Google Scholar

Taniguchi, S., Miyoshi, S., Tamaoki, D., Yamada, S., Tanaka, K., Uji, Y., et al. (2014). Isolation of Jasmonate-Induced Sesquiterpene Synthase of rice: Product of Which Has an Antifungal Activity against Magnaporthe Oryzae. J. Plant Physiol. 171, 625–632. doi:10.1016/J.JPLPH.2014.01.007

CrossRef Full Text | Google Scholar

Taylor, N. G., Laurie, S., and Turner, S. R. (2000). Multiple Cellulose Synthase Catalytic Subunits Are Required for Cellulose Synthesis in Arabidopsis. The Plant Cell 12, 2529. doi:10.1105/TPC.12.12.252910.2307/3871246

PubMed Abstract | CrossRef Full Text | Google Scholar

van Dam, S., Võsa, U., van der Graaf, A., Franke, L., and de Magalhães, J. P. (2018). Gene Co-expression Analysis for Functional Classification and Gene-Disease Predictions. Brief. Bioinform. 19, bbw139. doi:10.1093/BIB/BBW139

CrossRef Full Text | Google Scholar

Wasternack, C., and Hause, B. (2013). Jasmonates: Biosynthesis, Perception, Signal Transduction and Action in Plant Stress Response, Growth and Development. An Update to the 2007 Review in Annals of Botany. Ann. Bot. 111, 1021–1058. doi:10.1093/AOB/MCT067

PubMed Abstract | CrossRef Full Text | Google Scholar

Werner, T., Motyka, V., Laucou, V., Smets, R., Van Onckelen, H., and Schmülling, T. (2003). Cytokinin-deficient Transgenic Arabidopsis Plants Show Multiple Developmental Alterations Indicating Opposite Functions of Cytokinins in the Regulation of Shoot and Root Meristem Activity. Plant Cell 15, 2532–2550. doi:10.1105/TPC.014928

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, D., Zhou, M., Li, L., Leng, X., Zhang, Z., Wang, N., et al. (2019). Time Series Analysis of Gene Changes and Processes after Burn with Human Gene Expression Profiles. Burns 45, 387–397. doi:10.1016/j.burns.2018.08.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamada, S., Kano, A., Tamaoki, D., Miyamoto, A., Shishido, H., Miyoshi, S., et al. (2012). Involvement of OsJAZ8 in Jasmonate-Induced Resistance to Bacterial Blight in rice. Plant Cel Physiol 53, 2060–2072. doi:10.1093/PCP/PCS145

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Yang, Y., Wang, Y., Zhang, J., Wang, Z., Yin, M., et al. (2011). Identification of Hub Genes Related to the Recovery Phase of Irradiation Injury by Microarray and Integrated Gene Network Analysis. PLoS One 6, e24680. doi:10.1371/JOURNAL.PONE.0024680

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Zhao, W., Fu, R., Fu, C., Wang, L., Liu, H., et al. (2018). Comparison of Gene Co-networks Reveals the Molecular Mechanisms of the rice (Oryza Sativa L.) Response to Rhizoctonia solani AG1 IA Infection. Funct. Integr. Genomics 18, 545–557. doi:10.1007/S10142-018-0607-Y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Yu, S., Zuo, K., Luo, L., and Tang, K. (2012). Identification of Gene Modules Associated with Drought Response in rice by Network-Based Analysis. PLoS One 7, e33748. doi:10.1371/JOURNAL.PONE.0033748

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, D., and Oosterhuis, D. M. (2002). Cotton Carbon Exchange, Nonstructural Carbohydrates, and boron Distribution in Tissues during Development of boron Deficiency. Field Crops Res. 78, 75–87. doi:10.1016/S0378-4290(02)00095-3

CrossRef Full Text | Google Scholar

Zhao, S., Guo, Y., Sheng, Q., and Shyr, Y. (2014). Advanced Heat Map and Clustering Analysis Using Heatmap3. Biomed. Res. Int. 2014, 1–6. doi:10.1155/2014/986048

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: WGCNA, systems biology, co-expression, cytokinin, hub genes, QTLs

Citation: Mishra DC, Arora D, Budhlakoti N, Solanke AU, Mithra SVAC, Kumar A, Pandey PS, Srivastava S, Kumar S, Farooqi MS, Lal SB, Rai A and Chaturvedi KK (2022) Identification of Potential Cytokinin Responsive Key Genes in Rice Treated With Trans-Zeatin Through Systems Biology Approach. Front. Genet. 12:780599. doi: 10.3389/fgene.2021.780599

Received: 21 September 2021; Accepted: 18 November 2021;
Published: 07 February 2022.

Edited by:

Jitender Singh, Sardar Vallabhbhai Patel University of Agriculture and Technology, India

Reviewed by:

Bahman Panahi, Agricultural Biotechnology Research Institute of Iran, Iran
Vinay Kumar, National Institute of Biotic Stress Management, India
Ranjeet Ranjan kumar, Indian Agricultural Research Institute (ICAR), India

Copyright © 2022 Mishra, Arora, Budhlakoti, Solanke, Mithra, Kumar, Pandey, Srivastava, Kumar, Farooqi, Lal, Rai and Chaturvedi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: K. K. Chaturvedi, a2suY2hhdHVydmVkaUBpY2FyLmdvdi5pbg==

These authors have contributed equally to this work

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.