HvCslF6 is expressed across a range of tissues including grain and young/developing tissues such as the coleoptile, root and leaves (Burton et al., 2008). To study the transcriptional regulation of this gene in planta, transgenic barley lines carrying endoplasmic reticulum-targeted GFP (mGFP-ER) driven by a 3,000 bp 5′ sequence upstream of HvCslF6 (pHvCslF6) were generated and screened in grain and vegetative tissues. Transverse sections of transgenic mature grain showed strong mGFP-ER expression in the starchy endosperm (Figure 1A) compared to the negative control (empty vector) (Figure 1B). Similarly, increased mGFP-ER expression around vascular bundles was detected in transverse coleoptile sections of transgenic lines (Figure 1C) compared to the negative control (Figure 1D). Therefore, we confirmed that the 3,000 bp 5′ sequence upstream of HvCslF6 is sufficient to drive expression of this GFP gene in both grain and vegetative tissues.

As the 3,000 bp 5′ promotor region was sufficient to drive the temporal and spatial expression of HvCslF6 in developing and young tissues, we explored whether known transcription factor binding motifs were either present or enriched within this sequence. First, three TFBS prediction software packages (TRANSFAC^® v2014, JASPAR v2020 and PlantPAN v3.0) were compared to identify which produced the most detailed motif predictions with the least number of redundant results. TFBS prediction outputs were compared across the three programs, revealing differences in frequency (number of times a particular motif was predicted within the 3,000 bp pHvCslF6) and annotation (profile summary of each motif comprising TF class and origin species) of predicted TFBSs (Supplementary Figure 1 and Supplementary Table 1). It should be noted that a high rate of false positives in predicted TFBSs has been reported in most of these prediction tools (Shahmuradov et al., 2005). Therefore, choosing a comprehensive and experimentally validated database of TFBSs was key to obtain reliable in silico predictions. TRANSFAC^® was discarded due to the relatively low number of predicted motifs that were also identified using JASPAR and PlantPAN. A comparative motif analysis between PlantPAN and JASPAR identified seven classes of TFBSs with similar relative frequencies (allowing a ± 10% variation in motif prediction results across both programs). These were: MYB/SANT (12.5% in JASPAR and 14.0% in PlantPAN, respectively), DNA-binding with one zinc finger DOF (13.3% and 8.9%), APETALA2/Ethylene response factors AP2/ERF (9.3% and 8.4%), basic leucine zipper bZIP (4.7% and 6.5%), WRKY (5.4% and 3.2%), auxin response factors ARF (B3 domain; 3.1% and 5.2%), and NAC (4.6% and 3.6%). Notably, 52% of PlantPAN predicted motifs could not be classified because no information/description was available regarding their hit sequence.

Based on these results, JASPAR was chosen for subsequent in silico analyses of the proximal promoter region of HvCslF6 (1,000 bp upstream of the start codon) and 355 non-redundant motifs were identified. JASPAR TFBS scores are based on position frequency matrices (PFMs), defined as occurrences of each nucleotide at each position in a set of observed TF-DNA interactions (Fornes et al., 2020). We selected JASPAR for further analyses due to the higher number of non-redundant predicted TFBS and motif-associated information compared to TRANSFAC and PlantPAN. In addition, JASPAR is a regularly maintained open-access, manually curated and experimentally defined database for plant-specific TFBSs (Castro-Mondragon et al., 2022).

These 355 predicted motifs were filtered in accordance with an arbitrary prediction score threshold of ≥10 (maximum score was 15.9) to minimize the detection of false positives/redundant motifs in identical promoter locations (unfiltered TFBS predictions are shown in Supplementary Figure 2). By applying this filter, 107 unique predicted motifs with prediction scores between 0.7 and 15.9 were identified and grouped into eight classes of TFBSs including: MYB/SANT, AP2/ERF, basic helix-loop-helix (bHLH), ARF, DOF, Homeodomain-leucine zipper (HD-ZIP), Zinc Finger (ZF; SBP and SWIM type), NAC/NAM and a mixed class of TFBSs (low-represented putative TF classes with high prediction scores) (Figure 2A). Predicted motifs grouped into AP2/ERF, MYB/SANT and bHLH TF classes were the most abundant, hence over-represented within the HvCslF6 proximal promoter region (−1,000 bp) when compared to the total of 107 non-redundant motifs identified. JASPAR prediction scores for the three most abundant TF classes were similar to each other (11.6, MYB/SANT; 11.4, bHLH and 11.7, AP2/ERF; Shown on top of each column in orange, Figure 2A).

