A barley spike transcriptome atlas by tissue-specific RNA-seq of LCM tissues

To identify transcriptional signatures and molecular circuitries associated with meristem fate and organ differentiation along the spike axis, we used LCM to isolate IM, LR, and SR at the DR stage as well as IM, CS, and LS meristems at five initial stages of spike development ranging from triple mound (TM) to awn primordium (AP) (Fig. 1, A and B) from two-rowed barley cv. Bowman. In addition, IM at white anther (WA) was isolated, while the rachis (R1, provascular tissues) and complete spike sections (R2) from each developmental stage were selected as reference samples. Meristems of vegetative tissues, the vegetative apex (VA) containing the shoot apical meristem (SAM), the root tip (RT) including the root apical meristem (RAM), and the proliferating leaf blade (LB) from the primary leaf were collected as control samples for spike meristems (see list of enriched genes compared to R2 across all stages; table S1).

To interrogate changes induced upon carbon starvation of the plant, identical tissues (except vegetative, R1, and R2 tissues) were also isolated from shaded plants (∼70% light reduction). In total, four biological replicates per sample (with the exception of triplicate samples for VA, RAM, and LB) resulted in 209 highly tissue-specific and informative samples for RNA sequencing (RNA-seq). Merging all expressed genes from the six individual spike domains [IM, CS, LS (LR, SR), and R1] increased the number of detected genes to 34,540 [20,768 high confidence (HC) + 13,772 low confidence (LC)] compared to R2 alone [29,310; 18,610 HC + 10,700 LC; transcript per million (TPM) ≥ 1] (Fig. 1C). The increase of ~17%, i.e., 5230 more genes, shows the benefit of tissue-specific sampling in capturing rare transcripts that are otherwise not detectable in the entire organ samples (i.e., below threshold TPM ≥ 1 for expression). Principal component analysis (PCA) of all 134 samples under control conditions (Fig. 1D) shows a clear separation between spike and vegetative meristems (RT, LB), indicating distinct transcript profiles activated in spike tissues. The tissues R1 and IM at the WA stage clustered separately from the other spike tissues, revealing molecular programs probably associated with (pro-)vascular functions or tissue abortion, respectively. The PCA for spike tissues depicts predominant clustering according to tissue types, with largely distinct transcript profiles of IM, SR, and LR compared to CS and LS and broad separation of early and late spike developmental stages in the second component. Expression of the top 50 genes specifically up- or down-regulated in IM, CS, and LS across development demonstrates that the dataset is a highly valuable resource for extraction of domain-specific genes (Fig. 1E and table S2). Whole-mount in situ hybridization (WISH) on meristems at the SP/AP stage confirmed tissue specificity of candidate genes (Fig. 2). In silico expression profiles displayed by the tissue ePlant viewer for IM-specific homeobox-leucine zipper family/HvHD-zip (HORVU2Hr1G095210; Fig. 2A), CS-specific (late stages) HvMADS7 (HORVU7Hr1G054220; Fig. 2B), and LS-specific VRS5/HvTB1 (HORVU4Hr1G007040; Fig. 2C) coincide with signals from WISH.

The shading treatment had a notable low impact on transcriptional profiles with only marginal differences between growth conditions (except IM at WA) (fig. S1), and the low number of differentially expressed genes (DEGs) across common developmental stages (149 DEGs in at least two stages; table S3). Therefore, we based our analyses only on control conditions reflecting normal spike development.

Vegetative apical meristem versus IM—Transcriptional signatures underlie SAM maintenance and induction of flowering

