High expression of the MADS-box gene VRT2 increases the number of rudimentary basal spikelets in wheat

Abstract Spikelets are the fundamental building blocks of Poaceae inflorescences, and their development and branching patterns determine the various inflorescence architectures and grain yield of grasses. In wheat (Triticum aestivum), the central spikelets produce the most and largest grains, while spikelet size gradually decreases acropetally and basipetally, giving rise to the characteristic lanceolate shape of wheat spikes. The acropetal gradient corresponds with the developmental age of spikelets; however, the basal spikelets are developed first, and the cause of their small size and rudimentary development is unclear. Here, we adapted G&T-seq, a low-input transcriptomics approach, to characterize gene expression profiles within spatial sections of individual spikes before and after the establishment of the lanceolate shape. We observed larger differences in gene expression profiles between the apical, central, and basal sections of a single spike than between any section belonging to consecutive developmental time points. We found that SHORT VEGETATIVE PHASE MADS-box transcription factors, including VEGETATIVE TO REPRODUCTIVE TRANSITION 2 (VRT-A2), are expressed highest in the basal section of the wheat spike and display the opposite expression gradient to flowering E-class SEPALLATA1 genes. Based on multi-year field trials and transgenic lines, we show that higher expression of VRT-A2 in the basal sections of the spike is associated with increased numbers of rudimentary basal spikelets. Our results, supported by computational modeling, suggest that the delayed transition of basal spikelets from vegetative to floral developmental programs results in the lanceolate shape of wheat spikes. This study highlights the value of spatially resolved transcriptomics to gain insights into developmental genetics pathways of grass inflorescences.


Introduction
The arrangement of flowers in individual plants of the same species is highly conserved and follows a systematic and rhythmic pattern. This systematic appearance of flowers is not surprising, as floral architectures are determined by the regular initiation of flower primordia on the flanks of the apical meristem and their rate of initiation and developmental fate are under strong genetic control (Prusinkiewicz et al., 2007). The unifying feature of floral architecture in grasses (Poaceae) is the formation of all flowers (termed florets) within spikelets (Kellogg et al., 2013). Spikelets are the fundamental building blocks of grass inflorescences and their development and branching patterns determine the various inflorescence architectures of grasses (e.g. spikes or panicles). Wheat (Triticum aestivum) forms a spike-shaped inflorescence, in which sessile spikelets are directly attached to the inflorescence axis (or rachis) in a distichous phyllotaxis (Koppolu and Schnurbusch, 2019). Upon floral transition, the vegetative meristem ceases to initiate leaf primordia and transitions into the inflorescence meristem (IM). During the Double Ridge (DR) stage of wheat spike development, the IM initiates a lower leaf ridge and an upper spikelet ridge (or primordia) during each iteration. Within the inflorescence the upper ridges differentiate into spikelet meristems, while the lower ridges are suppressed upon flowering (Bommert and Whipple, 2018). DR initiation will continue at the IM until the terminal spikelet stage, when IM forms a final spikelet (Koppolu and Schnurbusch, 2019). Spikelet initiation and development have been extensively studied in wheat and other monocot crops, such as rice (Oryza sativa), maize (Zea mays), and barley (Hordeum vulgare), as the number of spikelets per spike is a major determining factor for grain number and thus yield per spike.
Not all spikelets across the wheat spike, however, produce the same amount of grain. The central spikelets produce the most and largest grains, while spikelet size gradually decreases acropetally and basipetally. Within a single spike, the most apical and basal spikelets might produce no or only one grain while the central spikelets of the same spike set 3-5 grains. Bonnett (1966) documented that this distinct lanceolate shape of the wheat spike is first established during the Glume Primordia (GP) stage (just after the DR stage). This asynchronous development among the spikelets is maintained throughout the development of the spike. The gradual decrease in spikelet size from the central to apical section of the spike can be explained by the continuous development of new spikelet ridges from the apical IM: the most apical spikelets are the youngest and had the least time to develop. However, basal spikelets are initiated first and it is unclear why they remain smaller than their central counterparts. In the mature spike the most basal one or two spikelets are often only formed in a rudimentary manner, with small glumes present but all floral structures remaining immature.
Efforts to understand the genetics of wheat spikelet initiation and development have focused on members of the MADS-box transcription factor (TF) family, which play central roles in the flowering gene models (Zhao et al., 2006). Li et al. (2019) showed that MADS-box genes of the SQUAMOSA-clade, VERNALISATION 1 (VRN1), FRUITFULL 2 (FUL2), and FUL3, have overlapping functions in controlling the timing of the transitions from the vegetative to IM as well as the formation of the terminal spikelet. In vrn1ful2null mutants, the IM remained indeterminate causing the mutants to form more spikelets per spike. However, all lateral spikelets were replaced by leafy shoots in the vrn1ful2 double and vrn1ful2ful3 triple mutants (Li et al., 2019). These mutants had increased expression of genes belonging to the SHORT VEGETATIVE PHASE (SVP) family of MADSbox genes, including VEGETATIVE TO REPRODUCTIVE TRANSITION 2 (VRT2). Subsequent studies determined that overexpression of VRT2 led to reversion of basal spikelets to spikes and the downregulation of other MADS-box genes required for floral development, including members of the SEPALLATA 1 (SEP1) clade . Together, these studies exemplify the importance of the temporal sequence of flowering gene expression for the correct development of the wheat spike.
Attempts to unravel the genetic network controlling wheat spike development have focused on these temporal changes in expression patterns across consecutive developmental stages. For example, Li et al. (2018) and Feng et al. (2017) performed transcriptome profiling using pooled samples of multiple complete spikes from six (vegetative to floret differentiation) and four (DR to young floret) developmental stages, respectively. In a few cases, studies have examined the expression patterns of individual genes (via reverse transcription-quantitative PCR [RT-qPCR]) and found gene expression gradients along the spike. For example, Debernardi et al. (2017) demonstrated that APETALA 2 (AP2) is expressed higher in the apical section of wheat spikes than in central or basal sections. This AP2 expression gradient was associated with morphological changes along the same spike. This study alongside work in barley (Youssef et al., 2017) suggests that gene expression gradients within individual developmental stages could be important to further unravel the genetic control of spike development. However, despite its potential biological importance, spatial transcriptome profiles along the spike have yet to be investigated in wheat.
In this study, we aimed to characterize gene expression profiles along the spike during the establishment of the lanceolate shape of the wheat spike from DR to GP. We adapted Genome and Transcriptome sequencing (G&T-seq), a low-input sequencing approach to sequence the transcriptome of the sections. Recently, Giolai et al. (2019) adapted the protocol to identify expression differences across single leaves of Arabidopsis (Arabidopsis thaliana; GaST-seq), demonstrating that the G&T-seq method can be readily used for sequencing of hand harvested, small input plant material without the need of previous tissue dissociation or treatment. G&T-seq is thus comparable to methods using laser-micro dissection followed by sequencing to achieve spatially resolved transcriptome-wide sequencing data. In comparison, the available transcriptome sequencing methods at higher resolutions (such as single-cell RNA-seq or fluorescence-activated cell sorting [FACS]) are not spatially resolved as the complete tissue is dissolved into single cells for barcoding or selection prior to sequencing (Rich-Griffin et al., 2020).
We sequenced the apical, central, and basal sections of individual spikes before (DR) and after (GP) the establishment of the lanceolate shape. Gene expression profiles differed most strongly between spatial sections of the same spike, as opposed to temporal sections (any two sections from different time points). Members of the SVP gene family were expressed most highly in the basal sections with expression decreasing upwards from the base (acropetally), while members of the SEP1 gene family showed the opposite expression pattern, that is, most highly expressed in apical sections, with expression decreasing toward the base (basipetally). The increased number of rudimentary basal spikelets (RBSs) due to VRT-A2 misexpression supports the hypothesis that high expression levels of SVPs in the basal section delays spikelet establishment, leading to their rudimentary shape in the mature spike. This study highlights that spikelets within the same spike experience significantly different flowering signals due to their consecutive development and spatial position within the spike. Acknowledging these differences can help us gain a better understanding of the genetic flowering pathway of grass inflorescences.

