Huirn Liu, Gng Li*, Xiujun Yng Hendrik N.J. Kuijer, Wnqi Ling,Ding Zhng*

aJoint International Research Laboratory of Metabolic & Developmental Sciences, Shanghai Jiao Tong University-University of Adelaide Joint Centre for Agriculture and Health, State Key Laboratory of Hybrid Rice, School of Life Sciences and Biotechnology, Shanghai Jiao Tong University, Shanghai 200240,China

bSchool of Agriculture,Food and Wine, The University of Adelaide,South Australia 5064,Australia

Keywords:Inflorescence meristem Transcriptome Gene expression Hormones Barley

ABSTRACT The shape of an inflorescence varies among cereals,ranging from a highly branched panicle in rice to a much more compact spike in barley (Hordeum vulgare L.) and wheat (Triticum aestivum L.). However, little is known about the molecular basis of cereal inflorescence architecture. We profiled transcriptomes at three developmental stages of the barley main shoot apex-spikelet initiation,floral organ differentiation,and floral organ growth-and compared them with those from vegetative seedling tissue. Transcript analyses identified 3688 genes differentially transcribed between the three meristem stages, with a further 1394 genes preferentially expressed in reproductive compared with vegetative tissue. Coexpression assembly and Gene Ontology analysis classified these 4888 genes into 28 clusters, revealing distinct patterns for genes such as transcription factors, histone modification, and cell-cycle progression specific for each stage of inflorescence development. We also compared expression patterns of VRS (SIX-ROWED SPIKE) genes and auxin-,gibberellic acid- and cytokinin-associated genes between two-rowed and six-rowed barley to describe regulators of lateral spikelet fertility. Our findings reveal barley inflorescence phase-specific gene expression, identify new candidate genes that regulate barley meristem activities and flower development, and provide a new genetic resource for further dissection of the molecular mechanisms of spike development.

1. Introduction

Reproductive development is essential for plant propagation and critical for crop yield. Agriculturally important cereals such as rice(Oryza sativa L.),barley(Hordeum vulgare L.),wheat(Triticum aestivum L.), maize (Zea mays L.), and sorghum(Sorghum bicolor L.) belong to the Poaceae family of grasses,one of the largest angiosperm families with some 11,000 extant species. These cereals, and grasses more broadly,display highly diverse inflorescence architecture in the arrangement of their spikelets and flowers, so that the molecular mechanisms that regulate inflorescence architecture and thereby grain yield are an important research focus[1,2].

In the grass main shoot apex (MSA), the organization of spikelet subunits within the inflorescence is determined by a complex hierarchy of specialized meristems [1]. The basal unit of a grass inflorescence is a spikelet,which is a short and condensed branch containing leaflike structures that enclose one or more florets (flowers). Inflorescences of the tribe Triticeae, including barley and wheat, exhibit a raceme-like branchless architecture (spike) with a shortened stem axis called the rachis. Barley develops three single spikelets at each rachis node and each spikelet bears only one floret. In two-rowed barley varieties, only the central spikelet in each triplet is fertile and can produce a seed, while the lateral spikelets are much reduced in size. In six-rowed barley varieties, both central and lateral spikelets can form seeds.Previous genetic studies have identified five independent SIXROWED SPIKE(VRS)genes associated with complete or partial changes in the fertility of the lateral spikelets of barley. VRS1 and VRS4 encode a homeodomain-leucine zipper (HD-Zip)class-I homeobox transcription factor and a lateral organ boundary (LOB) transcription factor, respectively. Loss-offunction mutations in both VRS1 and VRS4 lead to a full sixrowed phenotype, indicating that VRS1 and VRS4 act as key inhibitors of lateral spikelet fertility and development [3,4].

Increasing evidence has revealed that numerous other elements play key roles in determining inflorescence fate and patterning in plants. The CLAVATA-WUSCHEL (CLV-WUS)circuit is a negative-feedback loop that consists of the stem cell-promoting transcription factor WUS, the differentiationdriven peptide CLV3, and the ligand-receptor kinase of CLV1;mutations in CLV genes frequently lead to enlarged shoot apical meristems and additional organs in flowers [1,5]. A homeobox transcription factor, KNOTTED 1 (KN1), positively regulates shoot apical meristem establishment and inflorescence meristem maintenance in rice, maize, and Arabidopsis[1,6]. The plant hormones auxin, gibberellic acid (GA) and cytokinin (CK) play essential roles in MSA activity and inflorescence patterning by regulating cell proliferation and differentiation; disruptions of hormone homeostasis modulate inflorescence development and affect yield potential in crops [1,7]. For example, accumulation of cytokinin in IMs by perturbations in its degradative pathway leads to increased branch and spikelet numbers and grain yield[8].

MADS-box family transcription factors, which comprise the vast bulk of the ABCDE model of floral organ formation,work in combination with each other to regulate plant inflorescence, flower morphology and meristem activities.Key MADS-box genes include the SEPALLATA (SEP) clade and MADS-domain homeotic proteins APETALA1(AP1), APETALA3(AP3), PISTILLATA (PI), and AGAMOUS (AG) [1,9,10]. Overexpression of barley APETALA1 (AP1)-like VRN1 (VERNALIZATION1) (HvMADS14) and ectopic expression of SHORT VEGETATIVE PHASE (SVP)-like BM1 (HvMADS47) and BM10(HvMADS22) cause abnormal MSA activity and defects in spike pattern formation [11,12]. Expression of APELATA2(AP2)-like transcription factor is regulated by microRNA172 to determine the density of grains on the barley spike[13].

These previous results indicate essential roles for a large number of genes that regulate barley inflorescence development and morphogenesis.Agriculturally, barley grain yield is determined by three main components:number of spikes per square meter,kernel number per spike,and individual kernel weight. Kernel number per spike is determined early in reproductive development, so that the establishment phase is crucial for yield potential [14]. During establishment, MSA initiation(double ridge(DR)stage,Waddington stage 2.0) and floral differentiation(awn primordia(AP)stage,W3.5)are two key points that determine final grain setting [15]. The late reproductive phase, especially the green anther (GrA) stage(W8-W8.5),is also critical for grain yield,as this phase is vital for spikelet survival[15].