The conversion of SAM into IM represents the beginning of reproductive development. Here, we used the dataset for a comparison of VA (including SAM and adjacent LRs) and IM at the DR stage, two distinct meristematic tissues that are usually highly underrepresented in whole-organ samples because of diluting effects caused by their tiny size. In brief, morphometric changes during conversion of SAM into IM are associated with massive transcriptional reprogramming indicated by 1577 DEGs between VA and IM [log₂ fold change (FC) > 1, false discovery rate (FDR)–adjusted P < 0.05; table S4]. Gene Ontology (GO) term enrichment in the 1324 genes up-regulated in VA (Fig. 3A, I) suggests that processes related to adaxial/abaxial specification, axis polarity, meristem growth/maintenance, and circadian rhythm are activated (terms related to axis specification might also be attributed to the remaining LR activity). In contrast, photosynthesis, carbon fixation, light responses, as well as regulation of developmental processes and flower development stand out as enriched processes in IM (table S5). Intriguingly, transcripts involved in stress responses, such as senescence-associated protein (SAG), several jasmonate (JA)–/aluminum-induced proteins, and OXIDATIVE STRESS 3 (OXS3), belong to the most strongly up-regulated genes in VA (table S4), which coincides with expression of distinct transcription factor (TF) families (Fig. 3A, II), like AP2/ERFs, MYBs, bZIPs, and WRKYs, that are associated with defense responses and ethylene, abscisic acid (ABA), salicylic acid (SA), and JA signaling (15), thereby overlapping with accumulated hormone signaling elements in mature maize SAMs (16). Other TFs activated in VA belong to the AUX/IAA, bHLH, zinc finger, YABBY, GRAS, homeobox (HB), and PLATZ families. Five genes encoding YABBY proteins (YABBY3, YABBY4, and YABBY5), including an Arabidopsis thaliana (Arabidopsis) CRABS CLAW ortholog, are highly up-regulated in VA. C2C2-YABBY proteins are plant-specific TFs that are required for polarity specification in lateral organs (17) but are also proposed to restrict meristematic activity and growth to SAM by repression of KNOTTED homeobox genes (18). Up-regulated barley homeobox TFs, i.e., class I KNOX1-like, BEL1-like, WOX1, and PROTODERMAL FACTOR 2 (PDF2) (table S4) (19–22), indicate a regulatory role in SAM organization and maintenance. GRAS family TF (HORVU1Hr1G053510) shows a high similarity to LOST MERISTEM (LOM1+2) genes that were shown to promote cell differentiation at the periphery of Arabidopsis shoot meristems and are required for polar SAM organization (23). Three isoforms encoding SHORT VEGETATIVE PHASE MADS-box proteins and TERMINAL FLOWER 1 (TFL1/HORVU2Hr1G072750), previously shown to act as floral repressors (24–27), are concomitantly up-regulated in VA. Although only a few TFs are up-regulated in IM, some of these are indicative of initiated reproductive development. Three MADS-box genes, VRN-H1/HvBM5A, an APETALA1/FRUITFULL-like floral meristem identity gene orthologous to Arabidopsis APETALA 1 (AP1), and two SEPALLATA (SEP) orthologs (HvMADS 34/-WM19B) potentially involved in meristem growth and floral development (28–30) are distinctly induced in IM. Other TFs predominantly expressed in IM encode barley orthologs of LEAFY COTYLEDON 1 (LEC1-like), bHLH APRATAXIN (APTX)–like, and FLOWERING PROMOTING FACTOR 1 (FPF1), which modulate floral competency in Arabidopsis (31). A detailed summary of this tissue comparison can be obtained from table S4. Together, transcriptional profiles in VA show activation of factors involved in stem cell specification/SAM maintenance and prevention of flower induction, whereas the transition of SAM into an IM is characterized by genes that are described to promote reproductive development and flowering, particularly AP1/SEP orthologs.

SR versus LR—Distinct transcriptional signatures reflect opposing meristematic fates during early reproductive differentiation

SRs are formed as single axillary meristems that later partition into one CS and two LS—the spikelet triplet in barley. To identify regulatory pathways responsible for spikelet initiation and LR suppression, we compared the transcriptomes of SR and LR at the DR stage. Only 64 transcripts were differentially expressed between SR and LR (log₂ FC > 1, FDR-adjusted P < 0.05; table S6), enabling us to investigate fine-tuned differences in the transcriptional control of these juxtaposed tissues.

Despite the small number of DEGs, GO term enrichment shows that distinct biological processes are activated in each ridge type (Fig. 3B, I, and table S7). Regulation of transcription, RNA metabolism, gene expression, and cellular macromolecule synthesis are highly enriched in the SR. Opposingly, processes related to adaxial-abaxial specification (axis polarity), specification of organ identity, chromatin arrangement, and response to ABA are activated in the LR. Axis specification is one of the key features during leaf formation for optimizing photosynthetic capacity and confirms leaf identity of LR cells. Other regulatory genes are deemed to be involved in defense responses, senescence, and tissue abortion (Fig. 3B, I). The barley THIRD OUTER GLUME (TRD) gene, encoding a GATA TF involved in the suppression of LR development (13, 32, 33), showed a higher expression trend in the LR, signifying its role in LR suppression (fig. S2). A similar role could be attributed to several SQUAMOSA PROMOTER BINDING-LIKE (SPL) TFs characterized for vegetative to reproductive phase transition and leaf development such as orthologs of maize UNBRANCHED 2 (rice IDEAL PLANT ARCHITECTURE 1; HvSPL14/IPA1), TASSEL SHEATH 4, and TEOSINTE GLUME ARCHITECTURE 1 (fig. S3A) (34, 35). Our RNA-guided endonuclease knockout (RGEN KOs) lines of HvSPL14/IPA1 showed an extended vegetative phase and prolonged LR development with elongated LRs at DR and, consequently, a significant increase in total leaf number on the main culm (fig. S3, B and C).