To understand if there was conservation of binding site motifs across grass species, we carried out in-silico analysis of the putative CslF6 promoter sequences in wheat (Triticum aestivum), Brachypodium distachyon and rice (Oryza sativa), and compared the TFBS to those identified in barley (filtered by a ≥ 10 prediction accuracy score cut-off, as previously applied to HvCslF6 full length promoter analysis). This revealed 59 common motifs across all species (Figure 2B and Supplementary Table 2) corresponding to: MYBs/SANTs, ARFs, AP2/ERFs, AT-Hook DNA-binding factors, bHLHs, DREBs, and HD-ZIPs TF classes. The second largest group of common motifs (26) was shared between barley, wheat and Brachypodium CslF6 upstream sequences and corresponded to putative binding sites for MYB-related TFs, (predominantly MYBs/SANTs: MYB55/61, MYB59, MYB111 and MYB113; GARP/G2-LIKE: PHL11, PHL12; MYB-related: UIF1), and other less abundant bHLH and DOF-related motifs. The rice CslF6 upstream sequence contained the highest number of unique TFBS predictions (84) compared to other cereals (13, barley; 13, wheat; and 14, Brachypodium). The frequency, location, and description of predicted TFBS across all species is shown in Supplementary Table 2.

Next, the analysis of TFBSs was expanded to include three more members of the Cellulose synthase gene superfamily, namely HvCesA1, HvCesA2 and HvCesA6, which are involved in cellulose biosynthesis in the primary cell wall (Burn et al., 2002; Burton et al., 2006) and known to be co-expressed with HvCslF6 (Burton et al., 2004; Burton and Fincher, 2009). Common predicted motifs across HvCesA1, HvCesA2, HvCesA6, and HvCslF6 corresponded to MYB-related, AP2/ERF, HD-ZIP, DOF, ZF-HD, Homeo Domain (HD), and MADS box factors. The largest set of common motifs was found across HvCesA2 and HvCesA6 putative promoters (77), followed by HvCesA2 HvCesA6, and HvCslF6 (60) (Figure 2C). As a comparison three HvCesA genes (HvCesA4, HvCesA7, and HvCesA8) which are not co-expressed with HvCslF6 were also analyzed for TFBSs. This analysis identified one common predicted TFBS, a SQUAMOSA promoter-binding protein-like (SPL) motif across all putative promoters (SPL TFs are mainly involved in plant growth and development; Tripathi et al., 2018) and abundant gene-specific, not shared, motif subsets. The largest set of TFBSs was found in the HvCslF6 promoter (57) containing MYBs, AP2/ERFs HD-ZIP among other TF classes. (Supplementary Figure 3 and Supplementary Table 3).

To identify the HvCslF6 promoter region necessary for transcriptional activation, seven progressively deleted promoter constructs were generated from the 3,000 bp region upstream of the HvCslF6 start codon (−3,000, −1,846, −1,357, −858, −607, −382, and −199 bp) and tested in a barley leaf-derived protoplast system (previously optimized using GFP/YFP-expressing constructs; Supplementary Figure 4). The region of the promoter that produced the highest relative luciferase assay activity was a −382 bp region upstream of the HvCslF6 start codon (Figure 3A). Additionally, deletion constructs containing a −607 bp or −199 bp HvCslF6 promoter region showed a significant decrease in promoter activity based on the relative luciferase assay compared to the −382 bp region. Therefore, we had delimited a 183 bp sequence within the proximal HvCslF6 promoter between the two smallest deletion constructs (from −382 to −199 bp) that showed a relative increase in luciferase activity and could be responsible for upregulation of HvCslF6 expression. HvCslF6 promoter activity was comparatively low from −607 to −3,000 bp upstream of the HvCslF6 CDS, indicating that repressors of HvCslF6 expression could be binding to this −382 to −607 bp region. Furthermore, no significant differences in promoter activity were observed across constructs carrying the longest HvCslF6 promoter fragments (−3,000, −1,846, and −1,357 bp; Supplementary Figure 5).

We retrieved HvCslF6 proximal promoter regions (1,000 bp upstream from start codon) from cv. Morex and cv. Golden Promise and carried out a pairwise sequence alignment. Both promoters were identical based on the barley pangenome sequences available at GrainGenes¹ therefore subsequent in silico prediction of TF binding sites could also be extended to cv. Golden Promise.

