*3.2. Most of the Tiller Genes Are Not Drastically Respond to N Rate*

Recent progress in molecular genetics have deciphered that *MOC1* [9], *LAX1* and *LAX2* [12], *OsFC1*/*OsTB1* [7,8], *APC*/*CTAD1*/*TE* [11,12], and more than 60 other genes [13,14] are involved in tiller number control in rice. However, it is almost unknown which set of genes are involved in tiller regulation under general cultivation conditions, especially which ones are responsible to N rate. In an attempt to answer such a question, we turn to RNA-seq approach to compare the transcriptomic changes in the tissue around the SAM, where the tillers emerge.

Our RNA-seq data point to that drastic N enrichment caused more dramatic transcription change rather than mild N rate (Figure S12). The top categories of those DEGs fall into the metabolic and hormonal/signal transduction pathways. However, all changes of the tiller related genes did not reach the two-fold bar to be a DEG (Table 2). Furthermore, half of these tiller genes were showing opposite change between the two varieties in response to the N rate (Table 2). These suggest us that drastic N rate may trigger a dramatic global transcription shift, but those changes may not necessarily directly relate to tillering. More importantly, change in a gene expression level does not have to surpass the two-fold bar to impact on tiller number.

In addition, previous studies have shown that the tiller genes are mainly analyzed in a certain genetic background. Moreover, the interaction network or regulatory mechanisms of tiller genes at specific stages are not well understood, especially in the early stages of tillering occurrence around the SAM in rice. In this study, we use GWCNA methods to find co-expressed gene modules and explore the relationship between gene networks and phenotypes, as well as core genes in the network. Interestingly, we found that several tiller genes (*HTD2*/*D14*/*D88*/*qPPB3, D10*/*OsCCD8 and NRR*/*CRT*) fall in the green module, which is negatively correlated with tiller regulation (Figure 7B). *HTD2*/*D14*, which is a component of the stragolactones (SLs) signaling pathway, encodes an esterase that inhibits rice branching and negatively regulates rice tiller number [33]. *D88,* function through the MAX/RMS/D pathway, is expressed in most tissues of rice, including leaves, stems and roots and ultimately affect rice plant type by regulating cell growth and organ development. Mutations in *D88* affected the expression of genes involved in tiller formation, including *HTD1, OSH1, D10* and other genes were significantly up-regulated in *d88* mutants [34]. *D10* is a rice ortholog of MAX4/RMS1/DAD1, which encodes a carotenoid cleavage dioxygenase and is involved in the biosynthesis of levodolactone/levylactone derivative SLs [35]. *D10*/*OsCCD8* is involved in the synthesis of rice aboveground branching inhibitors, and transcription of *D10* may be a key step in regulating the branching inhibition pathway. The interaction among *D10*, auxin and cytokinin affects lateral bud elongation in rice [36]. *NRR*/*CCRT* can regulate the structure of rice roots, so that they can better absorb a large number of nutrients and play a negative regulation role in the growth of rice roots. It can respond to the level of photo-contracted compounds and coordinate the expression of genes related to starch synthesis [37,38]. Although *HTD2*/*D88* and *D10*/*OsCCD8* genes belong to the negative regulatory module and they are not key nodes of regulation network, these genes are regulated by the upstream hubgenes. Our results corroborate that these genes closely relate to tillering.

Meanwhile, our results showed that some genes related to tillering, N and other genes fall in other modules and these modules are not closely related to certain traits (tiller, N, etc.). We speculate that these genes may have some novel unknown molecular functions, which of course requires further experimental verification.

Obviously, these differences between varieties are the results of common regulation of many genes. Changes of certain pivotal gene expression will lead to changes in other genes. The difference in these expression patterns led to the difference in the response of N rate between *indica* and *japonica* rice.

In addition, consistent with previous reports on AS of pre-mRNA [39–43], the major category of AS events revealed in our current study is intron retention (RI), irrespective of rice varieties. Despite the trending differences in the percentage of AS events between the two varieties, we have not found significant distribution differences among the categories event-wise. However, as we did not compare the components of AS genes between the groups, it would be interesting to further analyze if there are differences in this aspects of AS events of certain genes between varieties or among N rates [41–43].

#### **4. Materials and Methods**

#### *4.1. Plant Materials*

Two representative rice varieties Nipponbare (NPB, *Oryza sativa* L. subsp. *japonica*) and Yangdao6 (YD6, *Oryza sativa* L. subsp. *indica*) were employed in this study. NPB and 9311 (a sister line of YD6) are widely being used as reference sequencing varieties, in molecular genetic analysis and practical rice production. NPB generally produces more tillers under lower N rate; whereas YD6, bears fewer tillers, yet produces more tillers under higher N rate.

### *4.2. Growth Conditions, N Rate Treatment and Measurement*

The experiments were conducted at the Yangzhou University, Jiangsu Province of China. Plants were grown in plastic pots (29 cm in diameter, 30 cm in height), which were filled with a mixture of soil and vermiculite at 5:1 (*v*/*v*) ratio. The soil type is sandy loam, containing 1.02 g·kg−<sup>1</sup> total N, 22.73 mg·kg−<sup>1</sup> available phosphorus, 49.24 mg·kg−<sup>1</sup> available potassium, and 13.98 g·kg−<sup>1</sup> total organic matter. The soil osmotic conductivity was 0.11 ms·cm−<sup>1</sup> and pH was at 7.66.