In contrast, SRs show transcriptional activation of regulatory genes that are associated with reproductive development (Fig. 3B, II). TFs, such as WD40 domain–containing transducin, SIX-ROWED SPIKE 2 (VRS2), homeodomain protein 40 (HB40), and bHLH APTX, may play a role in boundary and floral formation. VRS2 is required for floral organ patterning in barley, probably by directing auxin and gibberellin homeostasis along the spike axis (8). Orthologs of LITTLE zipper genes (ZPR3 and ZPR4) are highly expressed in the SR. ZPRs competitively inhibit HD-ZIPIII TFs involved in Arabidopsis SAM development (36). Intriguingly, several genes involved in auxin signaling are up-regulated in the SR, including HvAUX/IAA19, an ortholog of Arabidopsis IAA8, ARGONAUTE1, an RNA-slicer that regulates expression of ARF17, and orthologs of Arabidopsis ACC synthase11 and HB40 that were shown to respond to applied auxin. Other transcripts are associated with cytokinin metabolism/signaling: genes encoding two LONELY GUY enzymes and a type A response regulator (ARR6). Together, this is the first detailed global analysis of differential gene expression between SR and LR in grasses, thereby clearly reflecting opposing developmental fates of the two tissues. Our analysis showed that transcriptional signatures of axis specification that are typical for leaf development and cellular disintegration are characteristic for LRs, whereas auxin- and cytokinin-related genes may be involved in cell specification (cytokinesis, cell elongation/polarity) of SRs. The DR stage or a stage with a similar development (a diminutive leaf meristem subtending SM) is present in the vast majority of grass species (37, 38), reinforcing the broad applicability of these findings to other grasses.

CS versus LS fates—The gradual and multileveled shutdown of lateral floral organs

The fertility status of the two LSs determines the barley row-type (fig. S4), wherein the LSs of six-rowed barley are fertile and set grains in contrast to sterile LSs in two-rowed barley. At least five genes control the inhibitory nature of LSs in two-rowed spikes (negative regulators of LS development), including VRS1, VRS2, VRS3, VRS4, and VRS5/INTERMEDIUM-SPIKE C (INT-C) (2–5, 7, 8). The sterility in LSs is predominantly achieved through pistil and anther abortion (3). Although LSs are destined to be sterile, other spikelet and floret organs such as glumes, lemma, palea, and lodicules are still formed but reduced in size (figs. S4 and S5E). Quantitative variation in LS size indicates the presence of unknown regulators of LS development (39, 40).

For a comprehensive understanding of the regulatory landscape involved in LS fate, we analyzed the transcriptome of CS (destined to become fertile) and LS (destined to be inhibited) tissues at five developmental stages starting from TM until AP (Fig. 1A). A time course DEG analysis between CS and LS tissues at these stages identified in total 2047 DEGs (1225 HC and 822 LC; log₂ FC > 0.5, FDR-adjusted P < 0.05; Fig. 4A and table S8). Our data showed that significant transcriptomic changes between CS and LS start to appear from the GP stage that marks a clear separation between these tissues (Fig. 1A and fig. S6). It should be noted that in cv. Bowman LSs show a developmental delay compared to the CS, which is also reflected in the expression pattern of known floral homeotic genes and may potentially inflate the number of DEGs between CS and LS (fig. S5 and Supplementary Notes). Among developmental stages, LP and SP showed the highest number of DEGs (809 and 814, respectively; Fig. 4, A and B), revealing the main shift in gene expression indicative of a developmental retardation in LSs at these stages.

K-means clustering of DEGs over five CS and LS stages (HC genes) defined six coexpression clusters (C1 to C6). Two CS-specific clusters showed elevated expression from TM to SP (C1) and LP to AP (C2), whereas four clusters presented stage-specific or overall transcript up-regulation in LS (C3-GP; C4-GP to AP; C5-TM to AP; C6-SP) (Fig. 4C).

Genes in the C1 cluster (CS-specific) are enriched for processes involved in vascular development, auxin signaling and homeostasis, ethylene response, and two-component signaling (TCS) phosphorelays (Fig. 4D). CS-specific C2 cluster is predominantly enriched for processes involved in cell differentiation, reproductive organ development, response to light stimulus, and transcriptional regulation/TF binding (Fig. 4D). The C2 cluster also contained 15 MADS-box genes, including orthologs of various floral homeotic genes (SEP3, AGL6, AGL14, MADS1, MADS7, and MADS58 isoforms) particularly those involved in floral organ specification [(41); table S9], showing an elevated expression from LP to AP when floral organs are formed (see also Fig. 2B). Other enriched TFs belong to TCP and MYB families, while a suite of RAS–guanosine triphosphatase (GTPase), RHO-GTPase, and cyclin genes suggest activated cytokinesis and vesicle trafficking pathways that might be associated with establishment of cell polarity/specification during CS differentiation.

