Genome-wide transcriptome profiling provides insights into panicle development of rice (Oryza sativa L.)

Genome-wide transcriptome profiling provides insights into panicle development of rice (Oryza sativa L.)

Accepted Manuscript Genome-wide transcriptome profiling provides insights into panicle development of rice (Oryza sativa L.) Shanwen Ke, Xin-Jiang Li...

2MB Sizes 0 Downloads 59 Views

Accepted Manuscript Genome-wide transcriptome profiling provides insights into panicle development of rice (Oryza sativa L.)

Shanwen Ke, Xin-Jiang Liu, Xin Luan, Weifeng Yang, Haitao Zhu, Guifu Liu, Guiquan Zhang, Shaokui Wang PII: DOI: Reference:

S0378-1119(18)30765-0 doi:10.1016/j.gene.2018.06.105 GENE 43038

To appear in:

Gene

Received date: Revised date: Accepted date:

6 February 2018 26 April 2018 28 June 2018

Please cite this article as: Shanwen Ke, Xin-Jiang Liu, Xin Luan, Weifeng Yang, Haitao Zhu, Guifu Liu, Guiquan Zhang, Shaokui Wang , Genome-wide transcriptome profiling provides insights into panicle development of rice (Oryza sativa L.). Gene (2018), doi:10.1016/j.gene.2018.06.105

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Title: Genome-wide transcriptome profiling provides insights into panicle development of rice (Oryza sativa L.)

Author names: Shanwen Ke, Xin-Jiang Liu, Xin Luan, Weifeng Yang, Haitao Zhu, Guifu Liu,

Affiliation: Key

Laboratory

for

Conservation

and

Utilization

SC

State

RI

PT

Guiquan Zhang*, Shaokui Wang*

of

Subtropical

Agro-Bioresources, Guangdong Key Laboratory of Plant Molecular Breeding,

NU

College of Agriculture, South China Agricultural University, Guangzhou 510642,

MA

China.

Corresponding author: Corresponding author.

D

*

PT E

E-mail address: [email protected] (G. Zhang), [email protected] (S.

AC

CE

Wang).

1

ACCEPTED MANUSCRIPT Abstract Panicle architecture is an important component of agronomic trait in rice, which is also a key ingredient that could influence yield and quality of rice. In the panicle growth and development process, there are a series of complicated molecular and cellular events which are regulated by many interlinking genes. In this study, to

PT

explore the potential mechanism and identify genes and pathways involved in the formation of rice panicle, we compared the transcriptional profile of rice panicles

RI

(NIL-GW8 and NIL-gw8Amol) at three different stages of panicle development: In5 (formation of higher-order branches), In6 (differentiation of glumes) and In7

SC

(differentiation of floral organs). A range of 40.5 to 54.1 million clean reads was aligned to 31,209 genes in our RNA-Seq analysis. In addition, we investigated

NU

transcriptomic changes between the two rice lines during different stages. A total of 726, 1,121 and 2,584 differentially expressed genes (DEGs) were identified at stages

MA

1, 2 and 3, respectively. Based on an impact analysis of the DEGs, we hypothesize that MADS-box gene family, cytochrome P450 (CYP) and pentatricopeptide repeat

D

(PPR) protein and various transcription factors may be involved in regulation of

PT E

panicle development. Further, we also explored the functional properties of DEGs by Gene Ontology analysis, and the results showed that different numbers of DEGs genes were associated with 53 GO groups. In KEGG pathway enrichment analysis, many

CE

DEGs related to biosynthesis of secondary metabolites and plant hormone signal transduction, suggesting their important roles during panicle development. This study

AC

provides the first examination of changes in gene expression between different panicle development stages in rice. Our results of transcriptomic characterization provide important information to elucidate the complex molecular and cellular events about the panicle formation in rice or other cereal crops. Also, the findings will be helpful for the further identification of the genes related to panicle development.

Keywords Rice; Panicle; RNA-Seq; Transcriptome; GW8 2

ACCEPTED MANUSCRIPT Highlights Tanscriptome analysis was used to study the panicle differentiation of rice. Compare gene expression profile of two rice near isogenic lines during different panicle development stages. Obtain a list of candidate genes related to the formation of rice panicle involved in

RI

PT

OsSPL16.

Abbreviations

SC

COG, Cluster of orthologous groups of proteins; CYP, Cytochrome P450 monooxygenase; DEGs, Differentially expressed genes; FDR, False discovery rate;

NU

FPKM, Fragments per kilobase of transcript per million fragments mapped; GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; NIL, Near isogenic

MA

lines; Nr, Non-redundant protein database; PPR, pentatricopeptide repeat protein; qRT-PCR, Quantitative real-time PCR validation; QTL, Quantitative trait locus; SBP,

AC

CE

PT E

D

Squamosa promoter binding protein.

3

ACCEPTED MANUSCRIPT 1. Introduction Rice (Oryza sativa L.) is one of the most important cereal crops in the world. Benefit from the green revolution and the application of heterosis, the production of rice has achieved the double increase in the last 50 years, and even reached the increase of more than three times in some countries and regions (Xing and Zhang,

PT

2010). But it is predicted that the world food production must rise by 50% to meet the demand of the increasing population in 2050 (Ashikari et al., 2005). Panicle

RI

architecture is an important component of agronomic trait in rice, which determines yield of rice. The formation of panicle is a systemic process including the

SC

development of axillary meristems, the establishment of inflorescence and the formation of grains.

NU

Although researchers have studied the formation of grains for the whole last century, but the study of the molecular mechanism of grain size determination was

MA

only about 20 years, and the results mainly focus on model plants (Sundaresan, 2005). Grain affects not only the production but the appearance. Improvement of grain shape

D

is an effective way to increase yield and quality for rice. With the acceleration of

PT E

mutant research and the use of the advanced biological technologies, as well as the large-scale development of molecular markers, many quality trait genes associated with grain shape have been identified and cloned (Luo et al., 2013; Si et al., 2016).

CE

Some of them are negative regulators of grain length. For example, GS3, one major QTL for grain length, encodes a putative membrane protein that regulates grain

AC

growth by limiting cell proliferation (Fan et al., 2006; Mao et al., 2010). GL3.1 regulates grain size and yield by regulating the cyclin-T1;3 that influences the replication of chromosomes and promote cell division, leading to shorter rice grain (Qi et al., 2012). Studies have shown that GL3.1 influences the grain size and yield mainly through changing cell number of glume to regulates its length, on the other hand, it also affects the filling rate of rice (Zhang et al., 2012). GS2 functions as a positive regulator of grain length through regulating cell division and subsequent elongation (Duan et al., 2016). In addition, a few other genes have been studied as the regulators of grain width. GW2 is the first cloned gene which controls the grain width 4

ACCEPTED MANUSCRIPT and weight; it may be involved in the degradation of protein to promote the process of cell division, and then influence grain weight and yield by regulating the rice glume shell size (Song et al., 2007). The gene, qSW5, was discovered through gene mapping analysis; it regulates grain width by influencing the clever shell cell number (Shomura et al., 2008). Its allele GW5 was also cloned; it affect the grain width likely through

PT

participating in protein degradation pathways (Weng et al., 2008). GS3, GL3.1, GS2, GW2, qSW5 and GW5 reflect the negative relation of sequence. However,

RI

GS5 positively regulates rice grain weight; overexpression of GS5 increases grain size by promoting cell proliferation (Li et al., 2011). Another major QTL for rice weight

SC

and length (TGW6) promotes IAA-glucose hydrolyzed into free IAA and glucose, and participates in the carbohydrate transportation between source and sink (Ishimaru et

NU

al., 2013).

On the other hand, the number of branch and spikelets per panicle is one of the

MA

key factors affecting rice yield. The spikelet is the basal unit of inflorescence in grasses, and it is crucial for reproductive success and cereal yield. Although the

D

growth of plant organs, like panicles, often affected by the external environment

PT E

signals, there is still plenty of evidence to suggest that the morphology of plant organs is mainly affected by the internal signals which is controlled by many genes (Liu et al., 2002). Gn1a controls branch differentiation; the decreased expression level of Gn1a

CE

leads to the accumulation of cytokinins in inflorescence meristem, resulting in the increase of spikelets number per panicle (Ashikari et al., 2005). DEP1 is postulated to

AC

function as a negative regulator of branch differentiation, when in function-defect of DEP1, branch number and spikelets number increases significantly (Huang et al., 2009). PROG1 is also a down-regulator for the number of branch and spikelets, and it is mainly expressed in axillary bud organization (Tan et al., 2008). A few positive regulatory factors have also been found. For instance, FZP plays an important role in the forming process of spikelet meristems (SMs) from spikelet-pair meristems (SPMs) which may participate in the development of branch meristems (Komatsu et al., 2003). Besides, there are many positive regulatory factors in the process of the development of young panicles, like APO1, LAX1, LAX2, RCN1, RCN2 (Nakagawa et al., 2002; 5

ACCEPTED MANUSCRIPT Terao et al., 2010; Tabuchi et al., 2011). In rice, some SBP (Squamosa promoter Binding Protein) protein family members are tightly tied to the development of branch. Such as OsSPL14, its high-level expression leads to the increase of branch number and grain number (Jiao et al., 2010). OsLG1, encoding a SBP protein, regulates the panicle development by influencing the

PT

parenchyma cells of the gaps between the cobs and branches (Ishii et al., 2013). In our previous work, we have cloned a new gene GW8 which belongs to the family of SBP

RI

transcription factors. It positively regulates cell proliferation, improving not only the 1000-grain weight but also the rice quality (Wang et al., 2012). Based on the analysis

SC

of the single segment substitution lines (SSSLs), we found a GW8 allele from an indica rice (Amol3), named gw8Amol. Compared with the NIL-GW8, the grains

NU

showed high quality and slender shape without yield difference in NIL-gw8Amol because of the increase in the number of branch and grains per panicle (Fig. 1) (Wang

MA

et al., 2012). Our previous studies have shown that the expression of Gn1a, Lax1, FZP, OsCLV1, OsCLV2, OsCLV3, RCN1, RCN2, DEP1 and OsSPL14 have no significant

D

difference in panicle at the early stage of inflorescence development between

PT E

NIL-GW8 and NIL-gw8Amol plants (Wang et al., 2012). We suspect that GW8 regulates the panicle development and differentiation through some unknown metabolic pathways. It is interesting to identify and characterize the regulators of GW8 (ROGs)

CE

which is an essential gene affecting the branch differentiation and regulation of the grain number per panicle in rice.

AC

The completion of the sequencing of rice genome and the development of high-density molecular markers provide a new opportunity for the research of the genome of rice and other crops (Devos, 2005; Ashikari and Matsuoka, 2006; Xing and Zhang, 2010). In recent years, with the rapid development of next-generation sequencing technology, mRNA sequencing (RNA-seq) technology has become a useful means for the study of transcription profiles (Wang et al., 2010). In this study, a wide grain rice line (NIL-GW8) and a slender grain rice line (NIL-gw8Amol) were used to compare gene expression changes in panicles by deep RNA sequencing. We did the transcriptome analysis of three developmental stages of panicles in these two rice 6

ACCEPTED MANUSCRIPT lines. The identification and functional validation of ROGs may help us to better understand the working mechanism of GW8 in the regulation of cell proliferation in rice and eventually create new varieties with higher yield and better quality.

2. Materials and methods

PT

2.1. Plant materials Two indica rice (O. sativa L. ssp. indica) near-isogenic lines (NIL) NIL-GW8 and

RI

NIL-gw8Amol in the HJX74 genetic background were used for the transcriptome sequencing analysis in this research. NIL-gw8Amol contains a very short chromosome

SC

segment inherited from Amol3 (Wang et al., 2012), carrying an allele of GW8 (gw8Amol) which has a mutational OsmiR156 binding site. And NIL-gw8Amol produces