To further delimit putative TF binding sites within the region that showed increased luciferase activity, three additional deletion constructs were generated. Constructs containing −282 bp, −331 bp, and −357 bp of the HvCslF6 promoter, were tested alongside three previously tested deletion constructs (−199 bp, −382 bp, and −607 bp) to normalize dual luciferase measurements and compare promoter activity across deletion constructs (Figure 3B). A significant increase in promoter activity (p-Value 0.029, Tukey′s test) was detected in the −357 bp deletion construct. In-silico prediction of this promoter fragment indicates that it includes a putative TATA box element (5′–TATAAA–3′). This likely explains the increase in HvCslF6 promoter activity with the inclusion of this region of the promoter, and the low/zero levels of promoter activity in the first three constructs which based on our in-silico prediction lack a TATA box. Highest HvCslF6 promoter activity was narrowed down to a 51 bp region located −331 bp to −382 bp upstream of the start codon, although no significant differences (p-Value 0.321, Tukey’s test) were determined between the −357 bp and −382 bp deletion constructs. Moreover, PlantPAN and JASPAR in-silico motif analysis of the additional 25 bp included in the −382 bp fragment compared to the −357 bp fragment predicted a single MYB-related motif, AC-I, enriched in AC bases. Results also indicated low HvCslF6 promoter activity from −199 bp to −331 bp based on three deletion constructs analyzed (−199 bp, −282 bp, and −331 bp). This region contained two (three in the case of −331 bp construct) predicted calmodulin-binding NAC elements and a MYB-related motif, described as a secondary wall MYB-responsive element (SMRE) based on JASPAR analysis (Figure 3B).

Based on previous studies which suggest a key role for NAC and MYB TF families in SCW polysaccharide biosynthesis (Taylor-Teeples et al., 2015) and the involvement of MYB TFs in grass cell wall synthesis including hemicelluloses (Zhao et al., 2019), we mapped two putative MYB-related cis-elements (among other predicted motifs from our in-silico analysis) onto HvCslF6 promoter regions that conveyed increased promoter activity (Figure 3B and Table 1). A six-nucleotide, AC-I element (5′–ACCTAC–3′) was predicted from −364 to −358 bp upstream from the HvCslF6 start codon. Prouse and Campbell (2013), demonstrated that the R2R3-MYB AtMYB61 can bind to this AC-I element, therefore we were interested in the barley ortholog of this gene as a potential regulator of HvCslF6 expression. A putative binding site in construct −331 (5′–CGTTGGT–3′/3′–ACCAACG–5′) corresponds to another MYB consensus motif, also described as a secondary wall MYB-responsive element (SMRE) or a MYB46-responsive cis-regulatory element, M46RE (Zhong and Ye, 2012). This cis-element was detected in reverse orientation (3′–5′) and in Arabidopsis the same motif can be bound by AtMYB46 and AtMYB83 (Kim et al., 2012; Ko et al., 2014). Additionally, a secondary wall MYB-responsive element (SMRE) motif was predicted from −299 to −292 bp, downstream of the TATA box in the HvCslF6 promoter (see construct −331, Figure 3B). However, this had limited impact on promoter activity, likely due to the absence of the TATA box.

To assess protein sequence similarities in MYB46, MYB83, and MYB61 orthologs between barley (Table 2), Arabidopsis (Romano et al., 2012) and rice (Zhong et al., 2011), we constructed an unrooted phylogenetic tree of R2R3-type MYB transcription factors which includes those that could putatively bind to the −382 bp HvCslF6 promoter region (Supplementary Figure 6). The candidate MYB transcription factors that correspond to the predicted binding sites belong to the plant specific R2R3-MYB subfamily, which is characterized by the presence of two highly conserved MYB DNA-binding domains (Pfam PF00249; Bateman et al., 2017) and constitutes the largest MYB subgroup. HvMYB61 was grouped together with their orthologs from Arabidopsis (AtMYB61; Prouse and Campbell, 2013) and rice (OsMYB61a and OsMYB61b; Zhao et al., 2019). The best fitting barley ortholog for AtMYB46 and AtMYB83 was a single gene named HvMYB46/83, following previous studies in rice (OsMYB46/83; Rao and Dixon, 2018; Zhao et al., 2019) and grouped with its orthologs (Supplementary Figure 6).