Results
Low-input sequencing enables spatial analysis of the wheat spike transcriptome To investigate transcriptional differences between the apical, central, and basal section of developing wheat spikes, we adapted the low-input GG&T-seq method for RNA-seq of small plant tissue sections. G&T has been developed for single-cell RNA and DNA sequencing of mammalian systems (Macaulay et al., 2015) and was previously adapted for Arabidopsis thaliana (GaST-seq; Giolai et al., 2019). We collected four individual developing wheat (cv Paragon) spikes at both the DR and GP stage and hand-dissected them into apical, central, and basal sections ( Figure 1A).
On average, samples had 28,799,626 reads (coefficient of variation [CV] 43%), of which 90% (CV 8.5%) aligned to the genome post adaptor trimming ( Figure 1B and Table 1). Furthermore, the number of aligned reads and the number of expressed genes per library was largely homogeneous among the spatial sections and Waddington stages (Table 1). On average, 47,313 genes per library were expressed (410 read counts) and we found no difference (P 4 0.56, analysis of variance [ANOVA]) in the number of expressed genes across spatial (apical, central, and basal) or between temporal (DR, GP) conditions ( Figure 1C). We excluded three libraries with low average number of expressed genes (difference 45 times the standard deviation; Figure 1C; Supplemental Table S1) and two libraries because they were strong outliers in the principal component analysis (PCA; Supplemental Figure S1A). In total, 19 RNA-seq libraries (DR: 3 apical, 4 central, and 3 basal; GP: 2 apical, 3 central, and 4 basal) passed our selection criteria and were used in the subsequent analyses. We identified 91,646 genes being expressed across these 19 libraries.
Transcriptome-wide differences are largest between the apical and basal sections of the spike To investigate global differences among the 19 RNA-seq libraries, we performed a PCA (Figure 2A). The first two PCs explained 19% and 16% of the overall variance present in the libraries. We observed that the two PCs separated libraries by the spatial position (apical, central, and basal) rather than developmental stage (DR and GP). There was a clear separation between libraries originating from apical and basal spike sections, while libraries from central sections were dispersed between these two clusters ( Figure 2A). We investigated PC1-PC6 and found that none of these combinations clustered libraries by developmental stage (Supplemental Figure S1B). Given that we sequenced developing spike sections of single plants, as opposed to the more commonly employed pooling of multiple biological samples, we found as expected some degree of heterogeneity between samples from the same location and stage ( Figure 2A).
To investigate this variation further, we quantified changes in gene expression across biological replicates by calculating CVs for each gene (see "Materials and Methods"). The median CV for a gene across the biological replicates was 39% (Supplemental Figure S2) with a Q1-Q3 interquartile range between 24% and 62%. We also calculated the CV per gene for published datasets from Li et al. (2018) and Feng et al. (2017). Both studies sequenced developing wheat spikes at similar developmental stages, pooling many spikes per sample. Li et al. (2018) pooled between 100 and 200 spikes of winter wheat (KN9204) per sample, while Feng et al. (2017) reported pooling of 10-50 spikes (cv. Chinese Spring) per sample. In both studies, the median CV of a gene was lower (14% and 21%, respectively) than in our study. The larger CVs in our data could be explained by the biological variation that exists between individual plants, which may have been reduced by the pooling of many spikes in both Li et al. (2018) and Feng et al. (2017).
We first analyzed differentially expressed genes (DEGs) between the DR and GP stage and between apical, central and basal sections across the two Waddington developmental stages. The number of DEGs between DR and GP (215 genes) was smaller than the number of DEGs identified between the spatial positions, which ranged from 280 DEGs between central and apical sections to 914 DEGs between the apical and basal sections ( Figure 2B). Next, we compared the apical, basal, and central sections within each Waddington developmental stage. We identified more DEGs by comparing the spatial sections within either Waddington stage individually than in the combined analysis. The number of DEGs between apical and basal sections at each stage (DR: 2,438; GP: 1,315) was similar to the number of DEGs between central and basal sections (DR: 2,022; GP: 1,186). The number of DEGs between these sections at DR, however, was nearly double the number of DEGs at GP. In contrast, the number of DEGs between apical and central sections was similar at both stages (DR: 1,296; GP: 1,118), suggesting that the basal section of the spike is most different in the earlier developmental stage. Only 11% of the DEGs were shared between DR and GP in the apical to basal comparison, 7% between the central to basal DEGs, and 5% between apical to central DEGs. In total, we identified 5,353 unique genes as differentially expressed between any of the three sections at either Waddington stage (Supplemental Table S2). Overall, the number of DEGs was largest between the apical and basal sections, reflecting the strong spatial clustering observed in the PCA graph, but most genes that were differentially expressed across the spike did not maintain this gradient over the two developmental stages. In summary, despite the high biological variation in gene expression in our data compared to previous pooled whole-spike studies, we could detect transcriptome-wide differences between the spatial sections of developing wheat spikes.  Table S1. The SVP MADS-box TFs have opposing expression profiles to flowering E-class genes To further investigate the differences in expression across the spike and to identify genes with similar expression patterns, we performed hierarchical and k-means clustering. We restricted the clustering to the 5,353 genes identified as differentially expressed across the spike at either one or both Waddington stages ( Figure 3A). We identified seven nonredundant clusters, each containing between 8% and 21% of the 5,353 DEGs ( Figure 3A; Supplemental Figure 3A). Both hierarchical and k-means clustering produced highly similar results (Supplemental Figure S3B). We identified 1,894 genes (35% of DEGs) to be more highly expressed in the apical section, either across both time points (503 genes, Cluster 1), or only at DR (751 genes, Cluster 7) or GP (640 genes, Cluster 3). In the central section, 1,362 genes (25%) had higher relative expression at either DR (917 DEGs, Cluster 5) or GP (445 DEGs, Cluster 6). In the basal section, we observed 2,097 genes (39%) being more highly expressed. Cluster 4 contained the most DEGs (1,170) and was characterized by an upregulation of expression in the basal section at both Waddington stages, although this upregulation was higher at the DR stage. Another 927 genes were upregulated in the basal section, but only at the GP stage (Cluster 2). To further characterize the clusters, we independently tested for enrichment of TF families and gene ontology (GO) terms relating to biological processes (all GO terms and TF families in Supplemental Tables S3 and S4, respectively). Genes that were more highly expressed in the apical section were enriched for the GO terms "reproductive structure development" (GO:0048608) and "floral organ development" (GO:0048437; Cluster 1) as well as for HD-Zip_IV and SRS TF families (P 5 0.001; Cluster 1). In Cluster 3 (highly expressed in the apical section at GP; Figure 3B), we found no significant enrichment of GO terms relating to biological processes, but a significant enrichment of MADS_II TFs (P = 0.013). For clusters defined by an increased expression in the central sections (Clusters 5), we detected an enrichment for GO terms related to polyphosphate processes (GO:0006797/0006779; P 5 0.03), and a significant enrichment of the C2C2_CO-like TFs at DR (P 5 0.04). Genes with higher expression in the basal section of the spike (Cluster 4; Figure 3C) were enriched for a number of GO terms related to photosynthesis (e.g. GO:0015979; P 5 0.01) and "negative regulation of flower development" (GO:0009910; P 5 0.01), as well as for MADS_II TFs (P = 0.08). Cluster 2, which was characterized by higher expression in the basal section only at GP, was enriched for "Jasmonic acid response" (GO:0009753; P 5 0.001) and Tify (P 5 0.001) and C2C2_CO-like TFs (P 5 0.04).
We were interested in further characterizing the expression patterns of the MADS-box TFs as they were significantly enriched in two of the seven clusters and are important in floral transition and development (Becker and Theißen, 2003;Feng et al., 2017). We detected 14 differentially expressed MADS-box TFs in our study, 5 of which were more highly expressed in the apical section, 4 of these only at GP stage (Cluster 3), and 1 being consistently expressed across both Waddington stages (Cluster 1). In contrast, five MADS-box TFs were more highly expressed in the basal section at both Waddington stages (Cluster 4) and another two were more highly expressed in the basal section only at GP (Cluster 2). An additional two MADS-box genes were part of the remaining clusters (Supplemental Table S4).
In the apical/GP Cluster 3, we noticed that all MADS-box genes belonged to the T. aestivum SEP1 group ( Figure 3C). All three homoeologs of SEP1-4 (TraesCS7A02G122000, TraesCS7B02G020800, and TraesCS7D02G120500) and the Dgenome copy of SEP1-5 (TraesCS7D02G120600) were part of this cluster. The SEP genes were expressed at relatively low levels at DR (Supplemental Table S2), but were significantly upregulated at GP, with their transcript levels being highest in the apical section. The increased expression of SEP1-4 at GP was in agreement with their previously reported expression patterns in tetraploid wheat by Li et al. (2021).
In the contrasting Cluster 4 (upregulation in basal sections), we noticed the presence of multiple MADS-box genes belonging to the SVP family ( Figure 3B, right), which consists of three genes in wheat (SVP1, VRT2, and SVP3). Members of this family are important for the transition from vegetative to floral meristem identity in cereals (Trevaskis et al., 2007). All three homoeologs of VRT2 (TraesCS7A02G175200, TraesCS7 B02G080300, and TraesCS7D02G176700) and the B-genome copy of SVP1 (TraesCS6B02G343900) were present in Cluster 4. The cluster also contained FLOWERING LOCUS C (FLC-A1; TraesCS7A02G260900), although it was expressed higher in DR-central sections compared to the SVPs and had a linear expression gradient at GP. All SVPs had very similar expression patterns, being strongly expressed in basal sections only. Expression of SVPs was higher in all DR sections compared to the equivalent section in GP. Constitutive over-expression of SVP-family members in wheat and barley has been shown to delay or even reverse floral development (Trevaskis et al., 2007;Li et al., 2021). This led to hypothesis that the rudimentary development of basal spikelets was associated with an increase in VRT2 expression levels.

