*4.6. Quantitative Real Time RT-PCR Validation*

To validate the RNA-Seq data, twenty-one DEG at different expression level were selected to confirm their expression in the corresponding samples by quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) with a BioRad CFX-96 system (BioRad, Hercules, CA, USA), following the method in [50]. Briefly, one μg of total RNA from the same batch of RNA for high throughput RNA-Seq was used for the first strand of cDNA synthesis using iScript (Cat#1708891, BioRad, Hercules, CA, USA) according to the supplier's protocol. The PCR reaction was conducted in a total volume of 12 μL reaction mix, with one μL of cDNA template, 400 nM forward primer, 400 nM reverse primer and six μL of SsoFast EvaGreen Supermix (Cat. #1725200, Bio-Rad, Hercules, CA, USA), and three technical replicates. Two tubulin genes (LOC\_Os11g14220 and LOC\_Os03g51600) were deemed as the reference genes. Gene-specific primer sequences for qRT-PCR are listed in Table S2. Quantification was determined by BioRad CFX manager software (V3.1). Fold change relative to control level was determined by the 2−ΔΔC*<sup>t</sup>* method [50]. PCR amplifications of each sample were used in triplicate.

#### *4.7. Co-Expression Network Analysis for Construction of Modules*

WGCNA (v1.29) package in R was performed to construct the gene co-expression regulation network [20]. Detailed analysis procedures and methods were followed in accordance with Zhang et al. [51]. Among the 37,824 genes, 21,700 genes with an averaged FPKM from three replicates > 1 were used for the WGCNA unsigned co-expression network analysis. Through testing the independence and the average connectivity degree of different modules with different power value, the appropriate power value in this study was determined as seven. The modules were obtained by the automatic network construction function with default parameters in the WGCNA software package. The correlation between the modules and traits were calculated by the Pearson method using *blockwiseModules* function. The top ten genes with maximum intra-modular connectivity were considered as "highly connected gene" (hubgene) [52]. The top 20 genes including candidate hub gene network was visualized by the Cytoscape (version 3.7.2) [53].

#### *4.8. Statistical Analysis of Genes Expression Data in a Specific Pathway*

Multivariate analysis procedure of the general linear model (GLM) method from the IBM SPSS software (version 22, IBM, Armonk, NY, USA) was used in the ANOVA comparison of selected gene expression differences. Correlation analysis was made by the correlate procedure in the SPSS software.

#### **5. Conclusions**

This is a pioneer study to reveal the transcriptomic changes of SAM tissue in response to N rate between *indica* and *japonica* rice subspecies, especially for cultivars widely used in the production. N rate influenced tiller number in a different pattern between varieties, with NPB being more sensitive to N enrichment, and YD6 being more tolerant to N rate. Tiller number was positively related to N content in leaf, culm and root tissue, but negatively related to the soluble carbohydrate, irrespective of variety. Higher N enrichment brought more drastic transcription change than moderate N rate; however, varietal background dominated the differences. For the reported 65 tiller genes, less than half of them showed decent expression in SAM at tiller starting phase; among them only nineteen being significantly influenced by N rate, and two genes showing significant interaction between the N rate and variety. GO analysis revealed that the majority of these common DEGs are involved in general stress responses, stimulus responses, and hormonal signaling process. WGCNA network identified specific modules that are associated with the phenotypic traits and candidate hub genes for each module. Several genes associated with tillering and N content fall on certain most relevant modules. These results help us understand the complexity regulatory mechanisms involved in *indica* and *japonica* rice response to N rate.

**Supplementary Materials:** Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/23/ 5922/s1.

**Author Contributions:** Conceptualization, X.Z., J.Z., J.H. and Y.Y.; Data curation, N.H. and L.M.; Formal analysis, X.Z. and N.H.; Investigation, L.M., M.L., Y.G. and C.C.; Methodology, J.J.; Project administration, A.L. and Y.W.; Resources, Y.W., J.H. and Y.Y.; Software, Z.Y.; Validation, S.Y. and G.D.; Visualization, Y.Z.; Writing–original draft, X.Z.; Writing–review & editing, Y.Y.

**Funding:** This work was funded by the National Natural Science Foundation of China (31571608), Natural Science Foundation of Jiangsu Province (BK20151311), Major Basic Research Project of the Natural Science Foundation of the Jiangsu Higher Education Institutions (15KJA210003), the Open Funds of the Key Laboratory of Plant Functional Genomics of the Ministry of Education (ML201903), Jiangsu Agriculture Science and Technology Independent Innovation Fund (CX(19)3060), and the Social Development Project of Yangzhou (YZ2018067).

**Acknowledgments:** We gratefully thank Chenwu Xu (College of Agriculture, Yangzhou University, Yangzhou, China) for meaningful help.

**Conflicts of Interest:** The authors declare no conflict of interest.