Using gene expression data from 808 individual samples (Supplementary Table 4; Milne et al., 2021), we constructed a co-expression network to identify genes closely connected to HvCslF6. The samples covered a wide range of biotic and abiotic stresses, tissue types and cultivars, resulting in a robust dataset. Low expressed genes were removed (see section “Materials and methods”) and a weighted gene correlation network analysis (WGCNA) performed (Langfelder and Horvath, 2008). We identified 33 modules, with a color assigned to each module, plus one module (gray) which contains the unallocated genes. HvCslF6, along with the primary cell wall cellulose synthases (HvCesA1, HvCesA2, HvCesA6) were assigned to the yellow module (Supplementary Figure 7). HvMYB61 (one of the candidate HvCslF6 regulators) and HvTHX1, the barley ortholog of BdTHX1, which binds to the intronic region of BdCslF6 and acts as a positive regulator of expression of this gene (Fan et al., 2018), were all members of the yellow module. Visualization of the yellow network in Gephi showed a gene network comprising 992 genes and 54,090 gene connections (Supplementary Figure 8). To better investigate the connections, an intermediate stringency threshold of 0.188 (the correlation ranged from 0.15 to 0.335) was applied to define a subnetwork of the top 50 genes connected to HvCslF6.

This subnetwork based around HvCslF6 still included HvMYB61 (BART1_0-p01380 and the three cellulose synthases; HvCesA1, HvCesA2, HvCesA6, BART1_0-p60277, BART1_0-p40943, BART1_0-p44934) (Figure 4). All are either known to be involved in primary cell wall synthesis in barley, or in the case of HvMYB61 the ortholog of OsMYB61 has been shown to influence to cell wall synthesis (Zhao et al., 2019). Additionally, this module contained other proteins that have been linked to cell wall processes (Supplementary Table 5). COBRA-like protein 3 (BART1_0-p39134) located on chromosome 5H was shown to be co-expressed with HvCesA1, HvCesA2, and HvCesA6 (Houston et al., 2015). The other candidate regulator, HvMYB46/83 did not pass the filtering threshold when removing genes with low levels of expression or were assigned to the black instead of yellow module. Hence, this candidate regulator is unlikely to be a global regulator of HvCslF6.

Then we carried out a correlation analysis using transcript expression data from the eoRNA database (Milne et al., 2021), which contains 22 publicly available datasets and >800 RNA-seq datasets from a wide range of tissues, development stages and experimental treatments. We focused on two subsets of genes. The first included our candidate TF gene HvMYB61, a selection of genes from the network analysis and others known to be co-expressed with HvCslF6 (i.e., HvCesA1, HvCesA2, and HvCesA6). The second subset was comprised of secondary cell wall genes not co-expressed with HvCslF6 (i.e., HvCesA4, HvCesA7, HvCesA8), and included other members of the HvCslF/H/J gene family, and HvTHX1, the barley ortholog of BdTHX1, a positive regulator of BdCslF6 (Figure 5A and Supplementary Figure 9). We hypothesized that if HvMYB61 shows comparable or greater correlation with HvCslF6 compared to HvCesA1, HvCesA2, HvCesA6, this would provide further confidence in a relationship between these candidates. However, if strong positive correlation is observed between the expression of HvMYB61, our candidate gene, and the group considered to be secondary cell wall HvCesA genes (HvCesA4, HvCesA7, and HvCesA8), which show no correlation with HvCslF6, it might suggest that our candidate gene is unlikely to be involved in the broader regulation HvCslF6 expression.

For the transcripts of genes considered to be primary cell wall related CesA genes and HvCslF6, we observed Pearson correlation coefficients of between r = 0.9 to −0.1, while the secondary cell wall CesA genes had r = 0.4 to −0.2 (Supplementary Figure 9). The range of values is not only due to several genes being compared but also multiple transcripts per gene. For several transcripts there was a high degree of correlation between the expression profile for HvMYB61 and HvCslF6 (Person correlation coefficient r = 0.7 to 0.1) and HvTHX1 and HvCslF6 (Person correlation coefficient r = 0.6 to −0.1). Therefore, this analysis from an extensive range of samples and conditions provides additional evidence supporting the similar expression dynamics of HvMYB61 and HvCslF6.