The row-type genes VRS1, VRS2, VRS3, VRS4, and VRS5/INT-C (inhibitors of LS development) displayed characteristic higher expression in LSs (Figs. 2C and 4E). In addition, other genes proposed to regulate tillering, grain size, plant height, or grain number also showed up-regulation in LSs. These include DENSE AND ERECT PANICLE 1 (HvDEP1), UNICULM 4 (HvBOP1/2/CUL4–BLADE ON PETIOLE), and HvSPL14/IPA1 (table S8) (42–44). Along with SPL14, most of the SBP (SQUAMOSA promoter binding protein-like) family TFs are up-regulated in LS, implying their role during LS development (fig. S7 and Supplementary Notes). Notably, the LS up-regulated clusters are enriched for stress responses and growth repression processes, such as ABA response/signaling (C4-C6), GA (gibberellic acid) catabolism (GA2 oxidases, C5), JA response, and secondary metabolism (C6) (Fig. 4D and table S9), hinting at a prominent, yet undiscovered, role for ABA, JA, and GA in the multileveled abortion process of LSs (Supplementary Notes).

Identification of coexpression modules specifying fates of CS and LS

To identify gene coexpression modules, weighted gene co-expression network analysis (WGCNA) was conducted using the CS and LS expression data from TM to AP stages under control conditions. The WGCNA identified 16 coexpression modules [8227 genes; coefficient of variation (CV) ≤ 25%], showing CS- and LS-specific expression patterns (fig. S8 and table S10). Seven modules show specific higher expression in CS (CS-M1 to CS-M7), whereas eight modules are specifically up-regulated in LS (LS-M1 to LS-M8; Fig. 5A). Notably, most of the coexpressed genes belonged to LS-specific modules (6319 genes, 76.80%), further underlining the importance of tightly coordinated transcriptional modulation of LS developmental fate. Because one of our main interests is the understanding of the molecular processes underlying LS fate, we focused on modules preferentially up-regulated in LSs that were further classified into LS-up-early (peaking expression at GP; LS-M1 to LS-M3) and LS-up-late modules (peaking expression at AP; LS-M4 to LS-M8) (Fig. 5A).

The module LS-M1 harbors numerous genes involved in organ boundary formation/signaling and/or meristem maintenance (Fig. 5B). These include VRS4/HvRA2, HvSPL14/IPA1, HvSPL17, HvBOP1/2/CUL4, LATERAL ORGAN FUSION 1 (LOF1), AUXIN RESPONSE FACTOR 3 (ARF3; ETTIN), LIKE AMP 1 (LAMP1), HOTHEAD, GROWTH REGULATING FACTOR 9, FTSH PROTEASE 4, UDP-GLUCOSE DEHYDROGENASE 3, BRASSINOSTEROID KINASE 6, and the SCARECROW (SCR) genes, SHOOT GRAVITROPISM 1 and SCL8 (4, 42, 43, 45–50).

Organ boundary domains are characterized by low rates of cell division, proliferation, and asymmetric cell shape, probably directed by auxin/brassinosteroid activity (51). Notably, four hub genes from the LS-M1 module were implicated in cell proliferation control (HvDEP1) (52), asymmetric cell division (SCL8/GRAS) (53), negative regulation of auxin transport and cell growth [hydroxycinnamoyl CoA: shikimate hydroxycinnamoyl transferase (HCT)] (54), and membrane lipid catabolism (GDSL esterase/lipase) (55), pointing to a pivotal function in boundary formation between CS and LSs. Moreover, a cis-motif analysis of promoters from LS-M1 module genes shows a general enrichment of motifs bound by genes potentially regulating cell division, proliferation, and growth (cis-motifs identified from Arabidopsis/maize DNA affinity purification sequencing (DAP-Seq) or chromatin immunoprecipitation sequencing data; fig. S9). Spikes of our HvSPL14/IPA1 KO lines showed development of ectopic additional spikelets/florets (AS/AF; apart from the basal triple-spikelet) on the flanks of LSs (Fig. 6C, II), indicating that expression of these genes from TM-GP is essential for regular boundary definition between CS and LS. Similar boundary defective phenotypes are also evident in vrs4 mutants (another member of LS-M1 module; fig. S4). Along with members of the LS-M1 module (HvSPL14/IPA1 and VRS4), VRS3 showed similar expression maxima at TM and GP (Fig. 4E). vrs3 mutants also display ectopic AS/AF phenotypes (2, 9).