Recently, several studies have investigated gene expression in selected developmental stages of inflorescence development. The International Barley Genome Sequencing Consortium(IBGSC,2012)[16]reported the expression profiles of 26,159 high-confidence (HC) genes in inflorescence late developmental stages(W6-W7 and white anther stages)of the six-rowed barley cultivar Morex. Mascher et al. [17] updated the high-quality reference genome of Morex with 39,734 highconfidence annotated genes and showed that a large number of genes are undetected in the previous study.Digel et al.[18]described photoperiod- and PHOTOPERIOD1 (Ppd-H1)-dependent gene expression in inflorescence development and flower fertility of two-rowed barley. However, patterns of gene expression in earlier developmental stages critical for inflorescence patterning,and their co-expression in networks,remain largely unknown.

We performed high-throughput RNA-seq on six-rowed barley MSAs from spikelet initiation (DR stage) to flower organ differentiation(AP stage)and growth(GrA stage),and in the leaves of two-week-old seedlings (vegetative control), to investigate transcriptional progression and possible regulatory networks during inflorescence development. We reconstructed the expression network of stage-specific transcripts in barley MSAs and identified the involvement of pivotal genetic regulators, plant hormones, and cell-cycle genes in barley inflorescence development.

2. Materials and methods

2.1. Plant growth and tissue collection

Two spring barley cultivars,Morex and Scarlett,were grown in the greenhouse under 22 °C/17 °C day/night, 16 h/8 h light/dark,and 50%relative humidity(The Plant Accelerator at The University of Adelaide, Adelaide, Australia). For RNA-seq experiments in Morex and qRT-PCR in both cultivars, the developmental stages were determined by dissection under a stereomicroscope (S8 APO, Leica Microsystems, Germany).Four stages of tissues were collected: the double ridge meristem stage (DR), the awn primordium meristem stage(AP),the green anther stage(GrA),and two-week-old seedling leaves.The MSAs at the DR stage were collected when a plant had 3-4 leaves and the first stem node was visible(19-20 days after sowing). The AP stage was marked by the emergence of the second node and floret meristems (up to stamen primordia) [16,18]. In the GrA stage, the late reproductive stage, the spikes were 4-5 cm in length. The leaves of seedlings were collected from plants as vegetative materials two weeks after sowing seeds in soil.

All sample stages were immediately frozen in liquid nitrogen and stored at -80 °C until RNA extraction. For each MSA stage,15-40 spikes were pooled for each of two biological replicates.

2.2. RNA isolation and sequencing library construction

Total RNA was extracted from tissues using TRIzol (Life Technologies) and purified using a RNeasy Micro Kit (Qiagen,Germany) for RNA sequencing and qRT-PCR. Total mRNA of four tissues from Morex was enriched with oligo (dT) primer and fragmented by mixing with fragmentation buffer. The residual DNA was removed with a RNase-free kit (Ambion)and the RNA concentration and integrity were determined with a 2100 Bioanalyzer (Agilent) before RNA library preparation for RNA sequencing. Only transcripts longer than 200 nucleotides were further processed for RNA sequencing.

Stranded RNA-seq libraries were prepared with the NEXTflex Rapid Directional RNA-Seq Kit with NEXTflex DNA Barcodes (Bioo Scientific, manual version 14.09) using 1 μg total RNA of each sample to generate a library using NEXTfle beads. Libraries were normalized, pooled at equimolar concentrations, and diluted to 10 pmol L-1. Pools were clustered with the HiSeq Rapid PE Cluster Kit v2 (Illumina). During quality control (QC) steps, the quantity and quality of the sample libraries were assessed with an Agilent 2100 Bioanalyzer and an ABI StepOnePlus Real-Time PCR System.

2.3. RNA sequencing and data analysis