To further investigate the expression profile of candidate MYB transcription factors we surveyed several smaller transcript datasets from barley that are independent of those described above. An RNA-seq dataset from cv. Morex across 16 barley tissues showed that HvMYB61 is expressed abundantly across the entire barley plant (Figure 5B). HvMYB61 expression was detected in 10-day seedling, shoot (10 cm), tiller (3rd internode), rachis, developing grain at 5 days post anthesis (DPA) and embryo; these are tissues where (1,3;1,4)-β-glucan is present. Similarly, HvTHX1 is also expressed in multiple tissues including developing grain (5 and 15 DPA). Expanding the RNA-seq expression data across barley grain development, an independent expression dataset based on six different spring barley cultivars (cvs. Hopper, Sloop, Extract, Taphouse, Alabama, and Pewter; Aubert et al., 2018; Matros et al., 2021) was used to retrieve the expression profiles of the candidate genes from early to mid-late grain development. HvMYB61 was mainly expressed during early grain development, from 7 to 11 DPA, whereas HvMYB46/83 expression remained almost undetectable from 7 to 20 DPA. In the RNA-seq dataset (developing grain) Pearson correlation coefficients with HvCslF6 were: 0.95 for HvMYB61 and 0.98 for HvTHX1 (showing abundant expression from 7 to 20 DPA; Supplementary Table 6). Our results showed that the positive correlation between HvMYB61 and HvCslF6 across different tissues and datasets, in the context of HvTHX1 (ortholog of BdTHX1, known activator of BdCslF6; Fan et al., 2018) and other HvCesAs (co-expressed and not co-expressed), are consistent with a regulatory effect of HvMYB61.

Based on expression profiles, our luciferase assay results and data from gene network analyses, we focused subsequent assays on the HvMYB61 transcription factor, since evidence suggested this gene may influence HvCslF6 expression. The candidate regulator HvMYB46/83 was ruled out based on low levels of expression from gene co-expression analyses. The candidate transcription factor HvMYB61 was cloned and transfected into barley protoplasts plus an empty vector control, respectively. After 24 h, protoplasts were harvested and HvCslF6 relative transcript abundance was measured. We observed a significant increase in HvCslF6 expression in the presence of HvMYB61, representing a fold change of 2.63 ± 0.64 in HvCslF6 expression compared to protoplasts transfected with an equal amount of empty vector (p-Value 0.035, Student’s t-test; Figure 5C). Taken together, our results from the functional characterization of the HvCslF6 promoter, TF identification and gene co-expression analyses, indicate that the HvMYB61 transcription factor may be a positive regulator of HvCslF6 expression in barley.

Barley transgenic lines (cv. Golden Promise) carrying a 3,000 bp upstream region from HvCslF6 start codon fused to a endoplasmic reticulum-targeted GFP (pHvCslF6:mGFP-ER) construct were generated via Agrobacterium-mediated transformation at The University of Adelaide following the protocol outlined in Burton et al. (2011). Transgenic lines (T₂) were screened using a Zeiss Axioskop 2 Plus microscope (Carl Zeiss Microscopy GmbH, Jena, Germany) equipped with an external UV light source (HBO 100). mGFP was excited at 488 nm and emission collected at 500–530 nm. For coleoptile screening, lines were grown in a petri dish with filter paper for 3–4 days in the dark and screened 5 days after gemination. Transverse section images were captured with an AxioCam 512. Subsequent to collection, images were processed with Zen Software v6.0, utilizing global adjustment tools only.

A 3,000 bp region upstream of the HvCslF6 start codon (cv. Morex) was retrieved from the barley genome explorer (Mascher et al., 2021) and used for TFBS prediction comparison using TRANSFAC^® v2014 (Matys et al., 2006), JASPAR v2020 (Fornes et al., 2020) and PlantPAN v3.0 (Chow et al., 2019). Motif over-representation analysis was carried out using HvCslF6 TFBS predictions from JASPAR and filtered based on ≥ 10.0 prediction score (maximum score = 15) to minimize false positive motif predictions. Upstream CslF6 sequences from wheat (Triticum aestivum), Brachypodium distachyon and rice (Oryza sativa) were retrieved from Ensembl Plants². TFBS comparisons across species were carried out within a 1,000 bp CslF6 putative promoter region using JASPAR (Fornes et al., 2020) and predicted motifs filtered based on ≥ 10.0 prediction score. The same filter was also applied for the TFBS analysis of HvCesA upstream regions (1,000 bp). Venn diagrams were created using InteractiVenn (Heberle et al., 2015).