The expression patterns of selected genes involved in organ boundary formation/signaling and meristem maintenance in Arabidopsis, rice, maize, and barley (Supplementary Notes and table S11) revealed that most of the barley orthologs had expression maxima in LS, mainly at TM and GP (Fig. 5D and Supplementary Notes). Thus, a tight regulation of these genes at TM and GP stages seems to be decisive for the generation of three individual SMs derived from a single SR (at DR), a process potentially controlled by organ boundary genes.

LS-up-late modules (LS-M4 to LS-M8) show expression peak around AP (Fig. 5A), probably indicating their involvement in LS suppression. Module LS-M6 harbored key row-type genes (inhibitors of LS development) VRS1 (HvHOX1), its paralog HvHOX2, and VRS5/INT-C (HvTB1) (fig. S10). GO term enrichment in LS-up-late modules, in general, depicts activated cellular metabolic processes (nucleic acid/nitrogen metabolism) proposed to be associated with cellular energy status and resynthesis of cell components (figs. S11 to S14). The LS-M7 module is enriched for cell cycle control and cellular component biogenesis, revealing that growth of LS is not completely inhibited, although LSs are destined to become sterile (fig. S14 and table S12).

Genes in LS-M4, LS-M5, and LS-M6 modules are enriched for bZIP TF binding sites that are deemed to be involved in ABA signaling (fig. S15 and table S13). Up-regulation of ABA signaling genes in the LS may support a role for ABA as one hormonal trigger during the LS abortion process (figs. S16 and S17, table S14, and Supplementary Notes).

Novel regulators of LS fate

Our high-resolution dataset facilitated the identification of multiple genes with characteristic high expression in LSs typical for row-type genes (Fig. 4E). To further verify gene function during LS development, we performed in planta functional studies for hub genes from the LS-M1 module, i.e., HvDEP1, HvBOP1/2/CUL4, and HvSPL14/IPA1 (Fig. 5B and table S15). HvDEP1, a positive regulator of cell proliferation and organ size, has been implicated in plant height, grain size, and inflorescence architecture control (44, 52). Evaluation of LS morphometry in three barley dep1 mutant alleles—BW-ari-e.1, BW-ari-e.GP, and ari-e.119 (table S16)—elucidated an apparent reduction in LS size (Fig. 6A, I to III, and fig. S18, A and B), indicating that DEP1 positively affects cell proliferation and organ growth in LSs. Moreover, the dep1 mutants also displayed a reduction in CS area, spike length, and plant height, further corroborating the view that DEP1 can positively affect overall plant growth (fig. S19).

HvBOP1/2/CUL4 orthologs have previously been implicated in specifying polarity of the Arabidopsis leaf petiole (56), suppression of axillary branch outgrowth in maize (57), and specification of the LB/sheath boundary in barley and rice (43, 58). Analysis of two cul4 mutant alleles, BW-cul4.5 and cul4.24, exhibited a prominent increase in LS area (Fig. 6B, I to III, and fig. S18C), indicating that CUL4 negatively affects LS growth. Crucially, the florets of LSs of the cul4 mutants developed anthers with pollen, enlarged pistils (fig. S20), but similarly larger CS, further reinforcing the idea that CUL4 acts as an inhibitor of floral development (fig. S21, A and B). Notably, the barley LAXATUM-A (LAX-A; HvBOP1/2; HORVU5Hr1G043220) (59), a homolog of CUL4, is also preferentially expressed in LSs at TM/GP stages (Fig. 5D), similarly underlining a putative role during LS inhibition.

Barley spl14 mutants were generated using RGEN technology (figs. S22 and S23). Three independent spl14 mutants displayed quantitative variation for LS development ranging from enlarged LSs with awn elongation to complete fertility restoration and grain development (Fig. 6C and fig. S23), clearly demonstrating that HvSPL14 also functions as an inhibitor of LS development similar to row-type genes. The spl14 mutants showed defects in spikelet development with more than three spikelets/florets per rachis node (Fig. 6C, II), suggesting that defective organ boundaries between CS and LSs may promote allocation of cells to lateral domains, which, in turn, enhance LS growth and the development of AS/AF. This observation appears to be consistent with the known function of maize UNBRANCHED 2 and UNBRANCHED 3 (HvSPL14/IPA1 orthologs), limiting cell differentiation rate to lateral domains of the meristem (34).