The cDNA libraries were added to flow cells of an Illumina HiSeq 4000 for paired-end sequencing at the Beijing Genomics Institute. After sequencing, the quality of raw sequencing reads for all samples was evaluated with FastQC (version 0.11.4; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads with adaptors and having low quality (with the percentage of bases with a quality score of lower than 15 exceeding 20% in a read) were removed using Trimmomatic [19] software, and reads composed of ＞5% unknown bases (labeled as N) were discarded. The clean reads were mapped to the barley reference genome (the Morex genome released in 2017: http://webblast.ipk-gatersleben.de/barley_ibsc/) [18] using HISAT2 2.0.0 [20]. The mapped reads were sorted into two groups: singly and multiply mapped reads. HTSeq [21] was used to count reads mapped to each gene. Read counts were normalized to fragments per kilobase per million (FPKM) to obtain relative levels of transcript expression. Read counts and FPKMs of all barley transcripts are shown in Table S1. Correlation analysis of biological replicates was performed, and Pearson correlation coefficients were calculated among the eight samples to validate the quality of biological replicates. Each scatter plot was standardized to mean fold coverage of the FPKM, plotted as log (mean + 1), and a linear model was calculated and plotted. For a principal component analysis (PCA), regularized-logarithm transformation (rlog)[22]for count data of the samples was used.

2.4. Differential gene expression analysis

Analysis of differential gene expression between developmental stages was performed with the R package DESeq2 [22]using a negative binomial distribution model. Raw data counts of FPKM (per million mapped reads) were normalized and then transformed for estimation of their mean and variance. Results for pairwise comparisons of any two development stages were extracted to investigate the effects of genes. Benjamini-Hochberg (BH) [23] adjustment was performed to compute adjusted P value. Differentially expressed genes (DEGs) were considered significant if the fold change log2Ratio ≥2, adjusted P-value ＜0.05 and false discovery rate (FDR) ≤0.001 (Tables S2, S3). Hierarchical clustering for DEGs was performed with the heatmap package in R. Annotation of all DEGs was based on BLASTx alignment results against protein databases of Arabidopsis(TAIR10_peptide; http://www.arabidopsis.org/), rice(MSU7_peptide; http://rice.plantbiology.msu.edu/), and Brachypodium (Brachy3.1_peptide; https://phytozome.jgi.doe.gov/) with E ＜10-6(Fig. S1-A, B). Distribution and density of DEGs were plotted along the barley chromosomes using the chromPlot R package [24]. The genomic data of each chromosome were calculated to construct histograms representing numbers of genes in 1-Mb bins.Expression data of DR and AP in Scarlett were derived from Digel et al. [18] (Table S4).Sequences of DR and AP in Scarlett were extracted from the Barley project ftp address (ftp://ftpmips.helmholtzmuenchen.de/plants/barley/public_data/genes/, sequence version of March 23rd, 2012) and barley UniGenes (ftp://ftp.ncbi.nih.gov/repository/UniGene/Hordeum_vulgare/) using a python code. BLAST alignment of DR and AP in the Scarlett sequence against the Morex genome sequence was used to establish correspondences between genes in Scarlett and Morex.

2.5. Co-expression clusters

Co-expression analysis of 4888 DEGs was performed with the R package MBCluster.Seq version 1.0, based on a negative binomial distribution model [25]. DEGs were clustered based on the similarity of their expression patterns in MSAs and the optimal number of co-expression clusters was determined in the range of 5 to 80 with the Expectation-Maximization (EM)algorithm and negative binomial models. The EM algorithm was used to estimate the negative binomial model parameters and cluster membership. Negative binomial models were applied to implement hierarchical clustering and a hybrid tree. Based on the visualization of the hybrid tree, the final number of clusters was chosen to achieve a high average probability.

2.6. Gene Ontology classification

The DEGs annotations of Arabidopsis orthologs were subjected to BLASTx using BLAST+ (version 2.6.0) against the TAIR 10 protein database. Alignment results were used to assign GO annotations to barley DEGs. Barley gene-to-GO associations were captured by a python script from Arabidopsis Gene Ontology (GO) annotations (http://www.arabidopsis.org/download/index.jsp). The R package clusterProfiler [26] was used for Gene Ontology enrichment analysis of transcripts in MSAs_1394 (highly expressed genes in MSAs), SE_3731 (highly expressed genes in seedlings) and three tissue-specific CLUSTERs, I-III. Transcripts with FPKM＞50 in the DR, AP, GrA, and SE stages were selected for Gene Ontology enrichment analysis. GO terms with a corrected false discovery rate of ＜0.05 were considered to be significantly enriched. AgriGO2 (http://amigo.geneontology.org/visualize?mode=client_amigo) was used for GO classification analysis and the identification of pathways of phase-specific genes. The percentage of differential transcripts for each GO classification was calculated with reference to all the overexpressed transcripts in a given GO analysis. Clustering was performed in R using the k-means function, where k = 24 in the cluster package by Euclidean distance. REVIGO [27] was used to summarize and visualize the GO term results as tree maps of highly expressed genes. Using SimRel semantic similarity measures, terms were clustered at a specified similarity cutoff of 0.9 (a “l(fā)arge” REVIGO set). The detailed biological processes of representative terms are listed in Table S5.

2.7. Verification of gene expression by qRT-PCR analysis

cDNA was synthesized from 4 μg of total RNA using a reverse transcriptase kit (Bio-Rad). SYBR Green Mix (Kapa) was used for qRT-PCR on a QuantStudio Flex 6 (Life Technologies) machine as previously described [28]. The qRT-PCR data for each target gene are presented as average expression levels over three biological replicates, with two technical replicates per reaction, relative to the expression levels of the HvActin7 reference gene. Primer sequences are listed in Table S6.

3. Results

3.1. The yield of RNA-seq global data

To identify candidate genes regulating spike development in barley, we conducted whole-transcriptome expression profiles in MSAs throughout development in Morex barley, from the early DR stage(W2.0)(Waddington stage)[16]and AP stage(W3.5) to the late GrA (W8-W8.5) stage, and compared them with the transcriptome in vegetative tissues (2-week seedlings)(Fig.1-A).

We sequenced total messenger RNA from DR,AP,GrA,and seedling tissues, generating over 256 million reads of all libraries (Table 1; Fig. S1). Of the total reads from the four samples, each sample yielded at least 31.4 million reads, of which ＞95.2% could be mapped onto the barley reference genome,with over 70%mapping to a unique genomic location(Table 1). There was a high correlation between biological replicates for each developmental stage (Fig. S2-A), and distinct differences between the four developmental stages,as indicated by PCA(Fig.S2-B).

3.2. Global gene expression profiles of reproductive and vegetative tissues

Read counts from the raw data were normalized to fragments per kilobase per million(FPKM) to give the relative levels of transcript expression of the 39,734 high-confidence genes identified in the barley genome (Table S1). A total of 10,717 (26.97%) genes, whose locations were concentrated toward the distal ends of chromosome arms (Fig. S3), were differentially expressed among the four stages(Tables S2,S3).Generally, gene expression was most divergent between the reproductive and the vegetative stages (Fig. S4-A), with more DEGs between the three developing MSA stages and the seedling stage than between the reproductive stages (Fig. S4-B,C;Tables S2, S3).

Within the reproductive stages, GrA was most different from the other two, with genes generally showing upregulation in this stage (Fig. S4-B). A distinct subset of 3688 genes with differential expression between the three reproductive stages was selected,including 2764 transcripts whose expression increased as development progressed (from DP to AP to GrA), 770 whose expression decreased, and 154 whose expression either peaked or troughed in the middle AP stage(Fig. 1-A). Between the three reproductive phases, relatively few(703)genes were differentially expressed in the AP vs.DR set (Fig. 1-A), whereas DEGs between GrA vs. DR and GrA vs.AP showed a large proportion of overlap (1867 DEGs, 50.6%;Fig. 1-B), providing more evidence that the DR and AP stages are more similar to each other than to the GrA stage(see also Fig. S2-B). Only 106 genes (2.8%) showed significantly altered expression among comparisons between all three MSA stages(Fig. 1-B; Table S2), suggesting elaborate gene switches and tissue-specific gene expression patterns during barley inflorescence development.

There was a much larger difference between the reproductive tissues and the vegetative seedling stage.A total of 10,431 genes showed differential expression in the DR, AP, and GrA stages compared with seedlings, with 6564 showing reduced expression in MSAs(Fig.1-A).Almost half(5171,49.6%)of the DEGs were shared among all three comparison sets,of which 1394 genes were upregulated in all MSA stages (and named MSA_1394), and another pool of 3731 genes (named SE_3731)that were expressed more highly in seedlings (Fig. 1-C; Table S7). The inflorescence-enriched gene group MSA_1394, comprising transcription factors, kinases,and transporters (Table S8), showed a relatively stable expression pattern during the developmental process of the MSA.Only a small proportion of MSA_1394 genes showed differential expression between the three reproductive stages (Fig. S5-A-C), with only two genes also present in the 106 DEGs identified in all reproductive tissues(Fig.S5-D).

3.3. Expression of key regulators in barley developing inflorescence

For each tissue, we analyzed barley homologs of rice, maize and Arabidopsis genes with a known function in inflorescence development. A heat map of expression values showed that most of those orthologs were expressed at a higher level in reproductive than in vegetative tissues (Fig. 2-A; Table S9). These genes were categorized into three subgroups based on their roles in the plant meristem (Fig. 2-A).

Fig.1- Overview of differential gene expression analysis in barley developing inflorescence and two-week seedlings. DR,double-ridge stage(W2.0,Waddington stage);AP,awn primordium stage(W3.5);GrA, green anther stage.Transcripts with higher and lower abundance are indicated by red and green arrows,respectively.Transcripts with inconsistent expression are indicated by gray arrows.(A)DEGs between inflorescence tissues at three developmental stages and in seedlings of Morex.Bars in DR and AP,500 μm;bar in GrA, 1 cm.(B)Venn diagram showing the overlap of transcripts that were differentially expressed between the three comparison sets of the three reproductive stages.The box on the right indicates transcripts cochanged in all three comparisons.(C)Venn diagram representing the overlap of transcripts that were differentially expressed between the three comparison sets of each reproductive stage and the seedling stage.The box on the right indicates transcripts co-changed in all three comparisons.

Genes involved in meristem signaling, namely barley CLV1-like, CLV2-like, and WUS-like homologs, showed varied expression patterns. From DR to GrA, HvCLV2-like showed a consistent expression pattern, but the expression of HvCLV1 homologs and HvWUS homologs changed oppositely over time (Fig. 2-A; Table S9). To verify the RNA-seq data, we performed reverse transcription quantitative PCR (qRT-PCR)analysis for selected representative genes, confirming the negative correlation between HvCLV1-like and HvWUS-like expression, which suggests the possible regulation of these two genes by means of a negative feedback loop(Fig.2-B).

Several genes encode transcription factors, such as LFY1(LEAFY1, plant-specific), KN1 (homeobox), and ANT(AINTEGUMENTA, AP2-domain), that have been shown[1,29,30] to control inflorescence meristem identity and sizein plants.High expression levels of the three barley homologs of LFY, KN1, and ANT were detected in the DR stage in both RNA-seq and qRT-PCR results,and their transcript abundance generally decreased during barley spike development (Fig. 2-A, B), suggesting a conserved role in MSA transition from vegetative to reproductive meristem.

Table 1-Sequencing yields and alignments.Each library was sequenced in two lanes on an Illumina sequencer.

Regulatory factors comprised the largest subgroup of genes with known functions in inflorescence development. The expression of HvAP1 (HvMADS14) showed a slight increase from the DR to AP stages of the meristem, whereas the expression of HvLFY,a potential upstream regulator of HvAP1,decreased (Fig. 2-A, B; Table S9). This expression pattern suggests a negative effect of HvLFY on HvAP1 expression in barley, in direct contrast to Arabidopsis, where LFY acts as a positive regulator of AP1 expression [1,31]. Expression of another inflorescence regulator, HvTEL1 (TERMINAL EAR1)[32], encoding a putative RNA-binding protein, peaked at AP stage,suggesting a specific function of HvTEL1 in the spikelet meristem (Fig. 2-A; Table S9). The accumulation of HvFT-like(FLOWERING LOCUS T-like)from DR to GrA was an indicator of phase transition from vegetative stage to reproductive stage,given that FT homologs have been identified as direct stimulators of flowering in plants (Fig. 2-A, B) [33]. In maize,FEA3 (FASCIATED EAR3), encoding a leucine-rich-repeat receptor,and FEA4,encoding a bZIP transcription factor,control inflorescence meristem size [34,35]. Decreased expression of HvFEA4 during barley spike development suggested its function in early inflorescence formation, whereas increased expression of HvFEA3-like suggests that the regulatory relationship between FEA3 and WUS in maize is conserved in barley (Fig. 2-A, B) [33]. Overall, homologs of key genes controlling the inflorescence meristematic activity in model crops are likely to regulate similar developmental processes in barley.

3.4. MADS family genes may have conserved and divergent roles in early barley inflorescence development

Fig.2-Expression profiling of inflorescence development related regulators in barley.(A)Heat map of DEGs regulating barley spike development and meristem identity.Red represents overexpression and blue represents underexpression of genes in the first-named relative to the second-named tissue in each set.(B)qRT-PCR expression data of selected developmental genes in the DR,AP and GrA stages.All results were normalized to HvActin7 and then converted to a percentage of total expression across all three samples. All results are shown as mean ± SD of 3 replicates. DR,double-ridge stage;AP,awn primordium stage;GrA,green anther stage;SE,seedlings.CLV,CLAVATA;WUS,WUSCHEL;LFY,Leafy;KN1,Homeobox protein knotted-1;ANT,Aintegumenta;TEL1,Terminal EAR1-like 1;BOP2,BLADE-ON-PETIOLE2;SPL,SQUAMOSA PROMOTER BINDING PROTEIN-LIKE;CYP78A9,Cytochrome 78A9;bHLH,Basic helix-loop-helix;YSL,Yellow Stripe-Like;FT,FLOWERING LOCUS T;FEA,FASCIATED EAR.

Transcription factors, such as the MADS, KNOTTED-like homeobox (KNOX) and AP2 families, are critical for early inflorescence development in cereals [1]. Among the 10,431 DEGs between reproductive and vegetative tissues (Fig. 1-A),several families of transcription factors were identified (Fig.S6). A large number of transcripts of bHLH (basic helix-loophelix), MYB (myeloblastosis) and NAC (NAM, ATAF, and CUC)genes were recognized in barley MSAs, and the expression of most bHLHs and NACs decreased from DR to GrA, while the transcriptional levels of most MYBs increased in the GrA stage(Fig. S6-A-C). Inflorescence/floret meristem regulators KNs,BTB/POZs (Broad-Complex, Tramtrack and Bric-à-brac/Poxvirus and Zinc finger) [36] and AP2s were preferentially expressed at the DR stage (Fig. S6-D-F), reflecting the roles of these transcription factors in barley MSA initiation.

Homologs of MADS-box genes play pleiotropic roles, from reproductive phase initiation to floret organogenesis in the ABCDE genetic model in plants[1].We identified 31 MADS-box genes with altered expression by comparisons among DR,AP,GrA, and seedlings (Fig. 3-A). Both RNA-seq and qRT-PCR results showed that only SVP-like genes (HvMADS55,HvMADS47) and one AGL12-like gene (HvMADS26) were preferentially expressed in seedlings (Figs. 3-A, B, S7-A, F).ABCDE homologs were expressed strongly in barley MSAs(Figs. 3-A, B, S7-B-E), indicating the essential roles of these genes in barley MSA and spikelet/floret development.However, expression levels between different classes of MADS-box genes varied. B-class genes (HvMADS2, 4, 16) were expressed much more highly than B-sister genes(HvMADS29,30,31),and the expression of C-class genes(HvMADS3,58)was higher than that of D-class genes (HvMADS13, 21) in the AP and GrA stages(Fig.S7-C,D).

Fig.3- Expression characteristics of the MADS gene family in the barley developing inflorescence. (A)Heat map of MADS transcription factor expression profiling in barley MSA stages and seedlings.Red represents overexpression and blue represents underexpression of genes in the first-named relative to the second named tissue in each set.(B)qRT-PCR expression of selected developmental genes in DR,AP,and GrA stages. All results were normalized to HvActin7 and then converted to a percentage of total expression across all three samples.All results are shown as mean ± SD of 3 replicates.DR,double ridge;AP,awn primordium;GrA,green anther;SE,seedlings.

Many MADS family genes in rice and wheat, like MADS1,MADS5, MADS14(AP1),are highly and dynamically expressed in early stages of reproductive organs [10,11,37]. Among the identified MADS genes in barley, only SVP-like genes(HvMADS55, 47), two A-class genes (MADS15, 18), and one Eclass gene (HvMADS34) were expressed at DR stage (Fig. S7-B,E, F), while HvMADS14 (AP1) was consistently expressed (Fig.2-B), suggesting the divergent roles of MADS genes in early inflorescence initiation among grasses.

3.5. MSA-enriched genes control cell behaviors in the barley inflorescence

From RNA-seq data of 39,734 HC genes among DR, AP, and GrA and in seedlings (Table S1), we identified 23 genes highly expressed (FPKM ＞150) in at least one reproductive stage but undetectable in seedlings (Tables 2; S1). qRT-PCR validated our finding that the expression values of these genes were 3-1300-fold higher in MSAs than in seedlings (Table S10). Some of these genes function in inflorescence transition, spikelet specification, and meristem activity, e.g., squamosa promoter-binding-like (SPL) genes SPL17-L and SPL4-L, homeobox knotted-1-like gene KN1 (KNOX), and growth-regulating factor1-interacting factor (GIF1), and showed high expression in DR and AP stages (Table S10). HvMADS32, the homolog of rice OsMADS32 that plays roles in lodicule identity and stamen development [38], the only MADS-box gene unique to grasses, showed the highest expression in the barley MADS gene family, indicating the critical role of this gene in the early reproductive meristem.

3.6. Gene clustering reveals phase-specific genes responsible for inflorescence development

To assign the functions of highly expressed genes in MSAs, GO term enrichment was performed on the highly expressed genes in each tissue (FPKM ＞ 50) (Table S5). The tree maps of biological processes showed that DR, AP and GrA shared some GO terms, such as chromatin silencing, DNA/RNA, and protein modification. However, the AP stage showed a significant enrichment in GO terms of cell proliferation and meristem organization, such as meristem transition, floral meristem identity, and determinacy, compared with DR and GrA (Fig. S8). Unlike the reproductive tissues, most of the highly expressed genes in seedlings showed functions in chlorophyll biosynthesis and photosynthesis (Fig. S8).

To assess the co-expression of DEGs during spike development, all DEGs from DR, AP, GrA, and MSA_1394 were combined. These 4888 DEGs were grouped into 28 distinct expression clusters, which were further grouped into three phase-specific CLUSTERs, I-III (Figs. 4-A, S9; Tables S11, S12).The DEGs in CLUSTER I were expressed predominantly in the late reproductive phase, GrA (Fig.4-A; Table S11).Most MADS genes with this expression pattern were placed into the same subcluster, cluster 11 (Fig. 4-B). CLUSTER I was enriched for photosynthesis- and transport-related GO terms, which overlapped with GO terms of the DEGs highly expressed in seedlings(SE_3731) (Figs.5-B,C, S10-A;Table S13),suggesting active assimilation requirements and nutrient exchanges during the GrA phase. DEGs in CLUSTER II showed slightly higher expression levels at AP and DR than at GrA (Fig. 4-A;Table S11), with enrichment for genes with GO terms of cell behavior and cellular component organization.Notably,34 of 36 GO terms in CLUSTER II were the same as those of MSAs_1394 (Figs. 5-A, D, S10-B), supporting the consistent expression pattern of MSAs_1394 in reproductive stages.CLUSTER III, containing genes preferentially expressed at the DR phase (Fig. 4-A; Table S11), was enriched in genes associated with gene regulation and plant growth, whose GO terms were divergent from those of MSAs_1394 and SE_3731(Figs. 5-E, S10-C; Table S13), indicating a unique pattern of gene expression during barley early reproductive initiation.Thus, untargeted co-expression analysis and GO enrichment at defined developmental stages of barley MSA show that the major regulators of barley spike development are orchestrated in a stage-specific pattern.

3.7. Cell cycle activity is high in the AP stage of the barley inflorescence

Among DEGs in the developing MSA, a large number of expressed genes are involved in cell division and cell cycle, such as cell division control protein 48 (CDC48), cyclins and histones (Fig. 6-A, D, H; Table S8), likely required by the initiation progress and growth of spike and spikelet organs. Generally, the cell cycle comprises four phases: G0/G1 phase (cell growth), S phase (DNA synthesis), G2 phase (cell growth), and M phase (nuclear and cell division) [39]. Regulatory proteins (cyclins) are likely to be essential for all phase transitions and cyclin protein levels fluctuate throughout the cell cycle. Most genes of the cyclin family showed the highest expression level in the AP stage compared with other stages of the MSA (Fig. 6-A, E).

Histone genes,encoding a large family of ubiquitous proteins that are form the basic structural unit of chromatin and are involved in epigenetic regulation [39], were highly expressed at the AP stage(Figs.6-B,C,G;S11-A).Several other genes involved in cell cycle regulation, such as histone-lysine Nmethyltransferases and a type B histone acetyltransferase(HAT1) (Fig. 6-D, G), B-type cyclin-dependent kinases (CDKBs,Figs.6-F,S11-B),E2F transcription factors(Fig.S11-C)and CDC48 homologs(Figs.6-H,S11-D),also showed peak expression at the AP stage.Expression of histone and cyclin genes were assembled to co-expression cluster 10(Fig.4-A,B),a finding consistent with the dominant GO terms of cell proliferation at AP stage(Fig.5-D).Together, the synchronized expression profiling of cell-cycle genes indicated relatively high cell cycle activity in the AP phase during barley spike development.

3.8.Control of the six-rowed phenotype by VRS genes occurs at the transcriptional level

Determination of the two- versus six-rowed phenotype in barley spike is regulated by a series of VRS genes [3,4,8,40,41].We compared transcriptomes between the DR (W2.0) and AP(W3.5) stages in the two-rowed cultivar Scarlett and the sixrowed cultivar Morex.Using the raw RNA-seq data of Scarlett(Table S4)[18],we identified 703 DEGs in Morex and 1270 DEGs in Scarlett (Figs. 1-A, 7-A; Table S4). Of these, 200 DEGs were common to Morex and Scarlett,including some key transcription factors such as MADSs and bHLHs (Fig. 7-A; Table S4).Barley spike pattern regulatory genes VRS1, 2, 3, 4 and INT-C(VRS5),showed different expression levels between Morex and Scarlett in the DR to AP stages, with the expression of some genes declining in Morex and increasing in Scarlett (Fig. 7-B,C).The expression levels of VRS1 and INT-C were significantly higher in Scarlett MSAs than in Morex, confirming the essential role of VRS genes in repressing barley lateral spikelet fertility and determining barley inflorescence patterning (Fig.7-D).

3.9. Expression of genes controlling plant hormone levels changes dramatically during barley inflorescence development

Because the distribution and regulation of plant hormones play key roles in specifying inflorescence meristem formation and patterning [1,8], we investigated key genes in the auxin,GA and CK pathways (Figs. 8-A-C, S12). Auxin biosynthetic and distributive genes, Tryptophan Aminotransferase of Arabidopsis (TAAs) and YABBYs (YABs), generally showed increased expression from early reproductive stages (DR, AP)to late reproductive stage (GrA) (Figs. 8-A, S12-A). RNA-seq results also revealed the dynamic expression of auxin signaling repressors (AUX/IAAs), regulators (SAURs) and transporters, and the decreasing expression trends of most auxin response factors(ARFs)during inflorescence development(Fig.8-A).

Fig.4- Co-expression clustering of DEGs.(A)Heat map of co-expression clusters for 4888 DEGs.Colors represent log2-fold changes(log2-FC)in expression levels relative to the mean transcript abundance.Co-expression clusters(1-28)were assigned to three higher-level groups(CLUSTER I-III)following tissue-specific expression patterns.The number and assignment of DEGs to clusters are shown above and below the heat map.The similarity of co-expression clusters is indicated in the hierarchical tree structure above the heat map.(B)Representative co-expression clusters of DEGs during barley inflorescence development.The cluster sizes and example members are indicated above the co-expression plots.The expression levels of individual transcripts (gray lines)and the mean expression level across all transcripts within each cluster(red line)are plotted.The coexpression plots are shown as mean centered and scaled transcript levels(Z-scores).

Fig.5-Comparisons of Gene Ontology(GO)classifications of DEGs among MSAs and seedling.(A,B)Overrepresented GO terms assigned to biological processes,cellular components and molecular function for genes over expressed in MSAs(MSA_1394)(A)or in seedlings (SE_3731)(B).(C,D, E)Overrepresented GO terms assigned to biological processes,cellular components and molecular function for genes in CLUSTER I(C),CLUSTER II(D),and CLUSTERIII(E).

GA2OX,encoding a gibberellin 2-oxidase that catalyzes the deactivation of bioactive GA, was highly expressed at DR and AP stages, although biosynthetic genes GA20OX3 (Gibberellin 20 oxidase 3) and GA3OX1 (Gibberellin 3-oxidase 1) were highly expressed at AP and GrA phases, respectively (Figs. 8-B, S12-B),suggesting that higher gibberellin levels are present at later stages of MSA development. Consequently, gibberellin response factors (GASA) and receptor genes (GIBBERELLININSENSITIVE DWARF1, GID1) were expressed more highly at later stages(Fig.8-B).

Cytokinin biosynthetic genes (Isopentenyltransferase, IPT and Lonely Guy,LOG)also increased in expression levels over time (Figs. 8-C-E, S12-C), while cytokinin oxidase/dehydrogenase genes (CKXs) and type A cytokinin response regulators (ARRs) tended to peak at the AP stage (Fig. 8-C),suggesting active CK signaling in barley inflorescence development. IPT and LOGs were also expressed at higher levels and showed more significant expression changes in the developing MSA of Morex than in that of Scarlett (Fig. 8-D, E), implying that the CK metabolism or signaling is more active in six-rowed than in two-rowed barley. Thus, our results showed that plant hormones, including auxin,gibberellin and cytokinin, are active during barley spike development, indicating an essential role in barley inflorescence initiation, floral meristem differentiation and floral organ growth.

Fig.6- Cell cycle activity in the barley inflorescence. (A-D)RNA-seq expression of cell cycle-related genes differentially expressed in the developing barley inflorescence and seedlings.The normalized expression values of Cyclins(A),Histone H2A/2Bs(B),Histone H4s(C)and Histone-lysine N-methyltransferases and Histone Acetyltransferase 1(HAT1)(D)are reported in fragments per kilobase per million(FPKM).(E-G)qRT-PCR validation of cell cycle-associated genes in the developing barley inflorescence and seedling.The transcript levels of Cyclins(E), Cyclin-dependent kinase B (CDKB)and E2F(F),Histones and Histone Acetyltransferase 1(HAT1)(G)are shown relative to the transcript abundance of HvActin7.(H)qRT-PCR expression data of genes encoding cell division control protein 48(CDC48E)are shown relative to the transcript abundance of HvActin7.All qRTPCR results are shown as mean ± SD of three biological replicates.DR,double ridge;AP,awn primordium;GrA,green anther.

4. Discussion

Remarkable architectural diversity exists among flowerbearing inflorescences in cereals that ultimately produce grain for food and feed. Patterns of inflorescence meristems arise from positional cues and their formation and development contribute directly to grain yield [1,15]. The molecular networks that control inflorescence development in the economically important crop barley remain largely unexplored.A previous study [18] showed the effects of the daylength and photoperiod response gene Ppd-H1 on the inflorescence/floral transition in a two-rowed barley variety,but the key regulators and potential novel factors that control the pivotal phases of inflorescence development have not been thoroughly investigated. We conducted a systematic study of molecular regulatory networks for spike development in the six-rowed barley cultivar (Morex). The transcription profiles of three key developmental stages reveal phase-specific gene expression patterns, the identities of novel, highly-expressed genes specifically expressed in MSAs, and the dynamic progression of cell behaviors and plant hormone signaling.Gene expression network and co-expression assembly indicate the conserved and divergent functions of barley homologs of known developmental regulators from other plant species. These findings shed light on the molecular mechanisms of barley inflorescence development and for yield improvement.

Table 2-Protein-coding genes highly expressed at MSAs(FPKM ＞150) but undetectable in seedlings.

4.1. Key regulatory genes are expressed in phase-specific patterns during barley inflorescence development

Unlike reproductive meristems of Arabidopsis, in which the inflorescence and flower meristems clearly differ, the spike and spikelet/floret in grasses are determined by a more complex hierarchy of specialized meristems[1].Gene expression profiles suggest that some key regulatory genes for barley spikelet meristem development have conserved functions with homologs in other cereals. The homologs of WUS in barley show the highest expression level at DR stage, and a decreased expression trend at AP and GrA stages, consistent with those in early developmental stages of the shoot apical meristem (SAM) in Arabidopsis and MSAs in rice (Fig. 2) [1,6].The high expression of FEA3-like and CLV1-like genes at AP and GrA stages, and their expression levels relative to WUS,support a conserved function for these genes in regulating meristem size and activity as in other plants (Fig. 2) [6,33]. In rice,maize,and Arabidopsis,LFY1,KN1,and their homologs are involved in floral meristem identity and maintenance[1,7,28,42], and high expression of two barley homologs at the DR stage suggests that their functions are also conserved in barley inflorescence initiation (Fig. 2). As in wheat [36],barley LFY is expressed at lower levels in AP and GrA stages after peaking at DR,but rice LFY is not likely to be expressed at later reproductive stages [41], suggesting that LFY may have diverse functions in regulating grass floret meristem differentiation and floral organ growth. Also, the DR-specific expression pattern of transcription factors such as bHLHs and AP2s reveals the dominant transcriptional regulation events in the early reproductive stage of barley.

During MSA development, the cell cycle and cell division play a critical role in the differentiation of individual floral organs and spikelet number, which directly link reproductive activity with crop yield traits. For example, mRNA of Arabidopsis D-type Cyclin, CYCD3, is abundant in the inflorescence;overexpressing CYCD3 leads to the disturbed cell cycle of meristem cells and reduced inflorescence and flower number [38]. At the AP stage of barley MSA, genes encoding cell cycle-related regulators,cyclins,histones,CDKs,and E2Fs are preferentially expressed, indicating high levels of cell proliferation and differentiation activity to enable floral organ formation (Figs. 6, S11). At the GrA stage, DEGs are enriched with photosynthetic and cell wall organization components,showing a large overlap with DEGs of SE_3731 (cluster 1 and 19). Thus, GrA-specific gene expression reveals active organ growth in the late reproductive stage of barley.

MADS-box proteins play a major role in regulating floral growth and development. Most barley MADS homologs,including B-, C-, and D-class genes, displayed meristemstaged expression preference, consistent with that of the genes governing rice inflorescence architecture and Arabidopsis floret determination [10,30], implying conserved roles for many MADS genes in specifying inflorescence meristems of barley. However, there were some distinct differences. In rice, several regulators such as the SVP-like MADS-box genes(OsMADS22,OsMADS47,and OsMADS55),and E-class MADS-box gene (OsMADS34) control inflorescence branching and regulate spikelet meristem identity, which act as key players in yield improvement [10,11,30]. In Arabidopsis and tomato, the E-class MADS gene SEPALLATA 4 (SEP4)regulates inflorescence architecture [30,43]. We found that the barley E-class and SVP-like MADS-box genes show expression patterns similar to those of homologs in rice and wheat (Fig. 3) [11,36].However, ectopic expression of SVP-like genes in barley inhibits spike development and causes floral reversion rather than changes in inflorescence branching[13].Despite the conserved gene expression, the roles of barley Eclass and SVP-like MADS genes in spike branching require further elucidation, as these homologs likely have differing functions in grasses with different inflorescence structures.

Fig.7-Expression of spike patterning determinants in 6-rowed and 2-rowed barley inflorescence meristems.DR,double ridge;AP,awn primordium;GrA,green anther.(A)Venn diagram showing the overlap of transcripts that are differentially expressed between DR(W2.0)and AP(W3.5)stages in 6-rowed(Morex)and 2-rowed(Scarlett)barley.(B, C)RNA-seq expression data of selected spike patterning-associated genes differentially expressed in the developing inflorescences of 6-rowed(B)and 2-rowed(C)varieties.The normalized expression values of VRS1,VRS2,VRS3,VRS4,and INT-C(VRS5)are reported in fragments per kilobase per million(FPKM).(D)qRT-PCR expression of VRS1,VRS2,VRS3,VRS4,and INT-C in three MSA stages of 6-rowed(Morex)and 2-rowed(Scarlett)barley.Error bars mean ± SD of three biological replicates.

4.2.Novel high-abundance genes play potential roles in barley inflorescence identity

The developmental events and properties in vegetative and reproductive tissues are markedly different. DEGs expressed preferentially in the inflorescence, MSA_1394, displayed generally stable expression levels throughout spike development and showed significant difference in GO enrichment in comparison with genes highly expressed in vegetative tissues(SE_3731). Several regulators highly expressed in the inflorescence were identified. For example, HvMADS32, a monocotspecific MIKCc-type gene, is the most highly expressed of all barley MADS-box genes at any stage of MSA development. In rice, however, OsMADS32 (Chimeric Floral Organs1, CFO1) is expressed only in early stages of inflorescence and late kernel development, and loss of function of OsMADS32 leads to defective flower organs[37].Thus,HvMADS32 is likely to have different regulatory functions than its rice homolog, with essential roles throughout floret development.

SPL family transcription factors are also regulators of reproductive development and inflorescence architecture in plants [2,44]. SPL4-like and SPL17-like genes are highly expressed in barley MSAs, peaking in the AP and DR stages,respectively (Tables 2, S10). However, rice OsSPL17 is highly expressed in the early reproductive meristem (rachis and primary branch meristems),and only weakly expressed in the spikelet meristem;downregulation of OsSPL17 greatly reduces panicle branch numbers[45].In Arabidopsis,several SPLs have been reported [43] to control floral transition by binding directly to, and activating the transcription of, LFY and MADS-box genes SOC1 and AP1. Thus, the highly expressed barley SPLs may also be involved in transcriptional regulation.

Fig.8- Genes regulating auxin,gibberellin, and cytokinin biosynthesis and signaling during barley spike development.DR,double ridge;AP,awn primordium;GrA,green anther.(A,B,C)Relative expression of genes regulating auxin(A),gibberellin(B),and cytokinin(C)in three comparison sets of reproductive tissue.Red represents overexpression and blue represents underexpression of genes in the first-named relative to the second-named tissue in each set.(D)qRT-PCR validation of cytokinin biosynthesis IPT gene expression in three MSA stages of 6-rowed(Morex)and 2-rowed(Scarlett)barley.Results are shown as mean ± SD of three biological replicates. (E)Expression changes of genes encoding cytokinin-activating enzymes LOG3,LOG4,LOG7 between DR and AP of 6-rowed(Morex)and 2-rowed(Scarlett)varieties.The value represents fold changes of transformed FPKM counts on the log2 scale.

Other strongly MSA-enriched genes may be involved in cell behaviors (Table 2). For example, the GATA transcription factor 15 (GATA15) regulates cell differentiation [46], and cell division control 48-E-like genes (CDC48s) and E3 ubiquitin ligase gene (ORTH2) contribute to the activity of cell division[38,47]; Argonaute (AGO4), Nucleolar protein 58, and RNAbinding protein 1 (GRP4) are essential for chromatin events and DNA/RNA binding [48,49]; and pectin lyase-like gene and Expansin B2 (EXPB2) play key roles in cell wall biosynthesis and organization [50,51]. Two other genes, a gibberellin (GA)-regulated gene (HORVU2Hr1G081190) and bHLH093, are involved in plant hormone-mediated development [52]. Thus,these meristem-enriched expressed genes likely regulate key processes of cellular events during barley spike development,and represent potential new regulators for further functional analysis.

4.3. Genetic regulators and phytohormone levels determine barley spike patterning

The transition of the barley inflorescence from bract primordium to floral organ differentiation is the key step in spikelet formation [15,53], which is accompanied by the downregulation of meristem activity genes and a pronounced upregulation of genes associated with organ development. Barley homologs of various MADS-box genes may play central roles in this transition phase. Most barley MADS genes, such as ABCDE-class genes, are expressed at early MSAs stage and increase toward later stages of floral development, and coexpression analysis revealed the clustering of many MADS genes (Fig. 4). The expression patterns of these MADS genes are conserved in rice and wheat [36,44], suggesting essential roles of MADS genes in barley floral organ differentiation and spike morphology.

The row type of spike patterning is controlled by at least five identified loci, VRS1-5 [3,4,8,39,40]. Although van Esse et al.[54]found that the expression level of VRS1 plays a central role in determining the row type of barley,our results showed that the expression patterns of VRS1-5 are all different between six-rowed and two-rowed barley. Given that both VRS1 and VRS5(INT-C)are highly expressed at the AP and GrA stages (Fig. 7), transcriptional control of these two genes is a mechanism that likely regulates and represses barley lateral spikelet fertility.

In addition to the regulatory activities of various transcription factors, phytohormones are also involved in barley MSA formation and growth, with elevated expression of auxin-,GA-,and CK-associated biosynthetic and signaling genes from the MSA initial stage to floret organ differentiation stage.Notably,these genes are also regulated by other factors in the MSA.For example,the maize KN1 gene negatively controls the amount of GA in reproductive meristems by directly inducing the expression of the GA-inactivating enzyme gene GA2OX1.Rice and Arabidopsis KNOX (KN) genes promote the accumulation of cytokinin by inducing the cytokinin biosynthetic gene IPT [1]. In barley, the various expression patterns of the KN-like gene family are also likely to play a central role in the coordination of hormone signaling (Fig. 2). The expression patterns of plant hormone biosynthetic genes suggest variation in phytohormone levels between two-rowed and sixrowed barley.

5. Conclusions

Our findings reveal phase-specific gene expression patterns during barley inflorescence development,suggest new candidates for research, and provide resources for barley MSA analysis. When these findings are combined with technical advances such as genome editing by CRISPR/Cas9 and efficient transformation systems, we will be able to better dissect the regulatory pathways in barley inflorescence development. Such work should also provide valuable candidate genes that can be manipulated to improve the yield of barley.

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2019.04.005.

Declaration of Competing Interest

Authors declare that there are no conflicts of interest.

Acknowledgements

We thank Yang Zhang (King Abdullah University of Science and Technology, Saudi Arabia) for collecting materials, Jie Zong (Shanghai Jiao Tong University, China) for providing guidance in genome analysis, and Natalie Betts (The University of Adelaide)for valuable comments and for revision of the article. This work was supported by Australian Research Council (DP170103352); an Australia-China Science and Research Fund Joint Research Centre grant ACSRF48187;Start-up funding (Australia, 13114779, 62117250) for Dabing Zhang from the School of Agriculture,Food and Wine,The University of Adelaide; and the Innovative Research Team, Ministry of Education of China; the 111 Project(B14016).