Total N fertilizer (urea) enrichment rate was set up at 0 (LN, CK), 9 (MN) and 18 (HN) g N·m<sup>−</sup>2. Basal N fertilizer (50% of total rate, in the form of urea, dissolved into water first then further diluted and applied; the same practice was followed for all other top-dressings) was premixed into the top 10 cm of soil at 3 days prior to seeding. Top dressing was at 2nd and 4th leaf emerging stage, each at 25% of the respective N rate. Twelve sprouting seeds were planted in each pot at an even spacing. After germination, the pot was irrigated with a shallow layer (1–2 cm) of water. Four biological replicates were included, with 16 pots in each treatment.

Enumeration of tiller was conducted at 4th, 6th and 8th leaf emerging stage. Plants were sampled for dry matter measurement, and subsequently for soluble sugar, starch and N content determination. Total soluble sugar and starch detection were by anthrone reagent following Brooks et al. [44]. N determination was by Kjeldhal method with Kjeltec 8400 Analyzer Unit (Foss Analytical AB, Hoganas, Sweden) following the manufacturer's recommendation.

#### *4.3. Samples Preparation for RNA Extraction, cDNA Library Construction and Sequencing*

To investigate the transcriptome changes in response to N enrichment, tissues for RNA isolation were collected when 4th leaf started spreading out. After removing leaves and the out layers of sheaths with a surgical blade, only tissues surrounding the shoot apical meristem (SAM, about 5mm in length) were collected. SAM isolation operation was carried out on the ice, and the SAM tissues were wrapped in aluminum and snap frozen in liquid N2 before transferring to a deep freezer at −80 ◦C till use.

A total of 12 tissue samples were collected for RNA-Seq, representing two varieties (NPB and YD6) and three N levels (LN, MN and HN), with two biological replicates in each combination. Another set of 12 samples was collected for validation purpose. Total RNA was extracted using the RNAiso Plus Total RNA extraction reagent (Cat#9109, TAKARA, Kusatsu, Japan). Qualified total RNA was further purified by RNeasy micro kit (Cat#74004, QIAGEN, Duesseldorf, Germany). The rRNA was removed using Ribo-Zero rRNA Removal Kits (CAT# MRZMB126, Epicentre, Illumina, San Diego, CA, USA). RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Samples with RNA Integrity Number (RIN) ≥ 7 were subjected to the subsequent sequencing reactions. The libraries were constructed by using TruSeq Stranded mRNA LT Sample Prep Kit (CAT#15032612, Illumina, San Diego, CA, USA), and sequenced on the Illumina sequencing platform HiSeq 2500, consequentially 100–150 bp paired-end reads were generated.