Decisions about meristem fate in a plant’s life cycle determine fitness. This is of high importance for staple crop plants because, here, specific decisions on floral meristem fates finally decide about the number of grains and yield potential. Particularly in temperate cereals, such as wheat or barley, knowledge about molecular networks defining reproductive meristem identity and developmental fates remains scarce. To fill this gap, we generated a high-resolution transcriptome atlas of floral meristems in barley by using LCM-based RNA-seq. Expression data highlight key molecular switches associated with differential meristem and organ fate decisions. We show the transcriptional regulation of the transition of VA into IM and identified key factors that promote reproductive development or, in contrast, that are associated with SAM maintenance and prevention of floret induction. Tissue type–specific profiling identified molecular signatures that define opposing developmental fates of adjacent and juxtaposed primordia as exemplified for SR and LR, or LS and CS meristems, and have so far not been described in any grass species. The transcriptome of the SR is characterized by regulatory genes (TFs) probably giving rise to spike (reproductive) development and elucidated underlying hormonal regulation (i.e., auxin and cytokinin) determining SR initiation. On the basis of expression differences between CS and LS meristematic tissues and network modeling, we identified several genes affecting morphological traits and uncovered pathways, such as ABA signaling, homeostasis, and GA metabolism, which might contribute to LS growth and development. Functional studies of network-derived regulators confirm that transcriptome data enabled the discovery of novel regulators for LS morphometry and growth inhibition.

Our control samples of vegetative meristems and different growth conditions were also included in the dataset to provide broader applicability of this resource, allowing to perform multiple comparisons between tissue types on a highly resolved spatiotemporal level. In conclusion, our data provide an essential framework for barley inflorescence developmental studies but may similarly be indispensable for other comparative studies of monocot species in general. To support broad adoption of this dataset, these data are incorporated into an ePlant browser for user-friendly access of any gene of interest.

Plant material and growth conditions

For meristem sampling, barley plants (cv. Bowman) were grown in 20-cm (ø) pots (with three plants per pot) in a controlled phytochamber at 12°C/12-hour day and 8°C/12-hour night with relative humidity of 60 to 70% and light intensity of 270 μE. Plants under control conditions were grown over a period of 16 months with an approximate 4-month time interval between successive biological replicates. For IM at WA, VA, and LB tissues, three to four biological replicates were collected from a single experiment. The VA samples were collected under short day conditions, i.e., 12°C/8-hour day and 8°C/16-hour night. For sampling of the RAM, the grains were initially soaked in a petri dish with moist Whatman filter paper at 4°C for 2 days and transferred to room temperature, and 1.0-mm root apices were collected 5 days after germination (DAG). For sampling of LB, about 5.0-mm LB tissue from the first leaf (21 DAG) located exactly above the ligule region (junction between leaf sheath and LB) was collected. For the shading treatment, plants of the four biological replicates were grown under a shade net reducing light intensity to 30% (~90 μE) of control (270 μE) with all other conditions as described for control plants.

Tissue preparation and mRNA sequencing

All tissue samples used in this study were collected between 2 and 5 hours after the initiation of the light phase. While sampling, proper measures were taken not to expose the plants to ambient light. Sample processing for RNA-seq has been performed as described in (60). Briefly, barley spikes were hand-dissected at the corresponding developmental stages, immediately frozen in liquid nitrogen, and transferred to a cryostat (−20°C). Spikes were glued onto sample plates by O.C.T. medium, and serial cross sections (18 μm) of the complete spike were mounted on ribonuclease (RNase)–free polyethylene naphthalate (PEN) membrane slides (MMI, Eching, Germany). To distinguish CS and LS, a distinct orientation is needed and spikes were strictly sectioned from the dorsal/front view. PEN membrane slides were stored for 7 days at −20°C until completely dried. The MMI Laser Cell Cut system was used to capture distinct spike tissues in 0.5-ml Isolation Caps (MMI). Usually between 50 and 250 tissue regions from at least 10 individual spikes were pooled for each of the four biological replicates. RNA was isolated by using the Absolutely RNA Nanoprep Kit (Agilent), including deoxyribonuclease (DNase) I digestion, and mRNA was amplified by one round of T7-based in vitro transcription using the MessageAmp aRNA Kit (Invitrogen). Antisense RNA was used for library preparation with Illumina TruSeq RNA kit v2 (Illumina), and libraries were sequenced on HiSeq 2500 (Illumina) to generate 100–base pair (bp) paired-end reads. On average, more than 51 million reads with mean Q30 of 91.9% and Phred quality score (Q score) quality score (PF) 35.7 were generated per sample, yielding in total 1.05 tera bases of sequence information. Raw data of this study are deposited in the European Nucleotide Archive (ENA) under the study accession number PRJEB39672.

