*2.2. Tiller Number Is Related with N and Carbohydrate Content*

Analysis of the N and carbohydrate content of the tissues revealed that MN enrichment increased N content in all tissues (root, sheath, and leaf) at the 4th leaf emergence stage, irrespective of varieties (Figure 2). However, further enrichment of N rate from MN to HN did not alter the N content at this stage. At the 6th leaf stage, the N content showed consistently to be HN > MN > LN (CK), irrespective of tissues or varieties. Yet, the gap of N content diminished at the 8th leaf emerging stage. These indicate that MN could enhance N absorption at the 4th leaf emergence stage, but more N enrichment could not promote N absorption further. The differential gap in N content disappeared as the plant grew for two more leaf-age. It is worthy to notice that the N content dropped faster in YD6 than in NPB from 4th to 6th leaf stage, corroborating that YD6 requires more N input to maintain a similar N content than NPB. This suggests different N management strategy should be applied to different variety for keep N content for tillering.

**Figure 2.** Nitrogen (N) content in root, sheath and leaf organs in different stage of NPB and YD6. N content in root of NPB (**A**) and YD6 (**D**); N content in sheath of NPB (**B**) and YD6 (**E**); N content in leaf of NPB (**C**) and YD6 (**F**).

Water soluble carbohydrate (WSC) (i.e., soluble sugar) content was almost comparable between root and sheath, both being much lower than that in leaf (Figure 3A,C). MN enrichment treatment lowered WSC drastically, whereas HN caused slightly further decrease. Similarly, starch content reduced as N rate increased. Yet, the highest starch content was consistently observed in the sheath, with root being higher than that in leaf (Figure 3B,D), irrespective of growth stages. The correlation between N content and starch, soluble sugar are not significant, irrespective of varieties (Figure 3E,F).

**Figure 3.** Soluble sugar and starch under different N rates in NPB and YD6. Soluble sugar in NPB (**A**) and YD6 (**C**); Starch in NPB (**B**) and YD6 (**D**); Correlation between N content and starch (**E**), N content and soluble sugar (**F**).

The tiller number at the 4th leaf stage was significantly correlated with the tiller number at the 6th leaf stage, but not significantly correlated with those at the 8th leaf emerging stage, respectively (Figure S1A–C), suggesting the association was diluted along with the growth. The tiller number at the 6th leaf stage was positively correlated with the N content in the root (*R<sup>2</sup>* = 0.6303 ns), sheath (*R2* = 0.8783 \*\*) and leaf (*R2* = 0.0.8719 \*\*) at this stage (Figure S1D–F). It displayed a negative correlation with WSC in the sheath (*R<sup>2</sup>* = 0.3970 ns) and leaf (*R<sup>2</sup>* = 0.7956 \*\*) (Figure S1G–I), and also negative correlation with the starch content in the root (*R<sup>2</sup>* = 0.6677 \*), sheath (*R<sup>2</sup>* = 0.4848 ns) and leaf (*R<sup>2</sup>* = 0.7036 \*) (Figure S1J–L). These generally confirm that the tiller number correlated positively with the N content, but negatively with WSC and starch content.

#### *2.3. RNA-Seq Data Quality and Assembly*

To reveal the transcriptomic change responsive to N rates, cDNA libraries from tissues surrounding SAM were constructed and subjected to RNA-Seq analysis. A total of 624 million raw reads were generated from two subspecies *japonica* NPB and *indica* YD6 in combination with three levels of N rate (Table S1). After removal of adaptor sequences, short reads, low quality, and rice ribosome RNA reads, clean reads ratio proved to range from 91.6% to 97.1%, indicating that the sequencing reads were qualified for further analysis.

These clean reads were mapped to the NPB reference genome, the mapped ratio ranged from 97.3% to 98.1% in NPB samples, and from 90.9% to 93.5% in YD6. Since the most recent physical map covers about 97% of the NPB genome [16], our mapped ratio of the transcriptomes in SAM tissues in NPB was generally consistent with it. These also indicate that the alignments were successful. As expected, the transcriptome from YD6 was genetically more heterogeneous to NPB reference genome. Even when trying to map the clean reads of YD6 to a reference genome of 9311, a sister line of YD6, the successfully mapped ratio for YD6 ranged from 93.5% to 95.9%. For comparison consistency, we used the physical genomic map of NPB as a common reference genome in all of the following analysis.

Principal component analysis (PCA) projects the whole transcriptome profile onto a few variables to reflect the distribution across samples, which are often clustered in two or three-dimensional space [17]. PCA revealed that transcriptome profile between YD6 and NPB were more distant than that among N rates (Figure S2A), reflecting that the variety defined the most range of the differences. The variability of LN and MN were less than that of HN in NPB, while MN was less variable than that of LN and HN in YD6. Yet, the general variability in NPB was much less than that of YD6. These suggest that the transcriptomic response to N rate was more variable in YD6 than in NPB, the transcriptomic variability caused by N enrichment also depends on the plant's existing N rate condition, and the transcriptomic differences between varieties were more distant than their responses to N rate.

The correlation analysis of gene expression between the biological replicates reflects the experimental and sequencing reproducibility, whereas the correlation between the groups reflects the variety and N rate effects. The correlation coefficients were determined by the FPKM values of the individual genes to verify the general reproducibility and variability among replicates and groups, respectively (Figure S2B). The correlation coefficient *R* value between the biological replicates ranged from 0.96–0.97 in NPB and 0.95-0.96 in YD6, respectively, exhibiting a remarkably high consistency. This meant this set of RNA-Seq results had very high reproducibility in replicates. Meanwhile, the *R* value between NPB and YD6 at corresponding N rate ranged from 0.89–0.92, displaying certain consistency. Yet, these values were lower than the *R* between the biological replicates, indicating more variability between varieties, as expected. Furthermore, the heat map between the groups is consistent (Figure S2C,D). These ensure that the RNA-Seq generated a highly reproducible expression data among the replicates and a proper variability across the comparison groups, which corroborated with PCA results as well as phenotypic performance.