SVP expression is higher in basal and peduncle sections and increased across all sections in T. polonicum VRT-A2b isogenic lines
To validate the expression pattern of VRT2 in the individual spike analysis, we performed reverse transcription quantitative PCR (RT-qPCR) on independently collected, pooled  Supplemental Table S2. spike sections from cv Paragon, carrying the wild-type VRT-A2a allele ( Figure 4A (blue curves); Supplemental Table S5). We included a later time point, Terminal Spikelet (TS), which is about 10 d after GP to study how VRT2 expression changes in later stages. At TS stage, the central spikelets have developed multiple florets primordia. We also included a small part of the peduncle (stem) section just below the spike as an additional spatial section. We focused the expression analysis on the A-genome homoeolog, VRT-A2, as its role in spike, glume, and grain development of wheat was recently characterized (Adamski et al., 2021;Liu et al., 2021).
We identified a significant interaction effect between Waddington stage and spatial section (P 5 0.0001; Supplemental Table S6); we thus analyzed the three Waddington stages separately (Figure 4). At DR, we were limited to dissecting the spike into apical, basal, and peduncle sections, as the small size of the spike meristem did not allow precise dissection of the central section when using multiple (pooled) spikes. At DR, we found VRT-A2 marginally expressed, with significantly lower expression levels in the apical section compared to the basal (P = 0.003) and peduncle sections (P = 0.001). Although expression in the peduncle was higher than in the basal section at DR, this was not significant (P = 0.116). At GP, VRT-A2 expression was borderline detectable in the apical, central and basal sections, but expression was significantly higher in the peduncle with respect to the three spike tissues (P 5 0.001 for all three comparisons). Lastly, VRT-A2 expression at TS stage was just detectable and significantly different between all sections (P = 6.6E-06). Overall, expression decreased significantly from DR to GP/TS Waddington stages in the apical (P = 0.00015), basal (P = 0.0074), and peduncle (P = 0.012) sections consistent with the previously reported strong downregulation of VRT-A2 in the early wheat spike development (Adamski et al., 2021;Li et al., 2021;Liu et al., 2021). This is also consistent with the observed downregulation of VRT-A2 orthologs upon floral transition in barley (Trevaskis et al., 2007) and rice (Harrop et al., 2016). As observed in the low-input RNA-seq data, the RT-qPCR data confirmed the strong basipetal gradient in VRT-A2 expression across the spike at DR and revealed that its expression was even higher within the peduncle.
We hypothesized that the higher expression in the basal section of the wheat spike compared to the central and apical sections is associated with the rudimentary development of the basal spikelets. To test the effect of higher VRT-A2 expression on basal spikelet development, we analyzed the effect of the Triticum turgidum ssp. polonicum VRT-A2b allele on the expression gradient of VRT-A2 and spike morphology. Adamski et al. (2021) showed that VRT-A2 in T. polonicum, a tetraploid subspecies of wheat, carries a sequence re-arrangement in its first intron. This results in the higher expression of the T. polonicum VRT-A2b allele, with respect to the wild-type VRT-A2a allele, during early spike development. We performed RT-qPCR on a cv Paragon near isogenic line (NIL) carrying the VRT-A2b allele and compared VRT-A2 expression against the Paragon wild-type NIL Figure 4 Effect of VRT-A2b allele on VRT-A2 and SEP1-4 expression in the spike. Relative expression (2 ddCT ) of VRT-A2 (A) and SEP1-4 (B) in the different sections of the spike across three time points in NILs carrying either the wild-type VRT-A2a (blue) or the VRT-A2b allele from Triticum turgidum ssp. polonicum (orange). The data are shown as mean ± SE of gene expression compared with control gene Actin. N = 3 biological replicates (see Supplemental Table S5 for expression data and Supplemental Table S6 for statistical analysis of gene expression differences). described above (Figure 4; Supplemental Table S5). Consistent with the results of Adamski et al. (2021), we detected significantly higher expression of VRT-A2b compared to the wild-type allele across most of the tissue sections (see Supplemental Table S6 for individual comparisons), and a progressive decrease in VRT-A2b expression over time (P = 0.031). In contrast to the wild-type NILs, analysis of variance (ANOVA) did not identify a significant interaction effect between spatial section and Waddington stage in VRT-A2b NILs (P = 0.18). We thus examined the overall expression patterns and found that across all three developmental stages VRT-A2b expression differences were significant (P 5 0.0001). These results suggest that the basipetal expression gradient in the spike is maintained in the NILs with the T. polonicum VRT-A2b allele.
We also tested the effect of VRT-A2 expression levels on SEP1 expression in the tissue sections of VRT-A2 NILs. We confirmed that SEP1-4 expression is only marginally detectable at DR and differences in expression between the apical, basal, and peduncle section are hardly detectable with only marginally higher expression in the apical sections (P = 0.015) at this stage. At GP, SEP1-4 expression is significantly higher toward the tip of the spike (P 5 0.0001) consistent with the low-input RNA-seq data. Furthermore, SEP1-4 expression was significantly lower in VRT-A2b NILs compared to the wild-type allele across all spike sections (P = 0.008), confirming that higher VRT-A2 expression can negatively affect SEP1-4 expression. Similar trends were observed at TS, where expression was significantly lower in basal sections (P 5 0.0001) and the gradient across the spike was maintained in VRT-A2b NILs, but expression was overall lower (P 5 0.003). SEP1-4 is not expected to be expressed in vegetative tissue such as the peduncle; therefore, the lack of expression in this tissue across the three stages indicates that no floral tissue was accidently sampled as peduncle ( Figure 4B).