NU

slender grains and more branches and spikelets per panicle than NIL-GW8. The two rice lines were grown in Teaching and Experimental Farm in South China

differentiation period

of

MA

Agricultural University (113°36'E, 23°16'N). We sampled the panicles in the primary branch and secondary branch.

A total of

six

D

inflorescence samples (three stages: In5 (formation of higher-order branches), In6

PT E

(differentiation of glumes) and In7 (differentiation of floral organs) from the two NILs) were collected for RNA-seq according to Ikeda et al. (Ikeda et al., 2004). These collected samples were immediately frozen in liquid nitrogen and stored at -80℃

CE

until the total RNA extraction. The panicle organs which used in this study are very small, and the differentiation periods are relatively shorter with one to three days

AC

(Ikeda et al., 2004). For biological replicates, every sample of each developmental stage was a mixture of fifty panicles from at least ten individual plants. The growth and development status of these sampled panicles in a specific development period were highly consistent.

2.2. RNA isolation and cDNA library construction Total RNA was extracted from panicles and spikelets using TRIzol® RNA Purification Kit (Invitrogen) according to the manufacturer’s protocol, and electrophoresed in 1.0 % agarose gel. Nanodrop 2000 Spectrophotometer, Qubit® 2.0 7

ACCEPTED MANUSCRIPT Spectrophotometer

(Life

Technologies,

Carlsbad,

CA)

and

Aglient

2100

Bioanalyzer® (Agilent Technologies, Palo Alto, Calif.) were used to measure and test the purity, concentration and integrity of RNA samples, respectively. The enrichment of the poly (A)-containing mRNA was performed by using Oligo (dT) magnetic beads. The mRNA was fragmentated and used to synthesize double-stranded cDNA. The Agencourt AMPure XP beads (Beckmann Coul-ter, Pasadena, CA) were used to

PT

purify dsDNA. After that, the ends of purified ds-cDNA were repaired. Next, the ‘A’

RI

tail and paired-end adapters for sequencing were ligated to the end of cDNA. In the end, different sizes of fragments were selected using AMPure XP beads for PCR

SC

amplification to create the final RNA-seq library. When the cDNA library construction is finished, the concentration and insert size of the libraries were tested

NU

using Qubit2.0 and Agilent 2100, respectively. To ensure the quality of library, using the qPCR accurately quantify the effective concentrations of the libraries above. After

MA

complete the testing, the libraries were sequenced on Illumina HiSeq2500 Sequencer. These sequence data were deposited in an NCBI database and are accessible via

PT E

D

BioProject ID: PRJNA439157.

2.3. Quantitative real-time PCR validation The validation of the RNA-seq results was performed through ten genes with

CE

different expression profiles by quantitative RT-PCR. The read counts and the primer sequences of the ten genes were provided in Table S6 and Table S7, respectively.

AC

Total RNA was extracted references the RNA-seq experiments and first-strand cDNA was generated using a TransScript One-Step gDNA Removal and Cdna Synthesis SuperMix kit (TransGen Biotech). Quantitative RT-PCRs were carried out using the ABI Step One Plus Real-Time PCR System (Applied Biosystems) with the SYBR Premix Ex Taq™ RT-PCR kit (Takara). The rice gene Actin1 (LOC_Os03g50885) was used to standardize the relative expression levels.

2.4. Quality analysis and mapping reads to the reference genome The software Phred was used to evaluate the probability of incorrect base calling 8

ACCEPTED MANUSCRIPT according to phred quality score (Ewing et al., 1998). Before mapping reads to the reference genome, read quality control was conducted as follows: first, the sequencing adapters and primer sequences were removed; then, all the low-quality raw reads were removed to get the clean reads. Then all the clean reads were mapped to the rice reference

genome

(Os-Nipponbare-Reference-IRGSP-1.0,

PT

ftp://ftp.ensemblgenomes.org/ pub/plants/release-24/fasta/oryza_sativa/) (Kim et al.,

RI

2013) by software TopHat V2 (Langmead et al., 2009).

2.5. Analysis of gene expression

SC

Analysis of gene expression was performed based on Negative Binomial Distribution. The Cuffdiff components of Cufflinks software were used to do

NU

quantitative analysis of transcripts and calculate genes expression level according to reads mapped to specific genes. FPKM (Fragments Per Kilobase of transcript per

MA

Million fragments mapped) value was directly used for assessing the expression

D

abundance of each gene with high sensitivity across six orders (10-2 to 104).

PT E

2.6. Identification of differentially expressed genes We applied EBSeq package to identify differentially expressed genes (DEGs) among different panicle developmental stages of the two rice varieties. Genes with

CE

Fold Change ≥ 2 and FDR (false discovery rate) < 0.01 were selected as DEGs. The log2FC (fold changes) was used to indicate the significance of expression difference

AC

of each gene. The test corrections of P-values were performed using the Benjamini-Hochber False Discovery Rate method. The numbers of DEGs between sample pairs were counted and plotted in a Venn diagram. For Gene ontology (GO) enrichment analysis, BLAST2GO Ver. 2.7.2 was used for the detected DEGs according to biological process, molecular function and cellular component ontologies. GO terms with a FDR cut-off of 0.01 were considered as significantly enriched. Then the DEGs were subjected to search against COG (Cluster of Orthologous Groups of proteins) database to predict and classify their possible functions. After that, the pathway enrichment and annotation analysis of DEGs were 9

ACCEPTED MANUSCRIPT carried out by comparing to the KEGG(Kyoto Encyclopedia of Genes and Genomes) database using Cytoscape software. Finally, the DEGs significantly enriched metabolic pathways or signal transduction pathways were identified with ClueGO.

3. Results

PT

3.1. Illumina sequencing and alignment against the reference genome The sequencing of mRNA generated six transcriptome datasets from two NILs

RI

(NIL-GW8 and NIL-gw8Amol) at three different panicle developmental stages [W1: stage 1 (In5) in NIL-GW8; W2: stage 2 (In6) in NIL-GW8; W3: stage 3 (In7) in

SC

NIL-GW8; S1: stage 1 (In5) in NIL-gw8Amol; S2: stage 2 (In6) in NIL-gw8Amol; S3: stage 3 (In7) in NIL-gw8Amol]. W1, W2, W3, S1, S2 and S3 libraries generated around

NU

46.5, 53.1, 54.1, 40.5, 40.8 and 46.1 million high-quality clean reads, respectively, each of which was of a Q30 (an error probability of 1%) value of at least 89.48%

MA

(Table 1). All the reads were mapped to the rice reference genome by TopHat Ver. 2.0.11. A total of 227.1 million reads with a mean mapping rate of 80.88 % (ranging

D

from 79.88 to 82.35% for the six libraries) were mapped to the genome. Among of the

PT E

mapped reads, 77.55-80.01% were matched to unique genomic locations of the reference genome. 17.65-20.12% of the total reads were unmapped (Table 1).

CE

3.2. Global analysis of gene expression In the study, a total of 5,971 (W1), 5,899 (W2), 6,206 (W3), 5,737 (S1), 5,606 (S2)

AC

and 6,147 (S3) genes with different length (most of them were no less than 200 bp) were obtained from the six samples, respectively. Of the total 31,209 genes, 4,693 (15.03 %) genes longer than 3,000 bp accounts for the biggest proportion, and only 28 (0.08 %) genes with shorter than 200 bp (Fig. 2). More than half of those genes (51.45 %) were longer than 1,500 bp (Fig. 2). To get a better visualization of the gene expression in the six libraries, the Circos diagram was generated using the Circos software package. In each library, fragments per kilobase of transcript per million fragments mapped (FPKM) was used as the index of gene expression. As shown in Fig. 3, the expression of the genes had 10

ACCEPTED MANUSCRIPT apparent differences among different chromosomes. The genes with high expression were distributed mainly over Chr1, Chr2 and Chr3, while the genes on Chr11 and Chr12 showed low expression compared with other chromosomes in the six libraries. We obtained a large number of genes after removing partially overlapping sequences, and was able to get plenty of data for the analysis of rice panicle

characters.

Their

expression

in

the

three

panicle

PT

development. In the six libraries, the numbers of expressed genes displayed similar developmental

stages

RI

(2.0-millimetre-long panicle, 3.0-millimetre-long spikelet and 5.0-millimetre-long spikelet, respectively) were summarized in several Venn diagrams. Among these

SC

genes in NIL-GW8, besides the 24,133 genes expressed at all three panicle developmental stages, 563 genes were expressed in stage1 (W1) and stage2 (W2)

NU

simultaneously, 583 genes co-expressed in stage 2 and stage 3 (W3), and 563 genes co-expressed in stage1 and stage3. However, there were 313, 347 and 1,907genes

MA

specifically expressed in stage1, stage2 and stage3, respectively (Fig. 4). In addition, it also had some uniquely expressed genes in the two types of panicles from

D

NIL-GW8 and NIL-gw8Amol. The number of variety-specifically expressed genes was

PT E

560 (NIL-GW8) and 1,778 (NIL-gw8Amol) at stage 1, 710 (NIL-GW8) and 1,446 (NIL-gw8Amol) at stage 2, 960 (NIL-GW8) and 1,181 (NIL-gw8Amol) at stage 3, respectively (Fig. 4). These uniquely expressed genes may play an important role in

CE

the formation process of panicles.

AC

3.3. Identification of DEGs The growth and development of rice is a complicated procedure,there are many differentially expressed genes were identified between a rice line and its NIL (Moumeni et al, 2015; Tariq et al, 2018). Also, many genes involved in the development of rice panicle, especially the inflorescence and spikelet morphogenesis (Meng et al, 2017). The development period we chose involves the critical period of panicle morphogenesis, rapid changes in gene expression in rapidly differentiated organization (Ikeda et al., 2004; Itoh et al., 2005). On the other hand, every plant sample was from different individual plants, the individual differences were enriched 11

ACCEPTED MANUSCRIPT and amplified. So, we obtained a large number of differentially expressed genes (DEGs) in this study. Further, the differences in gene expression among the different developmental stages in NIL-GW8 and NIL-gw8Amol were analyzed. To identify DEGs, the false discovery rate (FDR) threshold was set to 0.01 and a threshold of 3-fold expression

PT

change was used as the selection criteria. All of the differentially expressed genes between two samples during different panicle developmental stages in NIL-GW8 and

RI

NIL-gw8Amol were listed in Table S1. As shown in Fig. 5, we found that a total of 858 (142 up-regulated genes, 716 down-regulated genes) DEGs in NIL-GW8 and 253 (79

SC

up-regulated genes, 174 down-regulated genes) DEGs in NIL-gw8Amol between stage 1 and stage 2, indicating that a lot more down-regulated genes were identified in

NU

NIL-GW8 than NIL-gw8Amol between stages 1 and 2. Between stage 2 and stage 3, 4,265 (2,835 up-regulated genes, 1,430 down-regulated genes) and 4,295 (2,933

MA

up-regulated genes, 1,362 down-regulated genes) DEGs were found in NIL-GW8 and NIL-gw8Amol, respectively. Compared to the early stage of the panicle development,

D

there were more significantly changed genes detected in the later period, suggesting

PT E

that the GW8 gene was specifically involved in the developmental processes at the stage of 2 and 3. There are a large amount of DEGs identified between the same stage of NIL-GW8 and NIL-gw8Amol. 726 (460 up-regulated genes, 266 down-regulated

CE

genes), 1,121 (1,082 up-regulated genes, 319 down-regulated genes) and 2,584 (1,160 up-regulated genes, 1,424 down-regulated genes) DEGs were identified for stages 1, 2,

AC