Barley plants (cv. Golden Promise) were grown in a growth chamber with a 16-h light/8-h dark photoperiod at the Functional Genomics (FUNGEN) facility, at The James Hutton Institute, United Kingdom. After 3 weeks, 0.25 g of primary leaves were harvested for protoplast assays. Leaves were cut lengthwise and crosswise (∼20 mm × 10 mm) and peeled (epidermis removal from the abaxial side) in sterile conditions using a sharp razor, scalpel and tweezers under a Leica MZ6 stereo microscope.

The barley protoplast isolation protocol was mainly based on the protocol for preparation of Arabidopsis mesophyll protoplasts (Yoo et al., 2007) with modifications described in Supplementary Method 1. Briefly, barley leaves without the epidermis were transferred into 5 mL of the enzyme solution [2% w/v Cellulase Onozuka R-10 (Duchefa, Haarlem, Netherlands), 0.1% w/v Macerozyme R-10 (Duchefa, Haarlem, Netherlands), 0.1% w/v BSA (Sigma-Aldrich, St. Louis, MO, United States), 0.55 M Mannitol (Sigma-Aldrich, St. Louis, MO, United States), pH 5.7]. Leaves were incubated in a 6-well cell culture plate (Sigma-Aldrich, St. Louis, MO, United States) in the dark for 2 h at 28°C. Released protoplasts were mixed with one volume of W5 solution [154 mM NaCl, 125 mM CaCl₂, 5 mM KCl, 2 mM MES (Sigma-Aldrich, St. Louis, MO, United States, all), pH 5.7], filtered through a 70 μM gauze (Sigma-Aldrich, St. Louis, MO, United States), washed with 2.5 mL of W5 solution and centrifuged for 3 min at 70 × g, RT (Centrifuge 5810 R, Sigma-Aldrich, St. Louis, MO, United States). Protoplasts were gently re-suspended in 5 mL W5 solution and kept on ice for 30 min to allow cell sedimentation by gravity. Protoplasts were re-suspended to a final concentration of approximately 2 × 10⁵ protoplasts/mL in MMG solution (0.6 M Mannitol, 4 mM MES, 15 mM MgCl_2, pH 5.7). Fluorescein diacetate (FDA) was used to determine cell viability under a using a UV-light Zeiss Axioskop 2 Plus microscope and Canon EOS D100 digital camera.

For each transfection assay, 100 μL of protoplasts were mixed with 7 μg of construct DNA and 110 μL of freshly prepared 40% PEG-solution [0.4 M mannitol, 40% w/v PEG-4000 (Sigma-Aldrich, St. Louis, MO, United States), 0.1 M CaCl₂ (Sigma-Aldrich, St. Louis, MO, United States), pH 5.7]. After a 10 min incubation in the dark, protoplasts were washed twice with 1.5 mL of W5 solution and centrifuged at 200 × g for 5 min (Centrifuge 5415 D, Sigma-Aldrich, St. Louis, MO, United States). Protoplasts were re-suspended in 80 μL of W5 solution and incubated in the dark for 24 h, at room temperature. Three independent transfection assays were performed per construct and averaged as biological replicates.

The HvCslF6 (HORVU.MOREX.r3.7HG0698110) promoter deletion constructs were generated by PCR amplification of a 3,000 bp region upstream HvCslF6 start codon as follows: 1 μL cv. Golden Promise gDNA, 1.25 μL each primer (forward and reverse listed in Supplementary Table 7) at 10 mM, 5 μL 5X Phusion HF Buffer, 0.2 μL Phusion HF DNA polymerase (Thermo Fisher Scientific, Waltham, MA, United States) and 11.3 μL sdH₂O in a total volume of 20 μL. Phusion PCR reaction conditions: 98°C for 30 s, 35 cycles (98°C for 10 s, 55°C for 30 s, 72°C for 30 s/kb), 72°C for 5 min and 4°C hold. Promoter fragments were ligated by NotI restriction enzyme cloning into pGreenII 0800-LUC vector which contains luciferase and renilla reporter genes, enabled to be used in Dual Luciferase Reporter Assays (Promega, Madison, WI, United States). A total of seven PCR amplicons corresponding to HvCslF6 promoter fragments were NotI digested and purified with the QIAquick Gel Extraction Kit (Qiagen GmbH, Hilden, Germany). pGreenII-0800 vector was linearized in a similar digestion reaction using NotI with the addition of 1 μL TSAP (Promega, Madison, WI, United States) and purified as described for the pHvCslF6 fragments. Following this, T4 DNA ligations (Promega, Madison, WI, United States) were set up for each construct and 1 μL of the ligation reaction was used for E. coli transformation. Plasmids were purified with QIAprep Spin Miniprep Kit (Qiagen GmbH, Hilden, Germany) and confirmed by Sanger sequencing using M13F primer at the Genome Technology facility, The James Hutton Institute.