Gene expression analysis

RNA-seq reads were processed using CASSAVA v1.8 and trimmed as recommended in the manufacturer’s instructions. RNA-seq quantification was done with Kallisto (61) using the barley genome [International Barley Sequencing Consortium (IBSC), 2016], and raw read counts were normalized to TPM expression levels. DESeq2 (62) was used for the identification of DEGs according to log₂ FC > 1 or log₂ FC > 0.5 (CS-LS comparison) and Benjamini-Hochberg (63) FDR-adjusted P values (<0.05). For the comparison between control and shading conditions (table S3), the threshold was lowered (adjusted P < 0.1) to identify DEGs. Barley sequence information is publicly available in the Barlex genome explorer (https://apex.ipk-gatersleben.de/apex/f?p=284:10), and TPM values of all samples are available under http://dx.doi.org/10.5447/ipk/2020/35 as dataset S1. Heatmaps for gene expression data were generated using the online tool ClustVis (64).

K-means clustering and GO enrichment analysis

A total of 1225 HC genes with a log₂ FC of ≥0.5 (FDR-adjusted P < 0.05) between CS and LS were selected for cluster analysis. To find the optimal number of clusters, figure of merit was calculated using median TPM expression data with a maximum number of clusters at 50 and 1000 iterations. For K-means clustering, the number of clusters was set to six based on the figure of merit calculated from Genesis 1.8.1 (65). K-means clustering was conducted with maximum iterations and number of clustering runs to 1000 in each case in Genesis 1.8.1. A total of 1000 randomizations were performed to test dependency of variables using a P value of ≤0.0001.

GO term enrichment in DEGs and K-means cluster genes was performed with G-profiler (66) using Arabidopsis orthologs (TAIR 10; e < 1.0 × 10^−5) as background for enrichment. Significance of enrichment was estimated by Benjamini-Hochberg FDR correction and/or Fisher’s exact test. For GO enrichment of WGCNA coexpression modules, the Cytoscape application BiNGO (67) was used with Arabidopsis as background reference using FDR correction (<0.05) for multiple testing.

Coexpression and cis-motif enrichment analyses

TPM data (HC and LC genes) of CS and LS tissues from 44 samples were used for construction of coexpression network, module detection, and identification of hub genes from the modules. The CV value for each gene among different replicates of each tissue was calculated, and genes with high expression variability within replicates (CV ≥ 0.25) were excluded from downstream analysis. The network was constructed using R package WGCNA (68) with 8227 genes. For this, a Pearson correlation matrix was calculated between all pairs of selected genes from which adjacency matrix was constructed by raising the correlation matrix to the power of 18. Average linkage hierarchical clustering was then carried out for grouping of genes owing to highly similar coexpression patterns. The modules were identified using Dynamic Hybrid Tree cut algorithm (69), keeping the minimum module size to 50. The expression profile of each module was summarized by module eigen-gene defined as its first principal component, and modules with highly correlated (>0.75) eigen-genes were merged. In this way, a total of 16 modules were identified. Furthermore, modules with expression patterns similar to genes of interest were selected for identification of hub genes and promoter enrichment analysis. For hub gene identification, the “degree” and “betweenness centrality” statistics were calculated separately for all the genes from each module. These values were used to rank the genes, and the top 10 genes from each module were considered as hub genes. The networks were visualized using Cytoscape 3.7.1 (70).

For cis-motif identification and enrichment analysis, 1500-bp upstream sequences from the predicted start codon (ATG) of all module genes were used. The analysis was carried out using “findMotif.pl” program from HOMER suite (71) that simultaneously performs de novo and known motif identification and enrichment analysis. The program was run to identify motifs of 4- to 10-bp sizes. The identified motifs were searched against 1500-bp background sequence defined by Homer database (based on GC content of the target sequences) for enrichment, and binomial distribution was used to calculate P values.

Phylogeny analysis

Barley SPL and SnRK proteins were extracted from Barlex genome explorer (https://apex.ipk-gatersleben.de/apex/f?p=284:10) using Arabidopsis SPL domain and kinase domains of SnRK1, SnRK2, and SnRK3 as query (e = 1.0 × 10^−5). Only full-length proteins were considered for phylogeny analysis, and sequences were aligned using MUSCLE. Maximum likelihood tree with bootstrap support inferred from 1000 replicates was calculated with MEGA X (72). The evolutionary distances were computed using the Jones-Taylor-Thornton (JTT) matrix-based method. Initial tree(s) for the heuristic search was obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model and then selecting the topology with superior log likelihood value. A discrete gamma distribution was used to model evolutionary rate differences among sites [five categories (+G, parameter = 0.8935)]. All positions containing gaps and missing data were eliminated from the analysis.

Whole-mount in situ hybridization

Barley spikes (cv. Bowman) were harvested at stamen and/or AP stages, immediately frozen in liquid nitrogen, and stored at −80°C until processing. Three to five spikes were used for each reaction. The reactions were generally carried out in 1.5-ml Eppendorf tubes with 125-μm nylon mesh baskets according to (73) for careful changing of solutions.

For preparation of RNA probes for HvHD-zip (HORVU2Hr1G095210), HvMADS7 (HORVU7Hr1G054220), and VRS5 (HORVU4Hr1G007040), ~300-bp sequences in the 3′ untranslated region were amplified by polymerase chain reaction (PCR) using gene-specific primers containing T3 and T7 promoter sequences (table S17). After purification, fragments were labeled with digoxigenin (DIG) by in vitro transcription using T3 and T7 polymerase according to the manufacturer’s instructions (Roche Diagnostics, Mannheim, Germany). After purification of riboprobes, efficiency of DIG labeling was verified by dot blotting and probes (100 ng of RNA) were denatured and treated with an RNase inhibitor for hybridization.

Hybridization was modified after (74). Briefly, spikes were fixed in 4% FAA (formaldehyde, acetic acid and ethanol) for 1 hour at room temperature after initial vacuum infiltration (3 min), washed, rehydrated, and treated with proteinase K (40 μg) for 20 min. Prehybridization (30 min) and hybridization (16 hours) with denatured probes were carried out at 52°C using hybridization buffer containing 2× SSC, 50% formamide, 2.5 g of dextran sulfate, 25 mg of transfer RNA (tRNA), and 5 mg of salmon sperm DNA. After washing, unbound RNA was digested for 30 min at 37°C using RNase A (20 μg/ml) (DNase free, Roche). Immunological detection using DIG antibodies (1:1000 in blocking solution) coupled with alkaline phosphatase and staining procedure with nitro blue tetrazolium (NBT)/bromochloroindolyl phosphate (BCIP) was done at room temperature according to the manufacturer’s instructions (Roche). Spikes were stained for 16 hours in darkness, and digital analysis was performed using a VHX-100K digital microscope (Keyence, Neu-Isenburg, Germany) or Zeiss Axioskop (Zeiss, Jena, Germany) and Axio Vision 4.7.2 software.

SPL14-specific mutagenesis by using RGEN

KO lines of the barley gene SPL14 (HORVU0Hr1g020810) were generated by targeted mutagenesis using RNA-guided Cas9 endonuclease (75). The SPL14 coding sequence was examined for potential target areas using the DeskGen selection tool (76), and guide RNA candidates were further validated in silico using the WU-CRISPR platform (77) and RNAfold (78) to model and evaluate their secondary structure. The guide RNA 2-4 targets, a motif in the second exon of SPL14 and the corresponding sequence, were integrated into the generic vector pSH121, resulting in pCH9-SPL-g2-4. The activity of guide RNA 2-4 was investigated using a yellow fluorescent protein restoration assay based on ballistic DNA transfer into barley leaf epidermis as previously described (79). The construct was transfected into Agrobacterium tumefaciens strain AGL1, and transformation was conducted using the barley cv. “Golden Promise” as previously described (80). The target region of SPL14 was amplified by PCR for Sanger sequencing, and the presence of the guide RNA and Cas9 transgenes in the plants was investigated (table S17). Progeny of three selected primary mutated plants was cultivated under greenhouse conditions (16 hours/18°C day, 8 hours/16°C night), and SPL14 modifications and segregation of the transgenes were confirmed in the next generation. Potential off-targets in the barley genes SPL17 (HORVU5Hr1G073440) and VRS4 (HORVU3Hr1G016690) were analyzed by Sanger sequencing (table S17).

Phenotypic analysis of mutants

The LSs from cul4 mutants, dep1 mutants, and SPL14 KO lines were extracted from main culm spikes approximately when their growth was arrested. The extracted LSs were imaged using a Zeiss Stemi-2000 stereo microscope (two independent experiments; n = 5 to 13 individuals). The morphological data, such as length, width, and area of LSs, were extracted from the image outlines using Axio vision SE64. Two-tailed Student’s t test was performed to find statistically significant differences between mutant and wild-type LS development. Phenotyping of spl14 mutant plants (six individuals) was carried out before harvest using the parameters plant height, fertile and nonfertile tillers, as well as the total number of tillers, the number and length of internodes and spikes, and the number of grains per plant. For significance test, one-way analysis of variance was performed with Holm-Sidak test by using SigmaStat 4 software.