and 3, respectively (Fig. 5). The numbers of DEGs in stage 2 and stage 3 were remarkably greater than stage 1 indicating the involvement of complex developmental events in the development. Besides, we also noticed that the number of down-regulated genes was higher than up-regulated genes only in the comparison of W3-S3 (stage 3). As shown in Table 2, some DEGs with significant expression difference were detected between NIL-GW8 and NIL-gw8Amol at different stages. Among these genes, some of them were regarded as new genes which could not be mapped to the rice reference genome. To obtain the annotation information of the new genes, these genes 12

ACCEPTED MANUSCRIPT were annotated to Nr, Swiss-Prot, GO, COG, KOG, Pfam and KEGG databases, respectively. Some genes, such as Os02g0127600, which belong to down-regulated genes at stages 1 and 2, showed significantly higher expression at stage3, and some up-regulated genes at stages 1 showed significantly lower expression at stages 2 and 3

PT

(Table 2).

3.4. Genes related to the development of panicle

RI

In the identified DEGs, some of them are important for the development of panicle branches and spikelets. APO2/RFL, playing an important role in controlling and

flower

development

through

interaction

SC

inflorescence

with

APO1

(Ikeda-Kawakatsu et al., 2012), can modulate the differentiation of axillary meristem

NU

by regulating the expression of LAX1 (Deshpande et al., 2015). They all exhibited different expression trends during stage transition from 2 to 3 in NIL-GW8 and

MA

NIL-gw8Amol panicles. Unexpectedly, APO2/RFL, APO1 and LAX1 were all down-regulated with -7.677, -5.396, -6.604 in NIL-GW8 and with -7.782, -5.963,

D

-4.901 in NIL-gw8Amol, respectively. Among them, APO2/RFL was the highest

PT E

down-regulated gene in NIL-GW8. TAW1, a special regulator of rice inflorescence architecture, was highly expressed in shoot apical meristems, inflorescence meristem and branch meristem (Yoshida et al., 2013). Its expression was declined with -6.542

CE

and -6.477 from stage 2 to 3 in NIL-GW8 and NIL-gw8Amol. However, the other genes associated with panicle development, like Gn1a, LAX2, OsREL2, Hd1, DEP1 and

AC

PROG1 (Ashikari et al., 2005; Tan et al., 2008; Huang et al., 2009; Endo-Higashi and Izawa, 2011; Tabuchi et al., 2011; Yoshida et al., 2012), they showed no significant expression difference among different stages.

3.4.1. MADS-box genes Previous research about rice has shown that MADS-box genes involved in the development of panicle branches and grains, like OsMADS34 (Zhang et al., 2016), regulated by SPL14 (Wang et al., 2015). In this study, OsMADS34 was significant differentially expressed between stage 2 and 3 in NIL-GW8 and NIL-gw8Amol (Table 13

ACCEPTED MANUSCRIPT 3). In addition, we found thirteen up-regulated and seven down-regulated MADS-box genes (|log2FC| > 2) in the two rice lines. Interestingly, the elevated level of eleven up-regulated genes was higher in NIL-GW8 than in NIL-gw8Amol. Among them, OsMADS13, OsMADS3, MADS6 and OsMADS58, the coordinative effects of which determine the development of ovule and the floral meristem identities (Li et al., 2011),

PT

were significantly up-regulated in this study. OsMADS58, known to play a vital role in early rice stamen development (Chen et al., 2015), was up-regulated in

RI

NIL-gw8Amol with 3.341, but was not increased obviously in NIL-GW8. Conversely, OsMADS29, a gene closely related to embryo and endosperm development in rice

SC

(Nayar et al., 2013), was up-regulated in NIL-GW8 with 3.396, but was not increased obviously in NIL-gw8Amol. In addition, OsMADS1 and OsMADS17 were up-regulated

NU

and OsMDP1 was down-regulated only in NIL-GW8 from stage 1 to 2. Previous studies have shown that OsMADS1 associated with spikelet development (Wang et al.,

MA

2017), and OsMADS17 is the direct target gene of OsMADS1 (Hu et al., 2015). OsMDP1, has a negative regulatory role in brassinolide (BR) signaling (Duan et al.,

D

2006). We speculate that the different expression patterns of these MADS-box genes

3.4.2. CYP genes

PT E

were mainly because of the different function of GW8 in NIL-GW8 and NIL-gw8Amol.

CE

The expressions of a series of genes named cytochrome P450 (CYP) were differentially expressed in NIL-GW8 and NIL-gw8Amol during stage transition from 2

AC

to 3 (Table 4). CYP comprises one of the largest families in plant enzyme protein, participating in a suite of developmental and physiological processes, including reproduction (Bolwell et al., 1994). In this research task, fourteen CYP genes were up-regulated (log2FC > 3.0) between stages 2 and 3 in the two rice lines. It should be noted that two CYPs(Os03g0332100 and Os05g0482400)were involved in the gibberellin (GA) biosynthesis and signal transduction in rice (Luo et al., 2006; Magome et al., 2013). However, the other gene PLA1 (Os10g0403000), encoding a cytochrome P450 CYP78A11 which acts downstream of GA signaling pathway (Mimura et al., 2012) and participates in the development of primary branch (Miyoshi 14

ACCEPTED MANUSCRIPT et al., 2004), was highly down-regulated in the two rice lines from stages 2 to 3. Furthermore, Os04g0570500 was up-regulated with 5.777 only in NIL-GW8, Os11g0289700 was up-regulated with 5.655 only in NIL- gw8Amol. During stage conversion from 1 to 2, there were no differently expressed CYP genes in NIL-GW8, whereas only one CYP gene (Os01g0227400, log2FC 4.742) was up-regulated in

PT

NIL-gw8Amol. The identification of these differentially expressed cytochrome P450

RI

genes, indicating that they may be correlated with the panicle formation.

3.4.3. PPR genes

SC

Pentatricopeptide repeat (PPR) protein is one of the largest families of proteins which are involved in each step of organelles gene expression in plant (Binder et al.,

NU

2011; Dahan and Mireau, 2013; Shikanai and Fujii, 2013). In Arabidopsis, two PPR genes GPR23 and PPR4 have already been linked to the embryo development (Ding

MA

et al., 2006; Kocabek et al., 2006). In rice, a PPR protein sped1-D may take part in the regulation of peduncle and secondary branch elongation by blocking the function of

D

GID1L2, APO2/RFL and WOX3 (Jiang et al., 2014). Even though the expression of

PT E

sped1-D there was not significantly different in NIL-GW8 and NIL-gw8Amol between different stages, but APO2/RFL and WOX3 were significant down-regulated in the two rice lines from stages 2 to 3. The APO2/RFL was the highest down-regulated

CE

gene in NIL-GW8 with -7.677 during the stage transition from 2 to 3. In this study, seventeen PPR genes were identified as DEGs with log2FC >2.0 between different

AC

stages in the two NILs (Table 5). Some of them were up-regulated in NIL-GW8 and down-regulated in NIL-gw8Amol, like Os01g0589900 and Os07g0244400. Some were down-regulated in NIL-GW8 and up-regulated in NIL-gw8Amol, like Os02g0127600. Besides that, Os08g0402600 was up-regulated from stages 1 to 2, but down-regulated from stages 2 to 3 in NIL-gw8Amol. Moreover, we found that some PPR genes were identified as DEGs only in NIL-GW8 or NIL-gw8Amol, like Os02g0759500, Os02g0750400, Os01g0815900, Os07g0635800, Os02g0170000 and Os07g0260000. Although the physiological function of these PPR genes mentioned above was still unknown, we hypothesized that they may play a role in the panicle cell proliferation 15

ACCEPTED MANUSCRIPT and differentiation in development.

3.5. Gene Ontology (GO) enrichment analysis of DEGs Next, we conducted GO enrichment analysis to classify the functions of the DEGs in the process of panicle development. To explore the relevant biological functions of

PT

the DEGs, the significantly enriched GO terms of DEGs were characterized. The identified DEGs could be categorized into three main GO categories (biological

RI

process, cellular component, and molecular function) which contain 53 functional groups, a range of 1 to 3,428 DEGs in each group includes (Table S2). In NIL-GW8,

SC

the DEGs were significantly enriched in 44, 50 and 48 GO groups in W1-W2, W2-W3 and W1-W3, respectively (Fig. S1). In the three comparison groups above,

NU

cellular process (GO: 0009987), cell part (GO: 0044464) and DNA binding (GO: 0005488) were the most enriched GO groups in biological process, cellular

MA

component, and molecular function, respectively. In addition, we found a high proportion of genes from functional groups of single-organism process (GO:

D

0044699), cell (GO: 005623), catalytic activity (GO: 0003824), metabolic process

PT E

(GO: 0008152), organelle (GO: 0043226), catalytic activity (GO: 0003824), response to stimulus (GO: 0050896) and biological regulation (GO: 0065007). Interestingly,

Fig. S2).

CE