The Dual-Luciferase^® Reporter Assay System (Promega, Madison, WI, United States) was used to determine the promoter activity of HvCslF6 deletion constructs in barley protoplasts, 24 h after transfection. Luminescence results were obtained using a Varioskan LUX Multimode Microplate Reader (Thermo Fisher Scientific, Waltham, MA, United States) following the dual luciferase kit instructions. Protoplasts were lysed with 80 μL of a 1X passive lysis buffer and incubated for 15 min at room temperature. In parallel, 100 μL of LAR II (Luciferase Assay Reagent II) were pre-dispensed in a Corning^® 96-well opaque flat bottom microplate (Thermo Fisher Scientific, Waltham, MA, United States). For all assays, 20 μL of lysed protoplasts were transferred to the 96-well plate and incubated for 1 min, room temperature. After Firefly luciferase (Fluc) activity was measured (3 s per sample with a 1.5 s delay between wells), 100 μL of Stop and Glo Reagent were added enabling Renilla luciferase (Rluc) detection (1.5 s with a 0.5 s delay between wells). For each construct tested, three technical replicates were performed per transfection assay. A negative control, which included protoplasts transformed with empty vector and dual luciferase reagents, was used in each batch of transfections. Differences in dual luciferase reporter assays across HvCslF6 promoter deletion constructs were determined by one-way ANOVA followed by Tukey’s test using GraphPad Prism 8.4.2 Software (CA, United States). Pairwise comparisons were determined by two-tailed Student’s t-tests.

The protein sequences of cv. Morex v3 genome³ were used to identify R2R3 MYB family members. The sequences were uploaded to a webserver⁴ for automatic identification of MYB gene family members (Pucker, 2022) using default parameters. Sequences which were not classified as R2R3 MYB were removed for the phylogenetic analysis. Protein sequences of selected Arabidopsis and rice MYB proteins (MYB61, including OsMYB61a and OsMYB61b and MYB46/83) were added and the sequences aligned using ClustalW in MEGA (Kumar et al., 2018). Unreliable positions in the alignment were removed using BMGE (Block Mapping and Gathering with Entropy; Criscuolo and Gribaldo, 2010) v2.0. Model selection for amino acid substitution was done in MEGA⁵ resulting in LG model plus gamma distribution with invariant sites (LG+G+I) as the best choice. The unrooted phylogenetic tree was build using Maximum Likelihood with 100 bootstrap replications to estimate bootstrap support. Only bootstrap values above 0.5 are displayed.

Gene expression values from 808 individual samples (Supplementary Table 4) were obtained and are available from Milne et al. (2021). In short, the reads were downloaded from the sequence read archive⁶ and mapped against BaRTv1 (Rapazote-Flores et al., 2019) using Salmon (Patro et al., 2017). Expression of the isoforms were added up to obtain gene expression. For network construction low expressed genes were removed by filtering for a TPM (transcript per million) of above 5 in at least a third of the samples. This condensed the number of expressed genes from 60,444 to 17,678. Weighted gene correlation network analysis (WGCNA; Langfelder and Horvath, 2008) was used for network construction. Euclidean clustering beforehand highlighted one sample from E15.MicrosporeEmbryogenic (SRR6433018) as outlier. Therefore, this sample was removed, leaving 807 samples for further analysis. The soft-thresholding power was set to 10 to achieve scale-free topology. Network analysis was carried out using the blockwise module method (Langfelder and Horvath, 2008), with a power of 10, TOMType = unsigned, corType = bicor, networkType = signed hybrid, minModuleSize = 30, maxPOutliers = 0.05, mergeCutHeight = 0.15, deepSplit = 3. Visualization of the HvCslF6 network was performed in Gephi⁷. Size and coloring of the nodes was determined by degree centrality (the number or proportion of other nodes linked to a specific node). Edge coloring was based on weight. The subnetwork was extracted by filtering for nodes connected to HvCslF6 (BART1_0-p53331). For this we applied an edge weight which is a threshold, whereby connections between gene represent by values greater or equal to 0.188 are correlated.