Misexpression of VRT-A2b in T. polonicum increases RBS numbers
To evaluate if the higher expression of VRT-A2 in basal spikelets affects their development, we examined the VRT-A2 NILs (BC 4 and BC 6 ) sown as winter crops in four environments. In each field trial, we evaluated the number of RBS, that is spikelets which are reduced in size and do not contain mature grains ( Figure 5A). The number of RBS was significantly increased in NILs carrying the VRT-A2b allele in all four environments (P 5 0.0001, except Morley, 2017 P 5 0.01; Figure 5B; Supplemental Table S7). The VRT-A2a NILs and the recurrent parent Paragon had on average 1.85 RBS, whereas VRT-A2b NILs produced on average 2.91 RBS. A similar difference in RBS between the NILs was observed in glasshouse conditions (VRT-A2b effect of + 1.6 RBS; Supplemental Table S8; Figure 5C). Furthermore, we also recorded the number of RBS in seven homozygous BC 6 recombinant lines used to fine-map VRT-A2b by Adamski et al. (2021). The RBS phenotype was mapped in complete linkage with the 50.3-kbp interval containing VRT-A2 ( Figure 5C; Supplemental Table S8). This genetic and phenotypic data suggests that the increase in RBS is a pleiotropic effect of the T. polonicum VRT-A2b allele and supports the hypothesis that misexpression of VRT-A2 negatively affects spikelet development in the base of the spike. In Paragon, the first (sometimes second) RBS fully develops the floral organs of florets one and two (e.g. lemma, palea, stamen, and ovary); however, these are severely reduced in size and delayed in development compared to the florets of central spikelets just before flowering ($Waddington stage 8-10; Supplemental Figure S4). At this stage, the further growth and development of these basal florets is stopped and in the mature spike only the glumes of RBS are visible. In NILs carrying the VRT-A2b allele, the development of the most basal spikelet is very similar to the wild-type. However, the second, third, and sometimes fourth spikelet also display similar signs of reduced development, leading to the larger number of RBSs ( Figure 5, B and C; Supplemental Figure S4).
To validate the phenotypic effect of VRT-A2, we analyzed transgenic wheat lines transformed with the complete genomic T. polonicum VRT-A2b sequence (including the native promoter and the intron 1 re-arrangement). Transgenic T 1 lines were classified based on the transgene copy number which was previously shown by Adamski et al. (2021) to be highly correlated with VRT-A2 expression levels in multiple tissues. We phenotyped lines with zero (n = 2 independent events; 5 plants each), low (1-5 transgene copies; n = 4 independent events; 5 plants each) and high (9-35 transgene copies, n = 2 independent events; 5 plants each) transgene copy number. We identified a significant and stepwise increase in the number of RBS with transgenic copy number, from 0.8 RBS (zero copy) to 1.6 RBS (low copy; P = 0.078 versus zero copy) and to 4.3 RBS (high copy; P 5 0.0001 versus zero copy) ( Figure 5D; Supplemental Table S9). The low copy number lines had an average increase of 0.8 RBS with respect to the zero copy number lines, equivalent to the average difference between the VRT-A2a and VRT-A2b NILs in the field (VRT2-A2b effect of + 1.1 RBS). The high copy number lines produced on average 4.3 RBS, which is higher than the VRT-A2b NILs and similar to the number of RBS observed in T. polonicum (3.75 ± 0.62 RBS; n = 16 spikes). The dosage-dependent effects observed in the transgenic lines provide further evidence that elevated expression of VRT-A2 leads to increased number of RBSs in polyploid wheat.