the most enriched GO groups was very similar to NIL-GW8 in NIL-gw8Amol (Table S2;

Moreover, to foster learning about the molecular mechanisms underlying the

AC

formation of spikelet in rice panicle, we put emphasis on GO functional classification analysis for up-regulated and down-regulated DEGs, by comparing the NIL-GW8 and NIL-gw8Amol at different stages (Table S3, Fig. 6). In comparison W1-S1 (stage 1), the 726 DEGs were allocated to 46 GO groups, 460 up-regulated genes and 266 down-regulated genes were classified into 46 and 42 groups, respectively. Among them, cell (GO: 0005623), cell part (GO: 0044464), organelle (GO: 0043226), cellular process (GO: 0042995), single-organism process (GO: 0044699), metabolic process (GO: 0008152) and response to stimulus (GO: 0050896) were the most enriched GO groups. As expected, the most enriched GO groups of W1-S1 (stage 1), W2-S2 (stage 16

ACCEPTED MANUSCRIPT 2) and W3-S3 (stage 3) were similar and all contained groups related to cellular activity. The GO groups “nutrient reservoir activity” (GO: 0045735), “extracellular matrix component” (GO: 0044420), “extracellular matrix” (GO: 0031012), “metallochaperone activity” (GO: 0016530) and “nucleoid” (GO: 0009295) were identified with a few DEGs at three stages. It might suggest that the genes relevant to these groups may not be the main factors which lead to the grain shape

PT

differences between NIL-GW8 and NIL-gw8Amol.

RI

To explore the function of GW8 and other DEGs related panicle development. We found GW8 and these genes were simultaneously assigned to single-organism process

SC

(GO: 0044699), developmental process (GO: 0032502), biological regulation (GO: 0065007), multicellular organismal process (GO: 0032501), reproductive process (GO:

NU

0022414), reproduction (GO: 0000003), DNA binding (GO: 0005488), organelle part (GO: 0044422), cell (GO: 0005623), cell part (GO: 0044464), organelle (GO:

MA

0043226).

D

3.6. Metabolic pathway by KEGG analysis of DEGs

PT E

In order to further understand the biological functions of the DEGs and summarize the biological pathways associated panicle development, the DEGs in each comparison were mapped to KEGG database (http://www.genome.ad.jp/kegg/). The

CE

DEGs significantly enriched pathways involved in metabolic or signal transduction were identified. Especially, the analysis of KEGG pathway were performed for the

AC

DEGs in different comparisons between NIL-GW8 and NIL-gw8Amol, allowing a rapid visualization and an overview over metabolism alterations about the biological pathways correlated with the formation of rice panicle. At stage 1, a total of 70 significant DEGs with pathway annotation were assigned to 57 KEGG pathways. Among the mapped pathways, nine pathways were significantly enriched (Pvalue < 0.05) in the comparison of W1-S1. The most enriched metabolic pathways were observed in taurine and hypotaurine metabolism (three genes, 4.29%), glycolysis / gluconeogenesis (nine genes, 12.86%), phenylpropanoid biosynthesis (six genes, 8.57%). The gene Os04g0404800 17

ACCEPTED MANUSCRIPT (AMP-binding enzyme) was detected in four pathways (Table S4; Fig. 7a). At stage 2 (W2-S2), 176 significant DEGs with pathway annotation were assigned to 75 KEGG pathways and 13 metabolism pathways were significantly enriched (Pvalue < 0.05). The carotenoid biosynthesis, plant hormone signal transduction, phenylpropanoid biosynthesis, phenylalanine metabolism, taurine and hypotaurine

PT

metabolism were identified (Pvalue < 0.001, Qvalue < 0.05) as the top five pathways. We calculated the gene counts in each pathway, there were 14, 27, 14, 10 and 5 DEGs

RI

in the five pathways, respectively. Compared with stage 1, the analysis showed taurine and hypotaurine metabolism and phenylpropanoid biosynthesis were

SC

significantly identified at stage 2 (Table S4; Fig. 7b).

At stage 3, 295 significant DEGs with pathway annotation were assigned to 100

NU

KEGG pathways. As a result, there were also a large amount of annotated DEGs were significantly assigned to thirty pathways (Pvalue < 0.05). Compared with W1-S1 and

MA

W2-S2, we found that more DEGs were detected in each pathway in W3-S3. Specific enrichment of genes was observed for pathways involved in phenylpropanoid

D

biosynthesis (25 genes, 25.47%), phenylalanine metabolism (16 genes, 5.42%) and

PT E

plant hormone signal transduction (27 genes, 9.15%) which have been identified as significant enrichment pathways in W1-S1 and W2-S2 as well (Table S4; Fig. 7c). The differences in metabolic and regulatory pathways indicate the different

AC

NIL-gw8Amol.

CE

molecular mechanisms underlying the panicle development between NIL-GW8 and

3.7. Analysis of DEGs identified as putative transcription factors Transcription factors play an important role in the process of plant gene expression and regulation by binding to the specific cis-acting element of target gene promoter region, as like GW8 which belongs to the SBP family. Given that transcription factors might have a major effect on the network of GW8, one of the objectives of our work was to identify transcription factors. We paid more attention to DEGs belonging to transcription factor families and summarized the family, gene description and GO annotation of these transcription factor genes in the comparisons 18

ACCEPTED MANUSCRIPT of W1-S1, W2-S2 and W3-S3 (Table S5; Fig. 7). The transcription factor database (http://drtf.cbi.pku.edu.cn/) was used to investigate transcription factors associated with panicle development. Together, a total of 345 DEGs encode known or putative transcription factors. As shown in Fig. 8, 46 putative transcription factors belonging to 15 transcription factor

PT

families were identified with significantly different expression levels in the W1-S1 comparison. Among of them, thirteen genes were characterized as AP2/ERF

RI

(APETALA2/ethylene responsive factor) gene family, constituting a large proportion of the identified transcription factors at stage 1. Based on FPKM values, twenty-five

SC

putative transcription factor genes were down-regulated and twenty-one were up-regulated.

NU

In the W2-S2 comparison (stage 2), 129 putative transcription factor genes belonging to 31 transcription factor families (e.g., CAMTA, WRKY, NAC and MYB)

MA

were differentially expressed between the panicles of the two rice lines. These genes included relative large number of AP2/ERF, bHLH (basic helix-loop-helix) and

D

WRKY transcription factors, suggesting that they may play important roles in the

PT E

genes network with GW8 correlation.

In the W3-S3 comparison (stage 3), 170 putative transcription factor genes belonging to 32 transcription factor families were significantly differentially

CE

expressed. Compared with stages 1 and 2, there were no DEGs were identified as the members of FAR1, GRF, MIKC, SRS, Trihelix and WOX transcription factor

AC

families. The ERF, bHLH and NAC families have the most genes which contained 22, 21 and 17 putative transcription factor DEGs, respectively. We hypothesized that these genes might be important regulators contribute to the structure difference of NIL-GW8 and NIL-gw8Amol panicles which induced by GW8 and gw8Amol. These

transcription

factor

genes

were

significantly

up-regulated

or

down-regulated, but the trend of expression of them varied at the three stages. For example, the expression of ZF-HD was only detected in W3-S3, and MIKC was only detected in W1-S1 and W2-S2. In short, a lot of transcription factor genes exhibited dynamic expression profiles during the panicle development of NIL-GW8 and 19

ACCEPTED MANUSCRIPT NIL-gw8Amol. Further analysis is needed to explore the involvement of these transcription factor families in the regulation of rice panicle formation.

3.8. Validation of RNA-Seq technology To validate the RNA-Seq results, ten genes were randomly selected for analysis

PT

by qRT-PCR referred to the same RNA extracts as for RNA-seq experiments (Fig. 9). The expression profiles of the ten genes determined by qRT-PCR were consistent

RI

with the RNA-Seq results, confirming a good correspondence between qRT-PCR and

SC

RNA-seq.

4. Discussion

NU

4.1. Biosynthesis of secondary metabolites may contribute to the formation of rice panicle

MA

According to the KEGG pathway analyses of DEGs among the three panicle development stages in NIL-GW8 and NIL-gw8Amol, we found that many DEGs were

D

significantly enriched in the pathways related to biosynthesis of secondary

PT E

metabolites, such as phenylpropanoid biosynthesis (ko00940), phenylalanine metabolism (ko00360), stilbenoid, diarylheptanoid and gingerol biosynthesis (ko00945) and flavonoid biosynthesis (ko00941). Among of them, the pathway of

CE

flavonoid biosynthesis is a branch of phenylpropanoid pathway which is associated with phenylalanine metabolism (Winkel-Shirley, 2001; Buer et al., 2010). We

AC

speculate that these secondary metabolites may contribute to the formation of rice panicle.

We all know that plant secondary metabolites constitute an unusually large group of structurally diversified compounds that can be produced in planta from various primary metabolites or their biosynthetic intermediates, constitutively or in response to different environmental stimuli (Piasecka et al., 2015), and were not generally regarded as indispensable components for the normal life activities. Decades of studies have suggested that the secondary metabolites were necessary in maintaining normal metabolism and completing life cycle of the plants and play an important role 20

ACCEPTED MANUSCRIPT in the life activities (Kutchan, 2001). For instance, gibberellin (GA3) and indole acetic acid (IAA) are involved in regulating many activities in plant cells; terpenes (chlorophyll and carotenoid) are involved in the processes of photosynthesis; flavonoids participates in auxin transport by regulating the activity of phosphatase and protein kinase in shoot and root apical meristem; carotenoids content influenced the

PT

drought tolerance in dendrobium moniliforme (Peer and Murphy, 2007; Pattanayak and Tripathy, 2016; Wu et al., 2016).

RI

The phenotypic difference of NIL-GW8 and NIL-gw8Amol is mainly because of the different gene function of GW8 and gw8Amol which belong to SPL (SBP-like)

SC

transcription factor family. In recent years, research found that many transcription factors participated in the regulation of plant secondary metabolites biosynthesis and

NU

accumulation (Appelhagen et al., 2011). Two Petunias transcription factor genes an2 and an11 have been shown to be effective in regulating gene expression in the

MA

flavonoids biosynthetic pathway (Quattrocchio et al., 1998). Overexpression of AtSPL9 and AtSPL10 remarkably negatively regulated the synthesis of anthocyanin,

D

and it show a negative relationship between the activity of miR156 and the contents of

PT E

flavonoids in Arabidopsis (Gou et al., 2011). Based on these findings, we focused on the expression of the genes related to secondary metabolites in our research. Just like what we thought, we found some

CE

genes were significantly enriched in the biosynthetic pathway of secondary metabolites. For example, Os04g0310900 (4-coumarate-CoA ligase) was enriched in

AC

phenylalanine metabolism pathway (Table 2). Despite the functions of secondary metabolites in the panicle are not entirely clear, these studies confirm that they do play important roles in the panicle formation. In summary, we speculate that GW8 regulate the growth and development of rice panicle by influencing the biosynthesis of secondary metabolites.

4.2. Genes related to plant hormone signal transduction are important for the formation of rice panicle In higher plants, the generation and establishment of meristem identity of 21

ACCEPTED MANUSCRIPT axillary meristem have been regard as the key factors of plant morphological and inflorescence architecture (Oikawa and Kyozuka, 2009). And he growth and development of axillary meristem are regulated by plant hormones (McSteen, 2009). In this study, a lot of DEGs in panicles between NIL-GW8 and NIL-gw8Amol were categorized into the pathway of plant hormone signal transduction, which was

PT

represented by many hormones-related genes, especially at stages 2 and 3, indicating that these genes are likely to be involved in the regulation of panicle development in

RI

rice by participating directly or indirectly in the plant hormone signal transduction. As a signaling molecule, auxin is one of most important hormones involved in a

SC

variety of development processes in plants. Studies in recent years have shown that auxin is involved in the modulating development of rice panicles and spikelets

NU

through direct induction of other growth regulating genes (Gao et al., 2016). And other research shows that auxin is necessary for the rice axillary meristem initiation

MA

and the second branch differentiation (Benjamins and Scheres, 2008; Mockaitis and Estelle, 2008; McSteen, 2009). Also, auxin is also involved in regulating the

D

inflorescence meristem development (Zhao, 2008). During the rice reproductive

PT E

process, the content of auxin in lateral meristem determined the number of primary branches and secondary branches (Komatsu et al., 2003). ASP1, a regulatory factor of auxin, controls the development of panicles in rice, and smaller number of branches

CE

and abnormal panicles were found in the mutant asp1 (Yoshida et al., 2012). In this research, many genes related to the auxin (IAA) synthesis (OsJRA2 and LC1), auxin

AC

transport (OsAUX1) and auxin response (OsSAURs) were significantly differentially expressed between the NIL-GW8 and NIL-gw8Amol panicles. It indicated that these genes might be involved in panicle developmental events. During the study many cytokinin-signalling and cytokinin-responsive genes were among the transcripts that were differentially expressed during panicle development of the two rice lines at different stages, such as a series of A-type and B-type response regulators (OsRR1, OsRR6, OsRR11, ORR1 and ORR4). Some theoretical and experimental evidence showed that cytokinin is tightly correlated with the panicle length and the number of branches and spikelets per panicle in rice (Ashikari et al., 22

ACCEPTED MANUSCRIPT 2005; Kurakawa et al., 2007; Li et al., 2013). The QTL GNP1, encoding GA20ox1, increases grain number and yield by increasing cytokinin activity in rice panicle meristems (Wu et al., 2016). Overexpression of OsRR6, a negative regulator of cytokinin signal transduction pathway, has resulted in insufficient development of panicles in rice (Hirose et al., 2007). In addition, a few other studies also

PT

demonstrated that A-type and B-type response regulator families have a vital role in the cytokinin signaling pathway (Zhao et al., 2015; Cho et al., 2016). Accordingly, we

RI

hypothesized that the expression variations of OsRRs and ORRs between NIL-GW8 and NIL-gw8Amol panicles may affect the signal transduction of cytokinin, resulting in

SC

the phenotypic differences of panicles.

In Arabidopsis, plant hormone jasmonic acid (JA) participates in the signal

NU

transduction pathway through activating the expression of JA-responsive genes by degrading the JAZ (jasmonate ZIM-domain) repressors (Chini et al., 2007). A

MA

previous study on spikelet in rice showed that JA was associated with the determination of rice spikelet morphogenesis, and the OsJAZ1 mutant exhibited

D

altered spikelet morphology with changed floral organ identity and number (Cai et al.,

PT E

2014). Decreased expression of OsCOI1 (Coronatine-insensitive 1), the receptor of JA, in rice RNAi lines resulted in long panicles compared with the wild type (Yang et al., 2012). In this study, we identified nine JAZ genes (OsJAZ3, OsJAZ4, OsJAZ6,

CE

OsJAZ7, OsJAZ9, OsJAZ10, OsJAZ11, OsJAZ12 and OsJAZ13) and a JAZ-interacting transcription factor (OsMYC2) in plant hormone signal transduction pathway.

AC

In addition to the above genes, according to the analysis of KEGG pathways, several genes mainly related to other hormone signal transduction pathways of plants were detected between the two rice lines, such as ABA-signalling genes (OsSAPK1, OsSAPK2 and OsSAPK9), ABA responsive element binding factor (OsABF1), ethylene-insensitive gene (OsEIL1) and ethylene receptor (ETR2). So, the expressions of these genes are likely to be affected by GW8 and play integral roles in the process of growth and development of panicles. In summary, our results provide important information for future functional studies of plant hormone signal transduction-related genes and a research foundation 23

ACCEPTED MANUSCRIPT for the mechanisms involved in rice panicle formation.

4.3. Difference in the regulation pattern of the formation of rice panicle between GW8 and GW7 Our previous research found that GW7 is a major QTL for both grain length and

PT

width in rice, increase of GW7 expression, longitudinal cell division of the grain was promoted and the transverse cell division was reduced, as a result of the elongation of

RI

grain (Wang et al., 2015). GW8 plays an important role in the development of panicle in rice, GW8 bound directly to the GW7 promoter and repressed its expression (Wang

SC

et al., 2012; Wang et al., 2015). But for this study, we found that there were no statistically significant differences in the expression levels of GW7 between NIL-GW8

NU

and NIL-gw8Amol at the three stages. So, we speculate that GW7 regulates the development of panicle in the latter stage which after the period of floral organs

MA

differentiation (In7). We compared the transcript levels of GW7 in more different panicle differentiation stages of rice. The result was as we’d expect, qRT-PCR

D

analysis identified differences in GW7 transcript abundance between NIL-GW8 and

PT E

NIL-gw8Amol, and the gene being transcribed more strongly in NIL-gw8Amol than in NIL-GW8 after the stage In7 (Fig. S3). Moreover, the expressions of GW7-related genes were analyzed among the

CE

different developmental stages in NIL-GW8 and NIL-gw8Amol. In the previous study, we have identified 18 candidate genes which interact with GW7 by yeast two-hybrid

AC

assay with GW7 as bait (Wang et al., 2015). For these 18 GW7-related genes, only three genes (Os01g0153800, Os02g0818000 and Os03g0135700) showed differences in expression between different stages in NIL-GW8 and NIL-gw8Amol with log2FC >2.0, but no significant difference between NIL-GW8 and NIL-gw8Amol at different stages in this study. Conclusion according to analysis on the expression differences of GW7 and its related genes, GW7 and GW7-related genes might not be the key regulation factors of spikelet formation as lateral branch of the inflorescence meristem. Our previous research has also shown that there are no differences in the number of primary and 24

ACCEPTED MANUSCRIPT secondary branches and spikelets per panicle between NIL-GW7TFA and NIL-gw7HJX74, but there exists a significant difference between NIL-GW8 and NIL-gw8Amol (Wang et al., 2012; Wang et al., 2015). Therefore, it would be assumed that the regulation pattern of GW7 different from GW8 in different stages of rice panicle development, although GW8 bound to the GW7 promoter and regulate its expression and they all

PT

regulate rice grain development.

RI

5. Conclusion

Our transcriptome analysis provided a comprehensive overview of the

SC

transcriptome changes that specifically occur in panicles of two NILs (NIL-GW8 and NIL-gw8Amol), which were characterized by distinct phenotypic differences of panicles.

NU

Based on the analysis of GO and KEGG, we found that gene expression profiles were dynamic during rice panicle development. Various types of stage-specifically

MA

expressed genes were identified from the three panicle developmental stages between NIL-GW8 and NIL-gw8Amol, and the identified DEGs were involved in diverse

D

pathways, including biosynthesis of secondary metabolites and plant hormone signal

PT E

transduction. The identification of DEGs including some MADS-box genes, cytochrome P450 (CYP) genes, pentatricopeptide repeat (PPR) protein genes and various transcription factors will enhance our understanding of the mechanism of

AC

CE

panicle formation affected by GW8.

25

ACCEPTED MANUSCRIPT Conflict of interests All authors declare that they have no conflict of interest.

Authors’ contributions S.K. performed the bioinformatics analysis and wrote the manuscript. X.J.L.

PT

conducted the sample collection and library construction. X.L. performed the RNA extraction. W.Y. performed the qRT-PCR. H.Z. participated in study coordination and

RI

project, and provided essential suggestions for this work. G.L. supervised and complemented the writing. G.Z. and S.W. conceived the study, participated in its

SC

design and coordination, and corrected the manuscript. All authors read and approved

NU

the final manuscript.

Acknowledgments

MA

This work was supported by the National Natural Science Foundation of China (91635302) and Foundation for Young Talents in Higher Education of Guangdong,

AC

CE

PT E

D

China (2014KQNCX037).

26

ACCEPTED MANUSCRIPT Figure legends Fig. 1. The phenotypic differences of the grain formed by NIL-GW8 and NIL-gw8Amol. Appearance of the grains (a). Grain length (b). Grain width (c). Bar = 5mm. Each column represents the mean ± s.e.m. (n = 120). Asterisks represent significant

PT

difference determined by Student's t-test at p-value < 0.01 (**).

Fig. 2. Length distribution of genes characterized from RNA-seq libraries in rice

SC

RI

developing panicles.

Fig. 3. The distribution of gene expression at three different panicle development

NU

stage in NIL-GW8 and NIL-gw8Amol. The twelve rice chromosomes are displayed by the first outer ring. Six concentric rings from outside to inside represent the

MA

corresponding six libraries (W1, W2, W3, S1, S2, and S3) in this figure, respectively. Colors represent (from green to red) the expression levels of all detected transcripts in

PT E

fragments mapped) value.

D

every library based on the FPKM (Fragments Per Kilobase of transcript per Million

Fig. 4. Venn diagram showing the genes expressed in each of six samples of of rice

CE

panicles development in NIL-GW8 and NIL-gw8Amol. (1) W1: stage 1 (In5) in NIL-GW8; (2) W2: stage 2 (In6) in NIL-GW8; (3) W3: stage 3 (In7) in NIL-GW8; (4) S1: stage 1 (In5) in NIL-gw8Amol; (5) S2: stage 2 (In6) in NIL-gw8Amol; (6) S3: stage 3

AC

(In7) in NIL-gw8Amol.

Fig. 5. Analysis of DEGs in each comparison among the different panicles development stages in NIL-GW8 and NIL-gw8Amol. The number of up-regulated and down-regulated genes in different comparisons is summarized.

Fig. 6. Functional analysis of differentially expressed genes (DEGs) based on RNA-Seq data in different comparisons. W1-S1 (a), W2-S2 (b), W3-S3 (c). 27

ACCEPTED MANUSCRIPT

Fig. 7. KEGG enrichment of DEGs at different panicle stages between NIL-GW8 and NIL-gw8Amol. The scatterplot of enriched GO terms between W1 and S1 (a); the scatterplot of enriched GO terms between W2 and S2 (b); the scatterplot of enriched

PT

GO terms between W3 and S3 (c).

Fig. 8. The number of DEGs belonging to different transcription factor families

RI

detected in the W1-S1, W2-S2, and W3-S3 comparisons. In total, we identified 345

SC

putative TF DEGs, which are classified into 38 families.

Fig. 9. qRT-PCR validation of the relative expression data of ten genes obtained in

NU

RNA-seq analysis. y axes represent the relative expression levels at different stages in NIL-GW8 and NIL-gw8Amol which are represented relative to W1 (assigned as a value

MA

of 1). OsActin1 was used as an internal control gene for transcript normalization.

PT E

Table legends

D

Error bars show standard deviations for three independent biological assays.

CE

Table 1 Summary of reads mapping to rice genome

AC

Table 2 Statistics of the main up-expressed and down-expressed genes

Table 3 Statistics of significant differently expressed MADS-box genes

Table 4 Statistics of significant differently expressed Cytochrome P450 genes

Table 5 Statistics of significant differently expressed PPR genes

28

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

CE

Fig. 1. The phenotypic differences of the grain formed by NIL-GW8 and NIL-gw8Amol. Appearance of the grains (a). Grain length (b). Grain width (c). Bar = 5mm. Each

AC

column represents the mean ± s.e.m. (n = 120). Asterisks represent significant difference determined by Student's t-test at p-value < 0.01 (**).

29

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 2. Length distribution of genes characterized from RNA-seq libraries in rice

AC

CE

PT E

D

developing panicles.

30

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 3. The distribution of gene expression at three different panicle development

CE

stage in NIL-GW8 and NIL-gw8Amol. The twelve rice chromosomes are displayed by the first outer ring. Six concentric rings from outside to inside represent the

AC

corresponding six libraries (W1, W2, W3, S1, S2, and S3) in this figure, respectively. Colors represent (from green to red) the expression levels of all detected transcripts in every library based on the FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) value.

31

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

MA

Fig. 4. Venn diagram showing the genes expressed in each of six samples of of rice panicles development in NIL-GW8 and NIL-gw8Amol. (1) W1: stage 1 (In5) in NIL-GW8; (2) W2: stage 2 (In6) in NIL-GW8; (3) W3: stage 3 (In7) in NIL-GW8; (4)

AC

CE

PT E

(In7) in NIL-gw8Amol.

D

S1: stage 1 (In5) in NIL-gw8Amol; (5) S2: stage 2 (In6) in NIL-gw8Amol; (6) S3: stage 3

32

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

CE

Fig. 5. Analysis of DEGs in each comparison among the different panicles development stages in NIL-GW8 and NIL-gw8Amol. The number of up-regulated and

AC

down-regulated genes in different comparisons is summarized.

33

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 6. Functional analysis of differentially expressed genes (DEGs) based on RNA-Seq data in different comparisons. W1-S1 (a), W2-S2 (b), W3-S3 (c).

34

ACCEPTED MANUSCRIPT

T P

I R

C S U

N A

D E

M

T P E

C C

Fig. 7. KEGG enrichment of DEGs at different panicle stages between NIL-GW8 and NIL-gw8Amol. The scatterplot of enriched GO terms

A

between W1 and S1 (a); the scatterplot of enriched GO terms between W2 and S2 (b); the scatterplot of enriched GO terms between W3 and S3 (c).

35

ACCEPTED MANUSCRIPT

T P

I R

C S U

N A

D E

M

T P E

C C

Fig. 8. The number of DEGs belonging to different transcription factor families detected in the W1-S1, W2-S2, and W3-S3 comparisons. In total,

A

we identified 345 putative TF DEGs, which are classified into 38 families.

36

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 9. qRT-PCR validation of the relative expression data of ten genes obtained in RNA-seq analysis. y axes represent the relative expression levels at different stages in NIL-GW8 and NIL-gw8Amol which are represented relative to W1 (assigned as a value of 1). OsActin1 was used as an internal control gene for transcript normalization. Error bars show standard deviations for three independent biological assays.

37

ACCEPTED MANUSCRIPT Table 1 Summary of reads mapping to rice genome. Total Reads

Mapped Reads

Mapped Ratio

Unique Mapped Reads

Unique Mapped Ratio

W1

46,479,240

37,345,492

80.35%

36,297,170

78.09%

W2

53,148,798

42,452,813

79.88%

41,216,953

77.55%

W3

54,120,416

43,417,595

80.22%

42,234,568

78.04%

S1

40,492,882

33,347,884

82.35%

32,398,176

80.01%

S2

40,822,488

33,533,503

82.14%

32,345,222

79.23%

S3

46,103,974

37,048,147

80.36%

35,813,220

77.68%

AC

CE

PT E

D

MA

NU

SC

RI

PT

Sample

38

ACCEPTED MANUSCRIPT Table 2 Statistics of the main up-expressed and down-expressed genes. Annotation

Regulated

NewGene_446

Hypothetical protein

NewGene_3734

Unknown protein

Os08g0527600

log2 FC W2 - S2

W3 - S3

up

9.17

9.10

6.04

up

3.78

5.24

7.53

Unnamed protein product

up

6.37

6.37

7.38

Os12g0210601

Hypothetical protein

up

7.23

7.36

7.12

Os08g0528100

Ubiquitin carboxyl-terminal

up

5.78

5.81

7.10

6.91

6.35

5.91

5.14

8.56

-

7.92

6.57

-

6.61

6.21

5.69

-5.43

-4.19

6.33

up

4.56

3.83

5.72

up/down

-6.78

-

5.62

up/down

6.13

-3.93

4.83

up/down

5.85

-6.17

-5.37

up

5.74

5.62

-

PT

W1 - S1

NU

Gene ID

Hypothetical protein

up

Os05g0469600

Pyruvate decarboxylase

up

Os03g0725200

Formin-like protein

up

Os10g0366400

Unknown protein

up

Os02g0127600

putative pentatricopeptide (PPR) repeat-containing protein hypothetical protein

Os02g0260500

Uncharacterized protein

Os01g0159000

Proline-rich receptor

Os01g0589900

putative pentatricopeptide (PPR)

up/down

MA

NewGene_3514

repeat-containing protein

SC

NewGene_3689

RI

hydrolase

Pectinesterase

Os03g0172100

Protein CHUP1, chloroplastic

up

5.59

5.44

Os05g0469701

Pyruvate decarboxylase

up

5.25

7.07

Os04g0310900

4-coumarate-CoA ligase

up

5.22

-

3.79

Os05g0162800

Transcription factor

up

-

6.45

4.87

PT E

D

Os02g0288100

Os06g0100550

Unknown protein

up

4.35

4.56

4.83

Hypothetical protein

down

-8.08

-6.85

-6.32

Hypothetical protein

down

-7.14

-5.19

-5.69

Unknown protein

down

-6.52

-6.21

-

Hypothetical protein

down

-6.30

-7.65

-5.95

Os09g0513900

Unnamed protein product

down

-5.95

-5.83

-4.38

Os12g0234000

Unknown protein

down

-5.91

-5.47

-2.05

Os01g0290800

ABC transporter B family

up/down

4.85

-5.72

1.55

down

-

4.91

-

Os10g0370100 Os07g0646200

AC

NewGene_3441

CE

NewGene_448

member NewGene_3073

Hypothetical protein

Os06g0361500

Unknown protein

down

-

-4.80

-4.30

Os03g0226200

Non-symbiotic hemoglobin

up/down

-

3.06

-7.50

Os12g0446500

metallopeptidase family M24

up/down

-

6.01

-7.09

down

-5.89

-

-5.56

down

-1.89

-2.30

-5.47

containing protein Os03g0229500

MATE efflux family protein

Os02g0586900

Glycine-rich cell wall structural protein (Precursor)

39

ACCEPTED MANUSCRIPT Table 3 Statistics of significant differently expressed MADS-box genes.

7.386 7.604 7.420 7.436 4.537 4.247 4.915 3.341 4.215 3.575 2.621 2.269 -3.252 -4.177 -2.024 -2.977 -3.903 -3.734 -7.073

PT

6.924 6.669 6.663 6.662 3.932 3.580 3.557 3.396 2.837 2.779 2.742 2.112 -2.467 -2.959 -3.153 -3.612 -4.549 -4.640 -5.823

AC

CE

PT E

D

S2-S3

RI

OsMADS13 OSMADS3 OsMADS7 OsMADS8 OsMADS17 OsMADS21 OsMADS1 OsMADS29 OsMADS58 OsMADS6 OsMADS4 OsMADS56 OsMADS26 OsMADS22 OsMDP1 OsMADS34 OsMADS50 OsMADS20 OsMADS32 OsMADS55

W2-W3

SC

Os12g0207000 Os01g0201700 Os08g0531700 Os09g0507200 Os04g0580700 Os01g0886200 Os03g0215400 Os02g0170300 Os05g0203800 Os02g0682200 Os05g0423400 Os10g0536100 Os08g0112700 Os02g0761000 Os03g0186600 Os03g0753100 Os03g0122600 Os12g0501700 Os01g0726400 Os06g0217300

log2 FC

NU

Gene Name

MA

Gene ID

40

ACCEPTED MANUSCRIPT Table 4 Statistics of significant differently expressed Cytochrome P450 genes. log2 FC

Gene Name

Os10g0525000 Os01g0804400 Os03g0332100 Os02g0186900 Os04g0570500 Os05g0482400 Os04g0573900 Os09g0441400 Os04g0407900 Os07g0635300 Os01g0227800 Os07g0603700 Os11g0289700 Os10g0486100 Os10g0403000

CYP714B2

S2-S3

6.510 6.179 5.856 5.782 5.777 4.952 4.704 4.548 4.035 4.022 3.993 3.935 3.003 -7.630

5.025 5.241 4.163 6.060 6.485 6.312 5.416 5.382 4.098 3.221 4.943 5.655 4.135 -6.811

SC

RI

EUI1

W2-W3

PT

Gene ID

NU

CYP78A13

AC

CE

PT E

D

MA

PLA1

41

ACCEPTED MANUSCRIPT Table 5 Statistics of significant differently expressed PPR genes. log2 FC

Gene ID

W1-W2

Os01g0589900 Os08g0402600 Os07g0244400 Os03g0363700 Os12g0552300 Os02g0759500 Os02g0750400 Os08g0340900 Os10g0116000 Os11g0433101 Os01g0815900 Os07g0635800 Os07g0621100 Os02g0127600 Os03g0283500 Os02g0170000 Os07g0260000

W2-W3

S1-S2

6.407

-6.071 4.894

-4.247 -2.258

PT

4.870 2.101 -4.793

S2-S3

RI

4.511 -4.545

NU

SC

4.740 4.144

MA

3.935 -3.599 -3.172

4.098 3.745 -3.604 -3.810 6.446 2.984

AC

CE

PT E

D

-3.148

42

ACCEPTED MANUSCRIPT References Appelhagen, I., Lu, G., Huep, G., Schmelzer, E., Weisshaar, B., Sagasser, M., 2011. TRANSPARENT TESTA1 interacts with R2R3-MYB factors and affects early and late steps of flavonoid biosynthesis in the endothelium of Arabidopsis thaliana seeds. Plant J. 67 (3), 406-419. https://doi.org/10.1111/j.1365-313X.2011.04603.x. Ashikari, M., Matsuoka, M., 2006. Identification, isolation and pyramiding of quantitative trait loci for rice breeding. Trends Plant Sci. 11 (7), 344-350. https://doi.org/10.1016/j.tplants.2006.05.008. Ashikari, M., Sakakibara, H., Lin, S.Y., Yamamoto, T., Takashi, T., Nishimura, A., Angeles, E.R., Qian, Q., Kitano, H., Matsuoka, M., 2005. Cytokinin oxidase regulates rice grain production.

PT

Science 309 (5735), 741-745. https://doi.org/10.1126/science.1113373.

Benjamins, R., Scheres, B., 2008. Auxin: The looping star in plant development. Annu. Rev. Plant Biol.

RI

59 (59), 443-465. https://doi.org/10.1146/annurev.arplant.58.032806.103805.

Binder, S., Holzle, A., Jonietz, C., 2011. RNA Processing and RNA Stability in Plant Mitochondria.

SC

Adv. Plant Biochem. Mol. Biol. 1, 107-130. https://doi.org/10.1007/978-0-387-89781-3_5. Bolwell, G.P., Bozak, K., Zimmerlin, A., 1994. Plant cytochrome P450. Phytochemistry 37, 1491-1506. https://doi.org/10.1016/S0031-9422(00)89567-9.

NU

Buer, C.S., Imin, N., Djordjevic, M.A., 2010. Flavonoids: New Roles for Old Molecules. J. Integr. Plant Biol. 52, 98-111. https://doi.org/10.1111/j.1744-7909.2010.00905.x. Cai, Q., Yuan, Z., Chen, M., Yin, C., Luo, Z., Zhao, X., Liang, W., Hu, J., Zhang, D., 2014. Jasmonic regulates

spikelet

development

in

rice.

Nat.

Commun.

5

(3),

319-333.

MA

acid

https://doi.org/10.1038/ncomms4476.

Chen, R., Shen, L., Wang, D., Wang, F., Zeng, H., Chen, Z., Peng, Y., Lin, Y., Tang, X., Deng, M., Yao, N., Luo, J., Xu, Z., Bai, S., 2015. A Gene Expression Profiling of Early Rice Stamen

D

Development that Reveals Inhibition of Photosynthetic Genes by OsMADS58. Mol. Plant. 8 (7), 1069-1089. https://doi.org/10.1016/j.molp.2015.02.004.

PT E

Chini, A., Fonseca, S., Fernandez, G., Adie, B., Chico, J.M., Lorenzo, O., Garcia-Casado, G., Lopez-Vidriero, I., Lozano, F.M., Ponce, M.R., Micol, J.L., Solano, R., 2007. The JAZ family of repressors is the missing link in jasmonate signalling. Nature 448 (7154), 666-671. https://doi.org/10.1038/nature06006.

CE

Cho, L., Yoon, J., Pasriga, R., An, G., 2016. Homodimerization of Ehd1 Is Required to Induce Flowering in Rice. Plant Physiol. 170 (4), 2159-2171. https://doi.org/10.1104/pp.15.01723. Dahan, J., Mireau, H., 2013. The Rf and Rf-like PPR in higher plants, a fast-evolving subclass of PPR

AC

genes. RNA Biol. 10 (9), 1469-1476. https://doi.org/10.4161/rna.25568. Deshpande, G.M., Ramakrishna, K., Chongloi, G.L., Vijayraghavan, U., 2015. Functions for rice RFL in vegetative axillary meristem specification and outgrowth. J. Exp. Bot. 66 (9), 2773-2784. https://doi.org/10.1093/jxb/erv092. Devos, K.M., 2005. Updating the 'Crop circle'. Curr. Opin. Plant Biol. 8 (2), 155-162. https://doi.org/10.1016/j.pbi.2005.01.005. Ding, Y.H., Liu, N.Y., Tang, Z.S., Liu, J., Yang, W.C., 2006. Arabidopsis GLUTAMINE-RICH PROTEIN23 is essential for early embryogenesis and encodes a novel nuclear PPR motif protein that interacts with RNA polymerase II subunit III. Plant Cell 18 (4), 815-830. https://doi.org/10.1105/tpc.105.039495. Duan, K., Li, L., Hu, P., Xu, S., Xu, Z., Xue, H., 2006. A brassinolide-suppressed rice MADS-box transcription factor, OsMDP1, has a negative regulatory role in BR signaling. Plant J. 47 (4), 43

ACCEPTED MANUSCRIPT 519-531. https://doi.org/10.1111/j.1365-313X.2006.02804.x. Duan, P., Ni, S., Wang, J., Zhang, B., Xu, R., Wang, Y., Chen, H., Zhu, X., Li, Y., 2016. Regulation of OsGRF4 by OsmiR396 controls grain size and yield in rice. Nat. Plants 2 (1), 15203. http://dx.doi.org/10.1038/nplants.2015.203. Endo-Higashi, N., Izawa, T., 2011. Time Genes Heading date 1 and Early heading date 1 Together Control

Panicle

Development

in

Rice.

Plant

Cell

Physiol.

52

(6),

1083-1094.

https://doi.org/10.1093/pcp/pcr059. Ewing, B., Hillier, L., Wendl, M.C., Green, P., 1998. Base-calling of automated sequencer traces using phred. I. Accuracy assessment. Genome Res. 8 (3), 175-185. https://doi.org/10.1101/gr.8.3.175.

PT

Fan, C.H., Xing, Y.Z., Mao, H.L., Lu, T.T., Han, B., Xu, C.G., Li, X.H., Zhang, Q.F., 2006. GS3, a major QTL for grain length and weight and minor QTL for grain width and thickness in rice,

RI

encodes a putative transmembrane protein. Theor. Appl. Genet. 112 (6), 1164-1171. https://doi.org/10.1007/s00122-006-0218-1.

SC

Gao, F., Wang, K., Liu, Y., Chen, Y., Chen, P., Shi, Z., Luo, J., Jiang, D., Fan, F., Zhu, Y., Li, S., 2016. Blocking miR396 increases rice yield by shaping inflorescence architecture. 15196. https://doi.org/10.1038/nplants.2015.196.

Nat. Plants 2 (1),

NU

Gou, J., Felippes, F.F., Liu, C., Weigel, D., Wang, J., 2011. Negative Regulation of Anthocyanin Biosynthesis in Arabidopsis by a miR156-Targeted SPL Transcription Factor. Plant Cell 23 (4), 1512-1522. https://doi.org/10.1105/tpc.111.084525.

MA

Hirose, N., Makita, N., Kojima, M., Kamada-Nobusada, T., Sakakibara, H., 2007. Overexpression of a type-A response regulator alters rice morphology and cytokinin metabolism. Plant Cell Physiol. 48 (3), 523-539. https://doi.org/10.1093/pcp/pcm022.

Hu, Y., Liang, W., Yin, C., Yang, X., Ping, B., Li, A., Jia, R., Chen, M., Luo, Z., Cai, Q., Zhao, X.,

D

Zhang, D., Yuan, Z., 2015. Interactions of OsMADS1 with Floral Homeotic Genes in Rice Flower Development. Mol. Plant 8 (9), 1366-1384. https://doi.org/10.1016/j.molp.2015.04.009.

PT E

Huang, X., Qian, Q., Liu, Z., Sun, H., He, S., Luo, D., Xia, G., Chu, C., Li, J., Fu, X., 2009. Natural variation at the DEP1 locus enhances grain yield in rice. Nat. Genet. 41 (4), 494-497. https://doi.org/10.1038/ng.352.

Ikeda, K., Sunohara, H. and Nagato, Y., 2004. Developmental course of inflorescence and spikelet in

CE

rice. Breed. Sci. 54 (2), 147-156. https://doi.org/10.1270/jsbbs.54.147. Ikeda-Kawakatsu, K., Maekawa, M., Izawa, T., Itoh, J., Nagato, Y., 2012. ABERRANT PANICLE ORGANIZATION 2/RFL, the rice ortholog of Arabidopsis LEAFY, suppresses the transition from

AC

inflorescence meristem to floral meristem through interaction with APO1. Plant J. 69 (1), 168-180. https://doi.org/10.1111/j.1365-313X.2011.04781.x. Ishii, T., Numaguchi, K., Miura, K., Yoshida, K., Pham, T.T., Than, M.H., Yamasaki, M., Komeda, N., Matsumoto, T., Terauchi, R., Ishikawa, R., Ashikari, M., 2013. OsLG1 regulates a closed panicle trait in domesticated rice. Nat. Genet. 45 (4), 462-465. https://doi.org/10.1038/ng.2567. Ishimaru, K., Hirotsu, N., Madoka, Y., Murakami, N., Hara, N., Onodera, H., Kashiwagi, T., Ujiie, K., Shimizu, B., Onishi, A., Miyagawa, H., Katoh, E., 2013. Loss of function of the IAA-glucose hydrolase gene TGW6 enhances rice grain weight and increases yield. Nat. Genet. 45 (6), 707-711. https://doi.org/10.1038/ng.2612. Itoh, J., Nonomura, K., Ikeda, K., Yamaki, S., Inukai, Y., Yamagishi, H., Kitano H., Nagato Y., 2005. Rice plant development: from zygote to spikelet. Plant Cell Physiol. 46(1), 23-47. https://doi.org/10.1093/pcp/pci501. 44

ACCEPTED MANUSCRIPT Jiang, G., Xiang, Y., Zhao, J., Yin, D., Zhao, X., Zhu, L., Zhai, W., 2014. Regulation of Inflorescence Branch Development in Rice Through a Novel Pathway Involving the Pentatricopeptide Repeat Protein sped1-D. Genetics 197 (4), 1395-1407. https://doi.org/10.1534/genetics.114.163931. Jiao, Y., Wang, Y., Xue, D., Wang, J., Yan, M., Liu, G., Dong, G., Zeng, D., Lu, Z., Zhu, X., Qian, Q., Li, J., 2010. Regulation of OsSPL14 by OsmiR156 defines ideal plant architecture in rice. Nat. Genet. 42 (6), 541-544. https://doi.org/10.1038/ng.591. Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., Salzberg, S.L., 2013. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14, R36. https://doi.org/10.1186/gb-2013-14-4-r36.

PT

Kocabek, T., Repkova, J., Dudova, M., Hoyerova, K., Vrba, L., 2006. Isolation and characterization of a novel semi-lethal Arabidopsis thaliana mutant of gene for pentatricopeptide (PPR)

RI

repeat-containing protein. Genetica 128 (1), 395-407. https://doi.org/10.1007/s10709-006-7518-x. Komatsu, M., Chujo, A., Nagato, Y., Shimamoto, K., Kyozuka, J., 2003. FRIZZY PANICLE is required

SC

to prevent the formation of axillary meristems and to establish floral meristem identity in rice spikelets. Development 130 (16), 3841-3850. https://doi.org/10.1242/dev.00564. Kurakawa, T., Ueda, N., Maekawa, M., Kobayashi, K., Kojima, M., Nagato, Y., Sakakibara, H.,

NU

Kyozuka, J., 2007. Direct control of shoot meristem activity by a cytokinin-activating enzyme. Nature 445 (7128), 652-655. https://doi.org/10.1038/nature05504. Kutchan, T.M., 2001. Ecological arsenal and developmental dispatcher. The paradigm of secondary

MA

metabolism. Plant Physiol. 125 (1), 58-60. https://doi.org/10.1104/pp.125.1.58. Langmead, B., Trapnell, C., Pop, M., Salzberg, S.L., 2009. Ultrafast and memory-efficient alignment of

short

DNA

sequences

to

the

human

genome.

Genome

Biol.

10

(3),

1-10.

https://doi.org/10.1186/gb-2009-10-3-r25.

D

Li, H., Liang, W., Hu, Y., Zhu, L., Yin, C., Xu, J., Dreni, L., Kater, M.M., Zhang, D., 2011. Rice MADS6 Interacts with the Floral Homeotic Genes SUPERWOMAN1, MADS3, MADS58, MADS13,

PT E

and DROOPING LEAF in Specifying Floral Organ Identities and Meristem Fate. Plant Cell 23 (7), 2536-2552. https://doi.org/10.1105/tpc.111.087262. Li, S., Zhao, B., Yuan, D., Duan, M., Qian, Q., Tang, L., Wang, B., Liu, X., Zhang, J., Wang, J., Sun, J., Liu, Z., Feng, Y., Yuan, L., Li, C., 2013. Rice zinc finger protein DST enhances grain production

CE

through controlling Gn1a/OsCKX2 expression. Proc. Natl. Acad. Sci. U. S. A. 110 (8), 3167-3172. https://doi.org/10.1073/pnas.1300359110. Li, Y., Fan, C., Xing, Y., Jiang, Y., Luo, L., Sun, L., Shao, D., Xu, C., Li, X., Xiao, J., He, Y., Zhang,

AC

Q., 2011. Natural variation in GS5 plays an important role in regulating grain size and yield in rice. Nat. Genet. 43 (12), 1266-1269. https://doi.org/10.1038/ng.977. Liu, J.P., Van Eck, J., Cong, B., Tanksley, S.D., 2002. A new class of regulatory genes underlying the cause of pear-shaped tomato fruit. Proc. Natl. Acad. Sci. U. S. A. 99 (20), 13302-13306. https://doi.org/10.1073/pnas.162485999. Luo, A.D., Qian, Q., Yin, H.F., Liu, X.Q., Yin, C.X., Lan, Y., Tang, J.Y., Tang, Z.S., Cao, S.Y., Wang, X.J., Xia, K., Fu, X.D., Luo, D., Chu, C.C., 2006. EUI1, encoding a putative cytochrome P450 monooxygenase, regulates internode elongation by modulating gibberellin responses in rice. Plant Cell Physiol. 47 (2), 181-191. https://doi.org/10.1093/pcp/pci233. Luo, J., Liu, H., Zhou, T., Gu, B., Huang, X., Shangguan, Y., Zhu, J., Li, Y., Zhao, Y., Wang, Y., Zhao, Q., Wang, A., Wang, Z., Sang, T., Wang, Z., Han, B., 2013. An-1 Encodes a Basic Helix-Loop-Helix Protein That Regulates Awn Development, Grain Size, and Grain Number in 45

ACCEPTED MANUSCRIPT Rice. Plant Cell 25 (9), 3360-3376. https://doi.org/10.1105/tpc.113.113589. Magome, H., Nomura, T., Hanada, A., Takeda-Kamiya, N., Ohnishi, T., Shinma, Y., Katsumata, T., Kawaide, H., Kamiya, Y., Yamaguchi, S., 2013. CYP714B1 and CYP714B2 encode gibberellin 13-oxidases that reduce gibberellin activity in rice. Proc. Natl. Acad. Sci. U. S. A. 110 (5), 1947-1952. https://doi.org/10.1073/pnas.1215788110. Mao, H., Sun, S., Yao, J., Wang, C., Yu, S., Xu, C., Li, X., Zhang, Q., 2010. Linking differential domain functions of the GS3 protein to natural variation of grain size in rice. Sci. U. S. A. 107 (45), 19579-19584. https://doi.org/10.1073/pnas.1014419107. McSteen, P., 2009. Hormonal Regulation of Branching in Grasses. Plant Physiol. 149 (1), 46-55.

PT

https://doi.org/10.1104/pp.108.129056.

Meng, Q., Li, X., Zhu, W., Yang, L., Liang, W., Dreni, L., Zhang, D., 2017. Regulatory network and

RI

genetic interactions established by OsMADS34 in rice inflorescence and spikelet morphogenesis. J. Integr. Plant Biol. 59(9), 693-707. https://doi.org/10.1111/jipb.12594.

SC

Mimura, M., Nagato, Y., Itoh, J., 2012. Rice PLASTOCHRON genes regulate leaf maturation downstream of the gibberellin signal transduction pathway. Planta 235 (5), 1081-1089. https://doi.org/10.1007/s00425-012-1639-5.

NU

Miyoshi, K., Ahn, B.O., Kawakatsu, T., Ito, Y., Itoh, J.I., Nagato, Y., Kurata, N., 2004. PLASTOCHRON1, a timekeeper of leaf initiation in rice, encodes cytochrome P450. Proc. Natl. Acad. Sci. U. S. A. 101 (4), 875-880. https://doi.org/10.1073/pnas.2636936100. Annu.

Rev.

Cell

MA

Mockaitis, K., Estelle, M., 2008. Auxin Receptors and Plant Development: A New Signaling Paradigm. Dev.

Biol.

24

(24),

55-80.

https://doi.org/10.1146/annurev.cellbio.23.090506.123214. Moumeni, A., Satoh, K., Venuprasad, R., Serraj, R., Kumar, A., Leung, H., Kikuchi S., 2015.

D

Transcriptional profiling of the leaves of near-isogenic rice lines with contrasting drought tolerance at the reproductive stage in response to water deficit. BMC Genomics, 16(1), 1110.

PT E

https://doi.org/10.1186/s12864-015-2335-1. Nakagawa, M., Shimamoto, K., Kyozuka, J., 2002. Overexpression of RCN1 and RCN2, rice TERMINAL FLOWER 1/CENTRORADIALIS homologs, confers delay of phase transition and altered

panicle

morphology

in

rice.

Plant

J.

29

(6),

743-750.

CE

https://doi.org/10.1046/j.1365-313X.2002.01255.x. Nayar, S., Sharma, R., Tyagi, A.K., Kapoor, S., 2013. Functional delineation of rice MADS29 reveals its role in embryo and endosperm development by affecting hormone homeostasis. J. Exp. Bot.

AC

64 (14), 4239-4253. https://doi.org/10.1093/jxb/ert231. Oikawa, T., Kyozuka, J., 2009. Two-Step Regulation of LAX PANICLE1 Protein Accumulation in Axillary

Meristem

Formation

in

Rice.

Plant

Cell

21

(4),

1095-1108.

https://doi.org/10.1105/tpc.108.065425. Pattanayak, G.K., Tripathy, B.C., 2016. Modulation of biosynthesis of photosynthetic pigments and light-harvesting complex in wild-type and gun5 mutant of Arabidopsis thaliana during impaired chloroplast

development.

Protoplasma

253

(3),

747-752.

https://doi.org/10.1007/s00709-016-0958-y. Peer, W.A. and Murphy, A.S., 2007. Flavonoids and auxin transport: modulators or regulators? Trends Plant Sci. 12 (12), 556-563. https://doi.org/10.1016/j.tplants.2007.10.003. Piasecka, A., Jedrzejczak-Rey, N., Bednarek, P., 2015. Secondary metabolites in plant innate immunity: conserved function of divergent chemicals.

New

Phytol.

206

(3),

948-964. 46

ACCEPTED MANUSCRIPT https://doi.org/10.1111/nph.13325. Qi, P., Lin, Y., Song, X., Shen, J., Huang, W., Shan, J., Zhu, M., Jiang, L., Gao, J., Lin, H., 2012. The novel quantitative trait locus GL3.1 controls rice grain size and yield by regulating Cyclin-T1;3. Cell Res. 22 (12), 1666-1680. https://doi.org/10.1038/cr.2012.151. Quattrocchio, F., Wing, J.F., Van Der Woude, K., Mol, J.N.M., Koes, R., 1998. Analysis of bHLH and MYB domain proteins: Species-specific regulatory differences are caused by divergent evolution of

target

anthocyanin

genes.

Plant

J.

13

(4),

475-488.

https://doi.org/10.1046/j.1365-313X.1998.00046.x. Shikanai, T., Fujii, S., 2013. Function of PPR proteins in plastid gene expression. RNA Biol. 10 (9),

PT

1446-1456. https://doi.org/10.4161/rna.25207.

Shomura, A., Izawa, T., Ebana, K., Ebitani, T., Kanegae, H., Konishi, S., Yano, M., 2008. Deletion in a

RI

gene associated with grain size increased yields during rice domestication. Nat. Genet. 40 (8), 1023-1028. https://doi.org/10.1038/ng.169.

SC

Si, L., Chen, J., Huang, X., Gong, H., Luo, J., Hou, Q., Zhou, T., Lu, T., Zhu, J., Shangguan, Y., Chen, E., Gong, C., Zhao, Q., Jing, Y., Zhao, Y., Li, Y., Cui, L., Fan, D., Lu, Y., Weng, Q., Wang, Y., Zhan, Q., Liu, K., Wei, X., An, K., An, G., Han, B., 2016. OsSPL13 controls grain size in

NU

cultivated rice. Nat. Genet. 48 (4), 447-456. https://doi.org/10.1038/ng.3518. Song, X., Huang, W., Shi, M., Zhu, M., Lin, H., 2007. A QTL for rice grain width and weight encodes a previously unknown RING-type E3 ubiquitin ligase. Nat. Genet. 39 (5), 623-630.

MA

https://doi.org/10.1038/ng2014.

Sundaresan, V., 2005. Control of seed size in plants. Proc. Natl. Acad. Sci. U. S. A. 102 (50), 17887-17888. https://doi.org/10.1073/pnas.0509021102. Tabuchi, H., Zhang, Y., Hattori, S., Omae, M., Shimizu-Sato, S., Oikawa, T., Qian, Q., Nishimura, M.,

D

Kitano, H., Xie, H., Fang, X., Yoshida, H., Kyozuka, J., Chen, F., Sato, Y., 2011. LAX PANICLE2 of Rice Encodes a Novel Nuclear Protein and Regulates the Formation of Axillary Meristems.

PT E

Plant Cell 23 (9), 3276-3287. https://doi.org/10.1105/tpc.111.088765. Tan, L., Li, X., Liu, F., Sun, X., Li, C., Zhu, Z., Fu, Y., Cai, H., Wang, X., Xie, D., Sun, C., 2008. Control of a key transition from prostrate to erect growth in rice domestication. Nat. Genet. 40 (11), 1360-1364. https://doi.org/10.1038/ng.197.

CE

Tariq, R., Wang, C., Qin, T., Xu, F., Tang, Y., & Gao, Y., Ji Z., Zhao K., 2018. Comparative transcriptome profiling of rice near-isogenic line carrying xa23 under infection of xanthomonas oryzae pv. oryzae. Int. J. Mol. Sci. 19(3), 717. https://doi.org/10.3390/ijms19030717.

AC

Terao, T., Nagata, K., Morino, K., Hirose, T., 2010. A gene controlling the number of primary rachis branches also controls the vascular bundle formation and hence is responsible to increase the harvest

index

and

grain

yield

in

rice.

Theor.

Appl.

Genet.

120

(5),

875-893.

https://doi.org/10.1007/s00122-009-1218-8. Wang, L., Li, P., Brutnell, T.P., 2010. Exploring plant transcriptomes using ultra high-throughput sequencing. Briefings Funct. Genomics 9 (2), 118-128. https://doi.org/10.1093/bfgp/elp057. Wang, L., Sun, S., Jin, J., Fu, D., Yang, X., Weng, X., Xu, C., Li, X., Xiao, J., Zhang, Q., 2015. Coordinated regulation of vegetative and reproductive branching in rice. Proc. Natl. Acad. Sci. U. S. A. 112 (50), 15504-15509. https://doi.org/10.1073/pnas.1521949112. Wang, L., Zeng, X., Zhuang, H., Shen, Y., Chen, H., Wang, Z., Long, J., Ling, Y., He, G., Li, Y., 2017. Ectopic expression of OsMADS1 caused dwarfism and spikelet alteration in rice. Plant Growth Regul. 81 (3), 433-442. https://doi.org/10.1007/s10725-016-0220-9. 47

ACCEPTED MANUSCRIPT Wang, S., Li S., Liu Q., Wu K., Zhang J., Wang S., Wang Y., Chen X., Zhang Y., Gao C., Wang F., Huang H., Fu X., 2015. The OsSPL16-GW7 regulatory module determines grain shape and simultaneously improves rice yield and grain quality. Nat. Genet. 47(8): 949-954. https://doi.org/10.1038/ng.3352. Wang, S., Wu, K., Yuan, Q., Liu, X., Liu, Z., Lin, X., Zeng, R., Zhu, H., Dong, G., Qian, Q., Zhang, G., Fu, X., 2012. Control of grain size, shape and quality by OsSPL16 in rice. Nat. Genet. 44 (8), 950-954. https://doi.org/10.1038/ng.2327. Weng, J., Gu, S., Wan, X., Gao, H., Guo, T., Su, N., Lei, C., Zhang, X., Cheng, Z., Guo, X., Wang, J., Jiang, L., Zhai, H., Wan, J., 2008. Isolation and initial characterization of GW5, a major QTL with

rice

grain

width

and

weight.

Cell

Res.

https://doi.org/10.1038/cr.2008.307.

18

(12),

PT

associated

1199-1209.

RI

Winkel-Shirley, B., 2001. Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology. Plant Physiol. 126 (2), 485-493. https://doi.org/10.1104/pp.126.2.485.

SC

Wu, Y., Wang, Y., Mi, X., Shan, J., Li, X., Xu, J., Lin, H., 2016. The QTL GNP1 Encodes GA20ox1, Which Increases Grain Number and Yield by Increasing Cytokinin Activity in Rice Panicle Meristems. PLoS Genet. 12 (10): e1006386. https://doi.org/10.1371/journal.pgen.1006386.

NU

Xing, Y., Zhang, Q., 2010. Genetic and molecular bases of rice yield. Plant Biology 61 (61), 421-442. https://doi.org/10.1146/annurev-arplant-042809-112209.

Yang, D., Yao, J., Mei, C., Tong, X., Zeng, L., Li, Q., Xiao, L., Sun, T., Li, J., Deng, X., Lee, C.M.,

MA

Thomashow, M.F., Yang, Y., He, Z., He, S.Y., 2012. Plant hormone jasmonate prioritizes defense over growth by interfering with gibberellin signaling cascade. Proc. Natl. Acad. Sci. U. S. A. 109 (19), 1192-1200. https://doi.org/10.1073/pnas.1201616109. Yoshida, A., Ohmori, Y., Kitano, H., Taguchi-Shiobara, F., Hirano, H., 2012. ABERRANT SPIKELET regulation

of

meristem

D

AND PANICLE1, encoding a TOPLESS-related transcriptional co-repressor, is involved in the fate

in

rice.

Plant

J.

70

(2),

327-339.

PT E

https://doi.org/10.1111/j.1365-313X.2011.04872.x. Yoshida, A., Sasao, M., Yasuno, N., Takagi, K., Daimon, Y., Chen, R., Yamazaki, R., Tokunaga, H., Kitaguchi, Y., Sato, Y., Nagamura, Y., Ushijima, T., Kumamaru, T., Iida, S., Maekawa, M., Kyozuka, J., 2013. TAWAWA1, a regulator of rice inflorescence architecture, functions through the

CE

suppression of meristem phase transition. Proc. Natl. Acad. Sci. U. S. A. 110 (2), 767-772. https://doi.org/10.1073/pnas.1216151110. Zhang, X., Wang, J., Huang, J., Lan, H., Wang, C., Yin, C., Wu, Y., Tang, H., Qian, Q., Li, J., Zhang,

AC

H., 2012. Rare allele of OsPPKL1 associated with grain length causes extra-large grain and a significant yield increase in rice. Proc. Natl. Acad. Sci. U. S. A. 109 (52), 21534-21539. https://doi.org/10.1073/pnas.1219776110. Zhang, Y., Yu, H., Liu, J., Wang, W., Sun, J., Gao, Q., Zhang, Y., Ma, D., Wang, J., Xu, Z., Chen, W., 2016. Loss of function of OsMADS34 leads to large sterile lemma and low grain yield in rice (Oryza sativa L.). Mol. Breed. 36, 147. https://doi.org/10.1007/s11032-016-0578-4. Zhao, Y., 2008. The role of local biosynthesis of auxin and cytokinin in plant development. Curr. Opin. Plant Biol. 11 (1), 16-22. https://doi.org/10.1016/j.pbi.2007.10.008. Zhao, Y., Cheng, S., Song, Y., Huang, Y., Zhou, S., Liu, X., Zhou, D., 2015. The Interaction between Rice ERF3 and WOX11 Promotes Crown Root Development by Regulating Gene Expression Involved

in

Cytokinin

Signaling.

Plant

Cell

27

(9),

2469-2483.

https://doi.org/10.1105/tpc.15.00227. 48