A nested PCR was used to amplify HvMYB61 (HORVU.MOREX.r3.1HG0018590) containing: 1 μL of cv. Golden Promise protoplast cDNA (1:10 diluted in sdH₂O), 1.25 μL of each primer (forward and reverse listed on Supplementary Table 7) at 10 mM, 5 μL of 5X Phusion HF Buffer, 0.2 μL of Phusion HF DNA polymerase (Thermo Fisher Scientific, Waltham, MA, United States) and 11.3 μL of sdH₂O in a total volume of 20 μL. Phusion PCR reaction conditions are described above. A 1 μL aliquot of the PCR mix was used as a template to attach Gateway^® attB sites, 8x His-tag and thrombin cleavage site at the C-terminus in a similar Phusion HF PCR reaction and gel purified. HvMYB61 was cloned into pDONR207™, linearized with EcoRI by a BP reaction containing: 3 μL of the purified attB-PCR product (34 ng/μL), 1 μL of linearized pDONR207™ (50 ng/μL), 2 μL 5X BP Clonase™ Reaction Buffer (Thermo Fisher Scientific, Waltham, MA, United States) and 4 μL of TE buffer, pH 8.0 in a total volume of 10 μL. The BP reaction was incubated at RT overnight and stopped by adding 2 μL of Proteinase K and incubated for 10 min at 37°C. 1 μL of the BP reaction was used for E. coli transformation. Positive colonies were transferred to 4 mL of LB medium with gentamycin (100 μg/μL) and incubated overnight, 37°C and shaking at 230 rpm. Plasmids were purified using the QIAprep Spin Miniprep Kit (Qiagen GmbH, Hilden, Germany) and confirmed by Sanger sequencing. HvMYB61 was transferred to pBract214m-HSPT plant expression vector by LR reaction containing: 2 μL of EcoRI linearized pDONR207™-HvMYB61 (40 ng/μL), 5 μL of pBract214m-HSPT (30 ng/μL), 1 μL TE Buffer, pH 8.0 and 2 μL LR Clonase™ II enzyme mix (Thermo Fisher Scientific, Waltham, MA, United States) in a total volume of 10 μL. The reaction was incubated at 25°C for 3 h and stopped by adding 1 μL of the Proteinase K (37°C for 10 min). After E. coli transformation, positive colonies were grown on LB medium with kanamycin (100 μg/μL) for plasmid purification. pBract214m-HSPT-HvMYB61 construct sequence was confirmed as described above using the ZmUbi forward primer (Supplementary Table 7).

Total mRNA was extracted from ∼2 × 10⁵ protoplasts (100 μL) using TRIzol^® Reagent (Thermo Fisher Scientific, Waltham, MA, United States) according to the manufacturer’s instructions. RNA concentration and purity were measured using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, United States). Total mRNA to cDNA conversion was performed using cDNA EcoDry™ Premix (Takara, Kyoto, Japan). For each reaction, 100 ng of total mRNA was added to the lyophilized master mix following the manufacturer’s instructions. Quantitative real-time PCR (qRT-PCR) was performed in a StepOne Real-Time PCR machine (Thermo Fisher Scientific, Waltham, MA, United States) using PowerUp SYBR Green Master Mix (Thermo Fisher Scientific, Waltham, MA, United States) to determine HvCslF6 relative transcript abundance. Three independent experiments were performed, each of them included: HvMYB61 PEG-mediated transfection, RNA extraction, cDNA synthesis and qPCRs, respectively. Three replicate qRT-PCR reactions were performed for each cDNA sample (technical replicates) including a negative control (sdH₂O as template). Each qRT-PCR reaction contained: 2 μL protoplast cDNA (1:2 dilution), 5 μL SYBR Green, 1 μL each forward and reverse primer at 4 mM and 1 μL of sdH₂O in a total volume of 10 μL. Primer sequences for HvCslF6 and qRT-PCR conditions were used as in Burton et al. (2008) and described in Supplementary Table 7. Relative HvCslF6 gene expression was normalized to α-Tubulin, Gapdh and Hsp70 housekeeping genes and calculated using the 2^–ΔΔCT method (Vandesompele et al., 2002) for multiple control genes. Statistical differences in HvCslF6 expression were determined by Student’s t-test compared to WT protoplasts (empty vector) using GraphPad Prism 9.1 Software (CA, United States).