Discussion
High-resolution spatial transcriptomics in crops and barley (VRS2; Youssef et al., 2017). However, currently available transcriptome data (e.g. Li et al., 2018 andFeng et al., 2017) lack the spatial resolution within each individual developmental stage to answer this question. This focus on "between stage" comparisons (as opposed to within a single stage) is perhaps related to the technical challenges of dissecting and sectioning young meristems. Given the relatively small size of these spike meristems (0.2-mm length at Transition Stage; 3-mm length at Terminal Spikelet stage), RNA-seq methods require bulking of multiple individuals (usually between 30 and 50 different plants) to accumulate enough tissue for a single RNA-seq sample. If one sought to further section each meristem, this would require even further bulking. While laborious, this is achievable; however, under this scenario, the challenge is to properly stage approximately 100 plants to an equivalent developmental stage. Furthermore, it can be technically challenging to section these young spikes each time into the exact same apical, central, and basal sections. Consequently, the spatial resolution in gene expression within a wheat spike at individual developmental stages has remained largely uncharacterized to date.
To address this challenge, we adapted the G&T method for micro-scale spatial-transcriptomics workflow (Macaulay et al., 2015;Giolai et al., 2019), to conduct RNA-seq of the apical, central, and basal sections of individual, hand-dissected wheat spikes. This highly automated workflow requires low tissue input and allowed us to combine 24 Nextera libraries into a single Illumina NovaSeq lane. For 19 out of the 24 samples, the method worked successfully, determined by more than 20,000 expressed genes per library and the clustering among biological replicates. We found that the number of expressed genes per library was on average similar to the number of genes reported for bulked whole spike RNA-seq samples (Feng et al., 2017;Li et al., 2018). This is consistent with the fact that the handdissected sections are composed of a large mixture of different tissues (e.g. rachis, spikelet, and floret primordia) and cell types, which in the equivalent maize ears have distinct expression profiles (Xu et al., 2021). Compared to previous bulk RNA-seq studies in developing wheat spikes, the variation observed here (measured as CV) was high among biological replicates (Supplemental Figure S2). This variation is likely caused by both biological variation (e.g. inherent variation of individual plants) and technical variation (e.g. inaccuracies in sectioning and in the developmental staging of the plant/spike) as well as the number of replicates in our analysis. A minimum of six replicates has been proposed for bulked RNA-seq (Schurch et al., 2016). Our results suggest that the RNA-seq from these small sections would benefit from a higher number of biological replicates, which should be feasible considering the high-throughput method employed for RNA extraction and library preparation, the low tissue input requirement, and the possibility to pool multiple biological replicates per sequencing lane. Despite some limitations, we could identify over 5,000 DEGs between the spatial sections for subsequent functional analysis.
In addition to G&T-Seq, several other technologies have been proposed for obtaining high-resolution transcriptional profiles of plant tissues, for example, single-cell RNA-seq (McFaline-Figueroa et al., 2020;Rich-Griffin et al., 2020), fluorescence-activated cell sorting (FACS), and the isolation of nuclei tagged in specific cell types. These methods, however, are not spatially resolved as the complete tissue is dissolved into single cells for barcoding or selection (Rich-Griffin et al., 2020). Thus, these current methodologies do not allow, for example, to investigate whether the cell-type composition of spikelets differs across the inflorescence. This would only be possible if spikelets were "harvested" individually, for example through laser capture microdissection (LCM) before dissolving the tissue further into individual cells. Thiel et al. (2021) recently combined LCM followed by RNA-seq of the distinct lower/leaf ridge and upper/spikelet ridge of barley spikes. This allowed them to identify precise spatio-temporal expression patterns of many genes related to architecture and yield in barley spikes with unprecedented resolution. Looking ahead, increased resolution of Spatial Transcriptomics (currently 100 mm; Giacomello et al., 2017), which quantifies full transcriptomes while maintaining tissue integrity, offers the true prospect of direct localization and quantification of gene expression. Our results argue strongly for the need of these transcriptome-wide and spatially resolved approaches to advance our biological understanding of fundamental developmental processes in plants.

The composite nature of spikes
Early morphological studies of wheat spike development described that the stronger elongation of central spikelets during their initial establishment (GP stage) first causes the lanceolate shape of the wheat spike (Bonnett, 1966). The continuous formation of primordia at the tip of the spike means that at any given growth stage, spikelets in different developmental stages will be present across the spike (Bonnett, 1966). In this study, we detected more DEGs between the three spatial sections of the spike (apical, central, and basal) than between the two investigated developmental stages (DR and GP). We identified 215 DEGs between the two developmental stages, consistent with Li et al. (2018) who identified 206 DEGs between consecutive stages across a time course of six inflorescence development stages. Feng et al. (2017) identified 753 DEGs between the DR and Floret Primordia stage, which are further apart in development than the stages used in this study. They also detected fewer DEGs when comparing early stages than between more developed spikes. In contrast, we identified 1,315 and 2,438 unique genes to be differentially expressed between the apical and basal section at DR and GP, respectively. The higher number of DEGs between spatial sections could be due to the developmental gradients occurring in the three spatial sections, which are revealed by the spatial sampling. These differences would be blurred when comparing whole inflorescences between stages due to the mixture of tissue types and spikelets at different developmental stages. A possible improvement for future transcriptome studies could be the collection of only central sections of the developing spikes or complete spatial sampling as conducted here.
The composite nature of the inflorescence tissues has been acknowledged by studies in maize (ears and tassels), where new meristems are initiated in a stepwise manner. Leiboff and Hake (2019) quantified the meristematic tissue composition of maize and sorghum (Sorghum bicolor) tassels. For example, maize tassels in the second stage are mainly composed of spikelet pair meristems, but also contain some meristems in spikelet and inflorescence state. They concluded that the changes in these tissue compositions over time corresponded well with the independently staged transcriptional changes of the tassels. Eveland et al. (2014) showed that the range of developmental ages across the maize ear, if acknowledged, can be used as an advantage in RNA-seq studies. They sequenced the tip, middle, and basal sections of 10-mm long ears independently, aiming to analyze the expression patterns in specific developmental meristem types enriched in these sections (inflorescence, spikelet, and floral meristems, respectively). The dissection of the ear, therefore, allowed them to study gene expression specifically for each meristematic tissue type rather than for all meristem types in intact ears. In this study, we observed that apically expressed genes are enriched for GO terms related to "shoot system development" and "maintenance of floral organ identity." This is consistent with the hypothesis that the apical part of the inflorescence is younger and undergoing early phases of spikelet development initiation compared to the central inflorescence section.

