*2.9. Construction of the Gene Co-Expression Networks and Identification of Candidate Hub Genes*

After screening for the highest significant correlated modules, we selected two modules (one positive correlation module and one negative module) to construct the gene expression networks (Figure 8, Figures S6–S8). Since the scale of N in stem and N in leaf, dry weight in stem and in leaf share the same selected modules, totally seven most relevant modules (light yellow, green, blue, tan, red, brown and pink module) were chosen for the further analysis. The association degree of eigengenes within the module was used to construct the gene expression network (Figures 8A,B, S5A,B, and S6A,B). The hub genes are selected by KME values through their closest connections in the gene network. The top 20 node genes of each module were screened to generate a network map (Figures 8C,D, S5C,D, and S6C,D).

**Figure 8.** Co-expression network analysis of tiller number related modules. (**A**,**B**) Gene co-expression networks of positively correlated lightyellow module (**A**) and negatively correlated green module (**B**) visualized using Cytoscape software platform. The circle size of and color depth indicate the degree of connectivity; (**C**,**D**) The correlation networks of top 20 nodes in lightyellow module (**C**) and green module (**D**). The color depth represents the number of associated nodes.

To screen for candidate hub genes, the top ten connected genes in the specific modules were picked out by their KME value (Table 1, Table S15). The tiller number relate positively with the lightyellow module (M18) and negatively with the green module (M3), respectively (Table 1). In the light yellow module, the top ten genes encode HEV3—Hevein family protein precursor, BBTI5, BBTI8, VP15, ribonuclease T2 family, pyridoxal-dependent decarboxylase protein, cysteine-rich receptor-like protein kinase, and kinesin motor domain containing protein. The correlation network of the light yellow module is shown in Figure 9A. Genes encoding the aforementioned proteins were identified as key candidate hub genes for the light yellow module (M18) (Figure 8C). In the green module (M3), the top ten genes encode MATE efflux protein, ethylene-responsive protein, protein kinase, GDSL-like lipase/acylhydrolase, aspartic proteinase oryzasin-1 precursor, eukaryotic translation initiation factor 1A, thaumatin, tetraspanin family protein, AGAP002737-PA, and RFC5. Genes encoding these proteins were identified as candidate hub genes for the green module (Figure 8D). Similarly, multiple functional proteins were enriched in this module, indicating that these regulatory networks may play a pivotal role in regulation of the scale of tillers.


**Table 1.** Candidate hub genes related to the tiller numbers, N rate in leaf and in stem.

The N content in stem and in leaf relate positively with the blue module and negatively to the green module, respectively (Table 1). In the blue module, the top ten genes encode RNA recognition motif containing protein, ribosomal L18p/L5e family protein, cytochrome P450, ribonuclease T2 family, heat shock 22 kDa protein, CCT motif family protein, NB-ARC domain containing protein, and transcription factor etc. The correlation network of the blue module is shown in Figure S5A. This indicates that these genes may play a crucial role in regulation of the scale of N in stem and in leaf. Interestingly, the negative module falls into the same module with the scale of tillers. In addition, the co-expression network and candidate hub genes in other modules are shown in Figures S6–S8 and Table S15.

#### *Int. J. Mol. Sci.* **2019**, *20*, 5922


**Figure 9.** Gene expression of top ten node hub genes of eight related modules. (**A**) Heatmap showing the expression profiles of the top ten node hub genes in each module; (**B**) Fold changes (log2) of the top ten node hub genes in lightyellow and green module. Red and blue color represent scale of up-and down-regulation, respectively.

It is noteworthy that we find four hub genes overlap with DEGs (between varieties) in the turquoise module (Table S15), which negatively correlated with the scale of dry weight in leaf. These four genes are LOC\_Os12g07830, LOC\_Os06g15570, LOC\_Os08g25050 and LOC\_Os06g12170. The results suggest that these genes in the regulatory network may play a core role in differential dry weight accumulation between varieties.

The heatmap and changes in the expression levels of the candidate hub genes from the selected modules are shown in Figures 9 and S9. The expression level of these hub genes are contrastingly different between the varieties, indicating that the hub genes may play a crucial role in their differential response to N rate between the varieties.