Delayed transition of basal spikelets from vegetative to floral developmental programs
We detected transcriptional gradients across the spike, with the basal section deviating most strongly from the rest of the spike. We noticed that both SVP and CENTRORADIALIS (CEN) genes remained highly expressed in the basal section of the spike, whereas their expression was lower in the central and apical sections. In situ hybridization of these genes also showed that their expression is strongest in vegetative tissue and basal spikelets in early spike development . In contrast, SEP1-4 and SEP1-5 genes were expressed in the opposite gradient and showed the strongest expression in apical and central sections of the spike at GP stage. Recent studies allow us to interpret these gradients in the context of the early steps of vegetative to floral growth transition. In wheat (Adamski et al., 2021;Li et al., 2021;Liu et al., 2021), rice (Sentoku et al., 2005;Lee et al., 2008), and barley (Trevaskis et al., 2007), SVPs have been characterized to be associated with vegetative growth and are downregulated upon floral transition.
In wheat, the double SVP mutant vrt2svp1 leads to the formation of axillary inflorescences . Similarly, overexpression of CEN-D2 (TERMINAL FLOWER 1) in wheat extends the duration of the DR stage (Wang et al., 2017), whereas loss-of-function mutations in barley CEN suggest they repress floral development under short-day conditions (Bi et al., 2019). Double knockout mutants of the MADSbox SQUAMOSA genes vrn1ful2 highlighted that these two genes act as transcriptional repressors of SVP and CEN genes in early wheat spike development (Li et al., 2019). Furthermore, through a series of genetic and biochemical studies, Li et al. (2021) showed that the downregulation of SVP genes is necessary for the formation of flowering promoting MADS-box protein complexes including VRN1, FUL2, and SEP proteins. Hence the coordinated downregulation of SVPs, and possibly CEN genes, along with the upregulation of SEP genes is required for normal floral transition and spikelet development in wheat. Previous studies in rice have found similar expression patterns, as well as mutant effects, of SVPs and SEPs suggesting a conserved function in flowering transition across the two species (Ren et al., 2016;Wu et al., 2018).
Based on our results, the floral developmental program across the wheat spike appears to be most advanced in its apical and central sections, while being delayed in the basal sections. We hypothesize that this is due to elevated VRT2 expression at the base of the spike, which hinders the progression of the flowering program via SEP class flowering genes. Likewise, the higher expression levels of the wheat CEN2 and CEN5 homologs at the base are consistent with a delay in floral transition that could interfere with the development of the spikelet primordia. Therefore, although the basal spikelet primordia are initiated first chronologically, their developmental age in terms of the floral program is delayed with respect to the more recently formed central and apical spikelet primordia. This could explain in part why the spikelet primordia in the basal region of the spike elongate less and develop slower than central spikelets despite being initiated first (Bonnett, 1966). Likewise, the less advanced floral developmental program could also explain why the overexpression of SVPs in barley (HvBM10) leads to complete floral reversion in basal but not apical spikelets (Trevaskis et al., 2007).
We hypothesize that SVPs need to be downregulated upon floral transition to allow timely establishment and progression of the early spikelet primordia. Failure to do so would delay their development and result in their final rudimentary shape in the mature spike. In line with this hypothesis, we observed increased RBS in genotypes with prolonged and increased VRT2 expression in a dosagedependent manner. In our RT-qPCR data, we also observe increased expression of SVPs alongside reduced expression of SEPs in VRT-A2b lines at DR and GP stage. However, we cannot exclude the possibility that the increase in VRT2 expression in the VRT-A2b lines could also affect basal spikelet development at a later stage of spike developmental (e.g. from Terminal Spikelet stage to anthesis). We are currently using quantitative live imaging to compare cellular growth dynamics of spikelets at different stages of spike development between the VRT2 NILs.
The finding that the expected downregulation of SVPs and CENs does not follow the chronological age of the tissues suggests that other gradients across the spike might influence spikelet development. Debernardi et al. (2017) showed that in tetraploid wheat AP2-5 and miR172 have consistent and opposing expression gradients across the spike at three consecutive developmental stages. The persistent expression gradient of AP2-5 supports the idea that expression patterns across the spike, beyond the ones caused by age differences of spikelets, exist. Furthermore, they proposed a model illustrating that the phenotypic effect of mutants across the spike differs due to the existing gradient of expression of this gene (Debernardi et al., 2017). Other examples of mutants with different phenotypic effects across an inflorescence include vrn1ful2 (Li et al., 2019) in wheat, tassel sheath1 (tsh1; Whipple et al., 2010), and ramosa2 (Bortiri et al., 2006) in maize, SEP double mutant Osmads5Osmads34 in rice (Zhu et al., 2022 ), as well as many noded dwarf1 (mnd1; Walla et al., 2020), frizzy panicle (Poursarebani et al., 2015), and vrs2 in barley, which was also found to be consistently differentially expressed across the spike (Youssef et al., 2017). VRS2 has been shown to maintain a basal to apical expression pattern across three, post awn initiation developmental time points in barley (Youssef et al., 2017). The study of vrs2 mutants revealed that VRS2 is furthermore engaged with the basal-apical patterns of auxin, cytokinin, and gibberellin across the spike. While hormonal gradients across the spike in early development have not been studied in great detail in wheat, they have been shown to play crucial roles in floral induction and development in Arabidopsis (Reinhardt et al., 2000). Their patterns across the spike should be investigated in future studies addressing developmental differences across the spike.
A model for the regulation of leaf and spikelet ridge outgrowth in the base of the spike Recently, Meir et al. (2021) proposed that in shoot apical meristems of tomato (Solanum lycopersicum), similar to processes during embryonic development, transient programs are required to inhibit a preceding setup (i.e. vegetative growth), before a new developmental program (flowering) can be initiated. We propose that the altered gene expression and development of the basal spikelets could be a consequence of their initiation during the transient phase between vegetative and floral network shifts and thus being exposed to mixed signals of development. Upon floral transition, the lower (leaf) ridge is suppressed, while the growth of spikelet ridges from the previously suppressed axillary meristems is activated. Development of lower ridges subtending all branching events is suppressed in grass inflorescences upon flowering transition (Whipple et al., 2010). Li et al. (2019) noticed that this suppression was disrupted in the double vrn1ful2 and triple vrn1ful2ful3-null mutants, which fail to downregulate SVP genes. In these mutants, the upper spikelet meristems generate vegetative structures resembling tillers that are subtended by bracts or leaves originating from the lower leaf ridge.
We observed that genes that were highly expressed in the basal section of the inflorescence (Cluster 4) have previously been shown to be expressed specifically in the lower/bract ridge and before or at vegetative to floral transition. This is also supported by the GO term enrichment of photosynthesis related terms in Cluster 4. Our tissue sections do not allow us to distinguish lower and upper ridge tissues; however, the two ridges have been separately collected and sequenced via LCM in barley (Thiel et al., 2021). In this barley dataset, we found a higher expression of HvVRT2 (HORVU7Hr1G036130) and FLOWERING LOCUS C (HvFLC; HORVU7Hr1G054320) in the lower ridge compared to the upper ridge, whereas HvSVP1 (HORVU6Hr1G077300) was also marginally more highly expressed (Supplemental Figure S5). Furthermore, the barley MND1 gene (HORVU7Hr1G113480) has recently been shown to be expressed in leaf primordia and during the DR stage in the basal region of the spike in barley (Walla et al., 2020), while it is most highly expressed in the vegetative meristem and lower/leaf ridge in the LCM data (Thiel et al., 2021). We observed that in our data, the wheat MND1 putative orthologs (TraesCS7A02G506400, TraesCS7B02G413900, and TraesCS7D02G494500) were significantly more highly expressed in the basal section than the apical section at both DR and GP stage (Supplemental Table S2). The suppressed leaf ridge (or bract) has been proposed to act as a signaling center, regulating the fate of the upper spikelet meristem ridge (Whipple, 2017). Insufficient bract suppression during the formation of the basal spikelets might therefore negatively affect initiation and development of spikelets.
At DR, the widest point of the spike is indeed as expected the base and not the central section ( Figure 1A). The lower ridge is, however, much less developed in the central section and can be hardly seen in the apical ridges. Interestingly, mutants failing to repress the lower ridge growth, such as third outer glume1 (trd1; Houston et al., 2012), the barley ortholog of maize tsh1, develop large bracts from the lower ridge in basal spikelets, unlike apical spikelets, which do not develop bracts from their lower ridges regardless of the absence of TRD. This is reminiscent of the gradient in the strength of the phenotypic effects observed from the top to the base of the inflorescence in multiple Poaceae mutants (discussed above). We, therefore, hypothesize that the basal meristems develop into smaller spikelets and larger bract primordia due to a slow suppression of "vegetative growth signals" (e.g. SVPs) and a concomitant slow upregulation of "floral growth signals" (e.g. SEPs) upon floral transition. To investigate how a change from vegetative to floral signaling might affect the development of individual meristems, we modeled the genetic interaction of SVPs and SEPs, as proposed by Li et al. (2021), in the spatial context of a growing spike (Supplemental File S1). Under the assumptions that SVP suppresses SEP expression, SVP expression is downregulated upon flowering, and that SEP promotes spikelet outgrowth, the model could recapitulate (1) the observed opposing gradients in expression of SVPs and SEPs along the spike and (2) the formation of a lanceolate shaped wheat spike with reduced spikelet elongation and stronger bract growth in the most basal spikelets. Thus, whilst this hypothesis will require further investigation and testing, modeling supports its plausibility.
For field experiments, we used between two to four sibling BC 4 /BC 6 NILs differing for the VRT-A2b allele.
Four individual spikes per developmental stage (DR and GP) were dissected into apical, central, and basal sections (1:1:1 ratio) using a stereo microscope (Leica MZ16). Sections were immediately placed into 96-well plates (on ice) containing 10 lL of RLT plus (Qiagen, Hilden, Germany). All instruments and surfaces were cleaned with 80% v/v ethanol, RNAse-free water and lastly RNAse-out solution after each sample to reduce cross-contamination and RNA degradation. Samples were stored at -80 C until cDNA preparation, using the G&T-seq method as previously described (Macaulay et al., 2015). cDNA was normalized to 0.2 ng/lL before Nextera (Illumina, San Diego, CA, USA) library preparation using a Mosquito HV liquid handler (STP, Royston, UK) in a total reaction volume of 4 lL as described in Mora-Castilla et al. (2016). Libraries were pooled by volume and sequenced on a single lane of a NovaSeq 6000 (NVS200S2 flow cell, 100-bp paired-end reads).
HT-read count normalization and differential expression analyses were performed using the DESeq2 version 1.28.1 R packages (https://bioconductor.org/packages/release/bioc/ html/DESeq2.html; Love et al., 2014;RStudio 1.2.5001). Genes with an average expression below 10 HT-count, and which were not expressed (i.e. 10 HT counts) in at least three libraries, were removed from the analysis. Correlation between expressed genes and Waddington stage and/or section was tested by ANOVA. Raw read data from Li et al. (2018) and Feng et al. (2017) were pseudo-aligned using Kallisto Sleuth pipeline (https://scilifelab.github.io/courses/ rnaseq/labs/kallisto) and the CV was calculated for each gene (by condition) using R (RStudio 1.2.5001) ddply (plyr 1.8.6). DEGs between the two Waddington stages (DR and GP) were calculated with the design "$plant + section," while DEGs between the three sections (apical, central, and basal) were determined using the design "$section + waddington:plant + waddington." DEGs among the sections within each Waddington stage were determined with the design "$plant + section." For each gene, an adjusted Pvalue was computed by DESeq2 (using the using Benjamini and Hochberg method; Benjamini and Hochberg, 1995), and those with an adjusted P-value of 0.05 were considered differentially expressed. DeSeq2 also computed Log2FoldChanges as well as the associated uncertainty (Log2FoldChanges (lfcSE), see Love et al. (2014) for further detail). The "contrast" function was used to determine pairwise comparison P-values. The full set of expression data and comparisons is presented in Supplemental Dataset S1. Enrichment of GO terms was performed using the online tool "PLAZA" (https://bioinformatics.psb.ugent.be/plaza; Van Bel et al., 2017) using the recommended settings, and all enriched GO terms of Biological function and Cellular Compartment were retained. In brief, PLAZA determines the overrepresentation of a certain GO term in a gene set compared to the genome-wide background frequency (= all expressed genes in this experiment; submitted manually).
The significance of over-or underrepresentation is determined using the hypergeometric distribution and the Bonferroni method is applied to correct for multiple testing. Note that enrichment folds are reported in log2 fold scale. Enrichment of TF families (Genes that were annotated as TFs were obtained from https://opendata.earlham.ac.uk/ wheat/under_license/toronto/Ramirez-Gonzalez_etal_2018-06025-Transcriptome-Landscape/data/data_tables/; Ramirez-Gonzalez et al., 2018) and MADS-box TFs (based on Schilling et al., 2020) was performed in R using the phyper() function from stats package version 4.0.1 to test for Hypergeometric Distribution. All DEGs were scaled and centered using R-base function "scale." All cluster analysis was performed on scaled data using R (stats) functions kmeans and hclust, followed by visualization through pheatmap version 1.0.12 (https://cran.rstudio.com/web/packages/pheat map/index.html). Correlation to centroid cluster shape of each gene expression pattern was calculated using the "cor" function from R stats.

RT-qPCR analysis
P1 NILs were grown in controlled growth chambers in 24cell seed trays under the same conditions as used in the low input RNA-seq experiment (see above). For each biological replicate, we pooled 30 inflorescences for DR stage, 15 for GP stage, and 9 for Terminal Spikelet stage (n = 4 biological replicates per stage). Inflorescences from NILs were dissected using a stereo microscope (Leica MZ16). Inflorescences were dissected into apical, central, basal, and peduncle sections (1:1:1:1 ratio). At DR stage, inflorescences were only dissected into apical, basal, and peduncle section as the inflorescences were too small to be accurately dissected into four sections for all 30 plants per biological replicate. Each section was immediately placed into 1.5-mL tubes on dry ice and tubes were snap frozen in liquid nitrogen as soon as all plants for the sample were collected. Samples were stored at -80 C until needed. Inflorescences were collected within 2-3 h, 9 h after the lights came on in the growth chamber. Tissue was homogenized in a TissueLyser II (Cat No.: 85300; QIAGEN, Hilden, Germany) using 3-mm steel beads (Cat No.: 69997; Qiagen); tubes were shaken for 20-s at 28 Hz with dry ice.
All RNA extractions were performed using the RNeasy Plant Mini Kit (Cat No.: 74904; Qiagen) with RLT buffer according to the manufacturer's protocol followed by RNA ethanol precipitation (https://projects.iq.harvard.edu/files/hla lab/files/ethanol-precipitation-of-rna_hla.pdf). DNA digestion was performed using the RQ1 RNase-free DNase set (Cat No.: M6101; Promega, Madison, WI, USA) according to the manufacturer's protocol. RNA was reverse transcribed using M-MLV reverse transcriptase (Cat No.: 28025013; Thermofisher, Waltham, MA, USA) according to the manufacturer's protocol. For the RT-qPCR reactions, LightCycler 480 SYBR Green I Master Mix (Roche Applied Science, UK) was used according to the manufacturer's protocol. The reactions were run in a LightCycler 480 instrument (Roche Applied Science, UK) under the following conditions: 10 min at 95 C; 40 cycles of 10 s at 95 C, 15 s at 62 C, 30 s at 72 C; dissociation curve from 60 C to 95 C to confirm primer specificity. All reactions were performed with three technical replicates per sample and using Actin as the reference gene (Uauy et al., 2006). Relative gene expression was calculated using the 2 -DDCt method (Livak and Schmittgen, 2001) with a common calibrator so that values are comparable across genes, tissues, and developmental stages. All primers used in RT-PqCR came from Adamski et al. (2021) can be found in Supplemental Table S10.
All RT-qPCR data were normalized using a log2 transformation. A three-way ANOVA including Waddington stage, section, and genotype yielded significant two-way interactions. The differences between sections of the genotypes were therefore further analyzed individually for each Waddington stage and genotype. For each of the two genotypes we individually performed Tukey's multiple comparison tests to determine differences between the sections within each developmental stage by Tukey's multiple comparison test. Differences between the two genotypes were also analyzed individually for each Waddington stage. Furthermore, the differences between the genotypes were investigated individually for each section within the Waddington stage if the interaction term was significant (in GP and TS). For all analysis see Supplemental Table S6.  (2019), we analyzed two BC 6 and one BC 4 line per VRT-A2 allele and in Morley and Church Farm (CF) 2020, we analyzed two BC 6 and two BC 4 lines for each VRT-A2 allele. All experiments were drilled as yield-scale plots (6 Â 1.2 m) and sown by grain number for comparable plant densities aiming for 275 seeds m -2 . The trials were arranged in a randomized complete block design with five replicates per sibling line per location. A 10-ear grab sample was collected from each plot preharvest for the assessment of RBS numbers and other phenotypes (recorded in Adamski et al., 2021). RBS was defined as spikelets carrying no grain at maturity and counted for each spike individually. To determine the differences between the P1 POL and P1 WT NILs, we performed ANOVA on the multiple field trials phenotypic data. For the analysis of individual trials, we used a two-way ANOVA including Genotype + Block performed in R ("car" package version 3.0-10; RStudio 1.2.5001).

Glasshouse phenotyping
We evaluated the BC 4 NILs and BC 6 F 3 recombinant lines, as well as T. polonicum accession T1100002, under standard glasshouse conditions. About 18-20 plants per genotype were grown in 1 L pots containing John Innes Cereal Mix under long-day conditions (16-h light, 8-h dark). The genotypes of all plants were confirmed using KASP marker SP1Pol (Adamski et al., 2021). We counted the number of RBSs for all tillers of all biological replicates at maturity. To evaluate the differences in RBS between genotypes, we performed a two-way ANOVA analysis and post-hoc multi-pairwise comparisons Sidak test ("car" package version 3.0-10; RStudio 1.2.5001).

Phenotyping of transgenic lines
T 1 lines from Adamski et al. (2021) differing for the copy number of the VRT-A2b transgenic construct (zero = 0 copies; low = 1-5 copies; high = 9-35 copies) were grown in 1 L pots with John Innes Cereal Mix under 16-h light at 20 C and 8-h dark at 15 C in controlled environment growth chambers. We measured RBS number for the main tiller of all plants at maturity. To determine differences in RBS between the three transgenic classes, we performed ANOVA (two-way ANOVA; "car" package version 3.0-10). We performed Dunnett tests to compare the low and high copy lines against the zero copy number controls (RStudio 1.2.5001).

Modeling
The computational model of wheat spike shape formation was developed using the multi-agent programming language and modeling environment, Netlogo (Wilensky, 1999). Gene interactions were modeled as previously described . The model can be accessed via the interactive webversion of the model (Supplemental File S1).
In brief, both spikelets and leaves are initiated with rates that depend on the levels of SEP. Leaf initiation rates are suppressed by SEP, whereas spikelet initiation requires SEP. The maximum initiation rates are the same for both spikelets and leaves but different before (r vegetative ) and after (r flowering ) flowering. Once initiated, the leaves and spikelets grow at a rate defined by the parameters r leaf and r spikelet , respectively. Leaf growth does not depend on SVP or SEP levels, whereas spikelets only increase in size every iteration if their SEP level is above a given threshold (SEP growth_threshold ). Expression of both SVP and SEP only occurs at meristem initiation. After this, the levels of SVP and SEP cannot increase, although SEP is degraded. SVP is not degraded, solely because at this point, nothing is dependent on SVP levels, whilst spikelet growth depends upon SEP levels.
SVP expression rates start to decrease, once flowering is triggered, according to: where r SVP t ð Þ is the rate of SVP expression at that time step, and f red is a rate reduction factor. SEP expression depends upon the levels of SVP in the meristem in which the initiation points are located, depending on a Hill function (Alon, 2007), where r SEP;max is the maximum rate of SEP expression, K D is the binding constant, and n is the Hill coefficient. The resulting curves for SVP and SEP expression are shown in Supplemental Figure S6. SVP levels are initiated with the current value of r SVP . SEP levels are initiated using r SEP , and reduce by degradation rate, d SEP , following dSEP dt ¼ SEP Á d SEP :

Data availability
The raw RNA-seq read libraries used in this study are available from NCBI BioProject PRJNA749586.

Accession numbers
Accession numbers for the sequence data are listed throughout the study and in Supplemental Table S2, and are based on gene nomenclature from the IWGSC RefSeq version 1.0 annotation.

Supplemental data
The following materials are available in the online version of this article. Supplemental Figure S1. PCA of RNA-seq libraries.
Supplemental Figure S2. Comparison of CV of gene expression in wheat RNA-seq data sets.
Supplemental Figure S3. Expression patterns of all DEG clusters (n = 7) as identified by k-means clustering.
Supplemental Figure S4. Dissected Floret 1 and 2 of basal and central spikelets of BC 6 NILs before anthesis.
Supplemental Figure S5. Expression of barley genes in Thiel et al. (2021), which are putative orthologs to wheat genes highly expressed in basal spike sections.
Supplemental Figure S6. Example of Netlogo simulation outcome with default parameters.
Supplemental Table S1. Quality control measurements of all 24 RNA-seq libraries.
Supplemental Table S2. Summary of normalized gene expression, statistical analyses, and k-means clustering for all DEGs (Padj 0.05).
Supplemental Table S3. Enrichment of GO terms in the seven identified clusters of DEGs.
Supplemental Table S4. Enrichment of TF families and MADS-box TF genes in the seven identified clusters of DEGs.
Supplemental Table S10. List of primers used in RT-qPCR in this study.
Supplemental Dataset S1. Summary of normalized gene expression, statistical analyses, and k-means clustering for all expressed genes.
Supplemental File S1. Interactive web-version of the wheat spike model.