Global analysis of SBP gene family in Brachypodiumdistachyon reveals its association with spike development

SQUAMOSA-promoter binding like proteins (SBPs/SPLs) are plant specific transcription factors targeted by miR156 and involved in various biological pathways, playing multi-faceted developmental roles. This gene family is not well characterized in Brachypodium. We identified a total of 18 SBP genes in B. distachyon genome. Phylogenetic analysis revealed that SBP gene family in Brachypodium expanded through large scale duplication. A total of 10 BdSBP genes were identified as targets of miR156. Transcript cleavage analysis of selected BdSBPs by miR156 confirmed their antagonistic connection. Alternative splicing was observed playing an important role in BdSBPs and miR156 interaction. Characterization of T-DNA Bdsbp9 mutant showed reduced plant growth and spike length, reflecting its involvement in the spike development. Expression of a majority of BdSBPs elevated during spikelet initiation. Specifically, BdSBP1 and BdSBP3 differentially expressed in response to vernalization. Differential transcript abundance of BdSBP1, BdSBP3, BdSBP8, BdSBP9, BdSBP14, BdSBP18 and BdSBP23 genes was observed during the spike development under high temperature. Co-expression network, protein–protein interaction and biological pathway analysis indicate that BdSBP genes mainly regulate transcription, hormone, RNA and transport pathways. Our work reveals the multi-layered control of SBP genes and demonstrates their association with spike development and temperature sensitivity in Brachypodium.


Results
Identification and characterization of SBP-box genes in B. distachyon. In this study, we identified 18 SBP genes in B. distachyon and designated as BdSBP. BdSBP family members were named according to the closest homologs present in wheat, barley or rice. Details of SBP gene family in Brachypodium are given in Table 1. Brachypodium SBP genes encode proteins ranging from 177 (SPL7) to 1,110 (SPLl4) amino acids (aa) in length and from 122 kda (SBP15) to 12 kda (SBP23A) in molecular weight. The number of exons ranged from 1 to 11 and isoelectric point (pI) was from 5 to 10. The 18 BdSBP genes were located on all 5 chromosomes (chr), with maximum number of BdSBP genes detected in chr 3 of B. distachyon (Table 1).

Phylogenetic analysis and gene duplication in BdSBP genes.
A phylogenetic tree was constructed using conserved SBP domain sequences of SBP proteins from Brachypodium, wheat, barley, rice and Arabidopsis (Fig. 1A,B). A total of 79 SBP proteins from different plant species including 18 from rice, 10 from wheat, 17 from barley and 16 from Arabidopsis were used for phylogenetic analysis. SBPs clustered into 8 groups (G1-G8), with AtSBP3/4/5/6 as ungrouped members. Each group contained at least one SBP protein from Brachypodium. As anticipated, BdSBPs exhibited closer relationship with the SBP proteins from barley and wheat as compared to rice and Arabidopsis. Group 1 and group 5 contained maximum number of BdSBPs, where SBP proteins from barley and wheat were also grouped. Moreover, gene duplication analysis among BdSBP genes identified 9 putative paralogous gene pairs in the Brachypodium genome ( Fig. 2A,B). Divergence time for duplicated BdSBP genes was estimated from Ka and Ks values and their ratios. The dates of duplication events (T) were calculated using Ks values through the formula T = Ks/2λ × 10 −6 (millions of year, Mya). The λ = 6.5 × 10 −9 substitutions per synonymous site per year was assumed as universal clock-like rate for Brachypodium distachyon. For BdSBP1 and BdSBP6 gene pair, Ka and Ks values were 0.60 and 2.14, respectively and their ratios 0.28 imply their evolution under purifying selection. Similarly, the ratio (0.30) of Ka and Ks values for BdSBP16 and BdSBP18 gene pair highlights purifying selection. Purifying selection also called negative selection, influence genomic diversity in natural populations. It eliminates the changes that produce deleterious effects on the fitness of the host. The frequency distributions indicate that SBP genes in Brachypodium went through a large-scale duplication event ranging from 55 to 164 million years ago (mya). SBP gene paralogs were located on same as well as different chromosomes, indicating that expansion of Brachypodium SBP genes was both, tandem as well as segmental/ block duplication during evolution. ) shared a similar type of motif structure. Some motifs were found to be specific to one or two groups of BdSBP proteins. Motif 6 that encodes miR156 target sequence was present in all miR156 targeted BdSBP proteins. Whereas motif 10 and motif 4 were found in group 1 BdSBP and group 3 and 4 BdSBP proteins respectively. To predict possible functions of BdSBP genes, we also performed gene ontology (GO) term enrichment analysis (Fig. 3C) Resul ts) from 9 different tissues and organs (leaf, early inflorescence, emerging inflorescence, anther, pistil, seed 5 days after pollination, seed 10 days after pollination, plant embryo and endosperm) was mined to understand the dynamics of BdSBP genes expression (Fig. 4). The expression profile of BdSBPs was grouped into three clusters. Higher expression of BdSBP genes of cluster 1 was observed in early inflorescence, emerging inflorescence and pistil tissues. Many BdSBP genes of cluster 1 were also expressed in anther, plant embryo, developing seeds (5 and 10 days after pollination), leaf, and endosperm tissues, implying their significant role throughout the Brachypodium plant development, especially in spike architecture. The BdSBP genes of cluster 2 (BdSBP7 and BdSBP1) either lacked expression in any tissue (BdSBP7) or poorly expressed (BdSBP1) in pistil, leaf, and developing seeds (seed 5 and 10 days after pollination), and endosperm. The BdSBP genes from cluster 3 were found to be expressed mainly in early and emerging inflorescence. Expression profiles of BdSBP genes indicate their involvement in the reproductive units of B. distachyon. post-transcription of BdSBP genes is regulated by miR156 and alternative splicing (AS). The cDNA sequences of BdSBP genes were searched for putative target sites of Brachypodium miRNAs (Fig. 5A). Ten www.nature.com/scientificreports/  www.nature.com/scientificreports/ BdSBP genes are found to be the target of miR156. Out of these, 8 contain miR156 complementary sequences in their coding regions. However, in other two genes, BdSBP1 and BdSBP13, miR156 target site was found in their 3′-UTR and 5′-UTR regions, respectively. The BdSBP gene family undergoes AS which specifically targets miR156 regulated BdSBP (Fig. 5B). The number of splice isoforms for each BdSBP genes was derived from plant Ensembl database. Splice variants from Ensembl gene are compared to generate an inclusive list of elementary alternative splicing events. The range of splice isoforms produced by BdSBPs was between 2 and 7. Most of the splice variants of BdSBP genes possess miR156 target site except BdSBP1, which has 4 splice variants and only one contains miR156 target site.
Organ specific differential accumulation of BdSBP genes is regulated by miR156. Three BdSBP genes (BdSBP3, BdSBP17 and BdSBP23) were analyzed for miR156 mediated transcript degradation by 5′-RLM-RACE (Fig. 5C,D). Additionally, to observe miR156 mediated cleavage pattern, we also constructed cDNA libraries from leaf and spike. Interestingly, BdSBP genes were highly degraded by miR156 in leaf as compared to spike tissue. We hypothesized that this differential degradation of BdSBP genes in leaf and spike tissues might be connected with their expression patterns in these tissues. To validate this, we performed semi-quantitative RT-PCR of several potential BdSBP genes on the basis of in silico expression data (Fig. 5E). Our data indicate that miR156 targeted BdSBP genes indeed expressed poorly in the leaf and abundantly in the spike, confirming our hypothesis. The miR156 non-targeted genes (BdSBP9 and BdSBP15) expression was constant in both the leaf and the spike, suggesting no effect of miR156 on these genes. Furthermore, to map the miR156 cleavage site in BdSBPs transcript, the 5′-RLM-RACE products were cloned and sequenced. Data indicates that miR156 cleaves between 9 and 10th nucleotide of 5′ site of BdSBPs transcript, except BdSBP3 where cleavage site was found BdSBP genes are involved in complex regulatory network and pathways. BdSBPs co-expressed genes were investigated using publicly available large-scale co-expression database (www.gene2 funct ion.de), and MapMan (https ://mapma n.gabip d.org) ontology term enrichment to study their roles in different biological pathways ( Fig. 6A-B). Around 710 co-expressed genes were found to be associated with 15 members of BdSBP family (Supplementary Table S5). The MapMan ontology of the co-expressed genes suggests that 22 out of 35 of the major biological classes have at least one of the BdSBP family members (Fig. 6B). Cell, development, transport, hormone metabolism, secondary metabolism, stress, lipid metabolism, cell wall, DNA, RNA, protein and signaling were major biological processes in which co-expressed genes of BdSBP family members were involved. Some other BdSBP family members and their co-expressed genes were enriched in photosynthesis, major CHO metabolism, fermentation, oxidative pentose pathway, mitochondrial electron transport and amino acid metabolism. However, BdSBP co-expressed genes were not augmented in the C1-metabolism, microRNA, polyamine metabolism, nucleotide metabolism, S-assimilation, N-metabolism, glycolysis and minor CHO metabolic pathways. In addition, a network of protein-protein interaction of B. distachyon proteins was developed using STRING database. This database predicts interactions based on experimentally determined, predicted, text mining, coexpression, gene fusion, gene neighbourhood etc. A total of 39 interactive proteins were found (confidence value = 0.5) for 9 of BdSBP proteins, which were based on either predicted interactions or text mining (Fig. 7A, Supplementary Table S6). Protein annotation reveals that BdSBP proteins might interact with MYB33, PHABU-LOSA (PHB), Homeobox TF family, Growth regulating factor 5 (GRF5), Heat shock TF, ZnF C2H2, F-box TF, NBS-LRR, protein kinase family protein, ankyrin repeat protein, protein kinase family protein, chlorophyll a-b binding protein, DCL1, DCL2 and DCL3 proteins. The BdSBP7 was the only B. distachyon protein that interacts with DCL2 and DCL3 proteins. The MapMan term ontology of interactive protein partners of the BdSBPs indicates that 10 out of 35 proteins of the major biological terms were enriched by at least one of the BdSBP protein network (Fig. 7B). These biological processes were linked to development, RNA, photosynthesis, cell wall, protein, transport, signaling, cell cycle and stress.
BdSBP genes express differentially during variable temperature conditions. In order to advance our knowledge about the molecular mechanism controlling heat stress in Brachypodium, we examined the tran-  BdSBP genes regulate spike development and flowering. In silico expression analysis revealed higher expression of BdSBP genes during spike emergence and in early inflorescence development (Fig. 4). Therefore, the transcript abundance of 9  (Fig. 8B). Expression of BdSBP9 was slightly lower at 15DAH as compared to 7 and 24DAHs. To ensure the reproductive success, flowering is the critical stage of plant reproduction, which is mainly regulated by gibberellin, vernalization, photoperiod and autonomous pathways. Vernalization promotes the flowering in alpine species and its molecular mechanism has been investigated in Arabis alpina and A. thaliana plants (Bergonzi et al. 51 ). Therefore, to understand the genetic control of vernalization response in grasses, we analysed the expression pattern of BdSBP genes in Brachypodium. Transcript abundance of several BdSBP genes was compared in the rapid flowering (Bd21) and delayed flowering (Bd1-1) accessions of Brachypodium under vernalization and non-vernalization conditions (Fig. 8C,D; Supplementary Fig. S4). We observed that Bd21 accession flowered rapidly under non-vernalised condition, whereas Bd1-1 lacked flowering until maturity.  To further confirm the function, T-DNA mutant for BdSBP9 gene was obtained from JGI ( Fig. 9A-C). The Bdsbp9 mutant has a T-DNA insertion in the first exon of BdSBP9. Electron microscopy indicated that different patterns of lignification in the wild type as compared to Bdsbp9. Wild-type patterns were straighter, with no circular patches whereas, Bdsbp9 patterns are less uniform, with some circular patches. Further, we investigated the promoter region (1000 bp upstream of initiation codon) of the co-expressed genes of BdSBP9 (Fig. 9D). Data indicate that 92% of the co-expressed genes contain GTAC motif, a specific binding site for SBP genes. The expression of one of the interacting partners matches with BdSBP9 expression pattern (Fig. 9E).

Discussion
SBP/miR156 genetic circuit controls the transition of vegetative to reproductive phase change in Arabidopsis 14,52 . Owing to the importance of SBP genes, we conducted the first-ever genome-wide identification of this gene family in Brachypodium and discovered 18 BdSBP genes ( Table 1). The number of SBP genes in Brachypodium were similar to the SBP genes in barley (17), B. luminifera (18), rice (19) and Arabidopsis (17), but was smaller in comparison to soybean (41), moso bamboo (32) and P. trichocarpa (28), suggesting that SBP genes were evolved in a species specific manner and underwent different gene duplication events. On the basis of phylogenetic analysis BdSBP genes were divided into eight (G1-G8) groups (Fig. 1A). BdSBP genes grouped closely with HvSPLs and TaSPLs, suggesting that these SBP genes possibly diverged from a common ancestor. The DNA binding SBP domain binds to the promoter regions of its target genes containing TNCGT ACA A consensus nucleotide sequence with GTAC as a core motif 53 . Two zinc ion binding motifs Cys3His1 and Cys2His1Cys1 at N terminus and a nuclear localization signal (NLS) at C terminus were found in BdSBP proteins (Fig. 1B). SBP genes share similar gene structures within their same phylogenetic group as mentioned previously in barley 13 , rice 11,54 , and tomato 55 . Gene duplication events are key to evolution and gene expansion which produce many paralogous gene pairs 56 . Additionally, gene duplication also assists organisms to cope up with different environmental conditions during growth and development 57 . In order to study gene duplication, we estimated the Ka/Ks ratio for each duplicated genes using Plant Genome Duplication Database and PLAZA 4.0 ( Fig. 2A,B), which suggests that BdSBP genes underwent duplication event ~ 74 to 164 mya. The Ka/Ks ratio of > 1 shows that the gene has experienced positive selection, = 1 indicate neutral selection and < 1 indicates purifying selection or negative, respectively. Based on Ka/Ks ratio 56 ; the BdSBP gene pairs which were ranged from 0.2 to 0.7 suggesting that these genes were duplicated under purifying selection. Also, BdSBP genes shared similar intron/exon structures within the same phylogenetic groups (Fig. 3A). Additionally, most of the BdSBPs from the same phylogenetic groups possess similar motifs (Fig. 3B). Consequently, the genes in the same phylogenetic group might have similar roles in Brachypodium, which have been supported by gene ontology terms of BdSBP genes (Fig. 3C). In addition to conserved BdSBPs motifs, several unique group-specific motifs were observed, such as motif 4, 9 www.nature.com/scientificreports/ and 10 in group 1 and motif 8 in group 8. These specific motifs might be important for specified roles of BdSBP genes, and their functional differentiation could arise during evolution of different lineages. All BdSBP genes expressed substantially during early and emerging inflorescence development except, BdSBP1 and BdSBP7, implying their role in inflorescence development of Brachypodium (Fig. 4). Most of the BdSBP genes except, BdSBP7/8/17/21, expressed significantly in pistil whereas BdSBP3/6/11/ 14/15/16 expressed highly in the anther, suggesting their role during reproduction. However, BdSBP6/9/15 were constitutively expressed in all the 9 tissues. Previously, differences in expression profiles of miR156 targeted and non-targeted SBP genes have been reported in barley 13 , Brassica napus 58 , Betula 59 and soybean 12 . Importantly, miR156 targeted BdSBP genes showed differential expression pattern and most of the miR156 non-targeted BdSBPs showed constitutive expression profiles in Brachypodium. Post-transcriptional regulation of SBP genes through miR156 has been considered as the key process for the functionality of these genes 27,29,60 . A total, 11 of 19 SBP genes in rice 11 , 10 of 17 SBPs in Arabidopsis 52 , 7 of 17 SBPs in barley 13 , and 18 of 28 SBPs in Populus 8 have been identified as targets of miR156. In our study, miRNA target prediction revealed that 10 BdSBP genes are regulated by miR156 (Fig. 5A). A total 8 (BdSBP3,-11,-13A,-14,-16,-17,-18, and 23) of 10 miR156 targeted BdSBP genes contained miR156 complementary sequence in their coding region whereas, BdSBP1 and BdSBP13 contained target site in 5′ and 3′UTRs, respectively. Thus, miR156 targets BdSBP1 and BdSBP13 along with other BdSBP genes will be unable to perform downstream roles. This phenomenon of binding of miRNAs to their complementary sequences in the coding sequences or un-translated regions of target genes to inhibit gene function either by transcript cleavage or deadenylating has also been reported elsewhere (Rhoades et al. 61 ).
In humans ~ 95% and in Arabidopsis > 60% of multi-exonic genes undergo AS 62 . Meanwhile, we noticed that miR156 targeted BdSBP genes produced different splice variants via. AS (Fig. 5B). AS generally produces transcripts with premature stop codon which are degraded in cytoplasm by non-sense-mediated decay (NMD) pathway 63 . Splice variants produced by AS generally exhibit spatiotemporal or environmental condition-specific expression patterns 63 . Our experiments in Brachypodium showed that BdSBP gene products are degraded at higher level by miR156 in the leaf as compared to the spike (Fig. 5C,D), resulting into higher transcript abundance of miR156 targeted BdSBP genes in young spike as compared to leaf ( Fig. 5E; Supplementary Fig. S1). However, expression of miR156 non-targeted BdSBP9 and BdSBP15 genes was constitutive in these tissues. This confirms that miR156 negatively regulates SBP genes in Brachypodium and is consistent with previous findings in Arabidopsis, rice, tomato and wheat 11,29,35,52,64 . Taken together, these results suggest that miR156 in conjunction with AS regulates the transcriptome dynamics. www.nature.com/scientificreports/ SBP-correlated gene network and interactome analysis revealed that SBP genes function by regulating other families of transcription factors and membrane transport proteins, and are involved in the metabolism of glucose, in-organic salts and ATP production in Arabidopsis 65 . Therefore, considering the significance of Brachypodium as a model plant for developmental biology of triticeae crops, we examined the co-expression and MapMan biological pathways (Fig. 6A,B; Supplementary Table S5). MapMan terms enrichment analysis showed that BdSBP genes perform their function by regulating transcription, protein, signalling, transport and development related biological pathways (Fig. 6B). The co-expression network contains mainly transcription factors, hormones (auxin, brassinosteroid, ethylene and gibberellin) responsive genes, cell wall biogenesis related genes and transporters, implying their roles in development as well as cell wall biogenesis of Brachypodium. Existence of CSLF3 and MYB TF indicate that BdSBP genes might be involved in secondary wall synthesis in Brachypodium 66 . Studying the protein-protein interaction network represents gene functions crucial to plant physiology, pathology, and growth 67 . Protein-protein interaction at the molecular level might be important in transcription regulation, post-transcriptional modification, cytoskeleton assembly, phosphorylation, acetylation, transporter activation and others 68 . Previously, it was found that IPA1 (OsSPL14), an important factor which controls plant architecture interacts with D53 protein (DWARF53) in-vivo and in-vitro 69 . Recently, OsSPL14 protein has been shown to be associated with disease and yield in rice by phosphorylation and non-phosphorylation of Ser 163 amino acid respectively during Magnaporthe oryzae fungal infection 5 . In our study 39 interacting proteins with 9 BdSBP proteins were identified (Fig. 7A,B; Supplementary Table S6). These interacting proteins mainly belonged to bZIP, Homeobox, MYB33, ZnF_C2H2, F-box and heat shock transcription factor families, Dicer-like proteins and protein kinases. These interacting protein partners have been involved in the regulation of the biological pathways including development, RNA, protein, stress, photosynthesis and cell wall, implying the diverse roles of BdSBP proteins in Brachypodium growth and development.
The grain development and filling of Brachypodium spikelet are completed (dry) in 50 days and has been classified into three stages namely-embryo and endosperm development [0-14 days after fertilization (DAF)]; maturation (14-36 DAF) and desiccation (36-50 DAF) stages 70 . Higher expression of BdSBP1/-3/-16/-17/-18/-21 and 23 at spikelet initiation stage as compared to the maturation stage, might be key to early spikelet development in Brachypodium (Fig. 8A,B). Further, BdSBP9 and BdSBP15 genes exhibited constitutive expression pattern during embryogenesis and maturation stages, suggesting their importance for these stages. Plants bear flowers at a certain time of reproductive phase which is mainly regulated by SBP/miR156 pathway 29 . As plants grow older, the level of SBP genes increases while miR156 abundance declines. Previously, it was reported that higher production of SBP genes ensures flowering in response to cold in the model perennial Arabis alpina accession Pajares 51 . It has been reported in Cardamine flexuosa that SBP/miR156 pathway plays a key role in flowering through integrating age and vernalization pathway 71 . The SPL/miR156 module has been known to be a key component for flowering phases 1,52 . Involvement of SBP genes in the control of flowering time of B. distachyon accessions Bd21 and Bd1-1 under vernalization condition (Fig. 8C,D) suggest that BdSBP1 and BdSBP3 potentially involved in this. This result positively supports the previous study about sensitivity of certain SBP genes to vernalization in older plants 51 . Cereals inflorescence (spike) architecture is one of the main determinants of their yield. In rice, SBP genes have been reported as an important regulator of plant architecture. Overexpression of OsSPL14, present on the IPA1 (ideal plant architecture)/WFP (wealthy farmer's panicle) QTL, decreased tiller branching but increased panicle branching and grain weight 18,36 . Likewise, OsSPL7, OsSPL13 OsSPL16 and OsSPL17 also regulate grain size, shape and yield in rice 16,37,38 . In our study, the Bdsbp9 mutant showed abnormal spike and delayed flowering, implying its role in spike development (Fig. 9A-E). In Arabidopsis, miR156/ SPL module confers thermotolerance at reproductive stage [ 24 , 72 ). Our study also indicates that BdSBP genes contribute thermotolerance during spike development in Brachypodium. Interestingly, differential expression of BdSBP genes in the developing spike under variable temperatures was not been associated with ploidy level in Brachypodium genome as described previously 49 . However, specific expression of these genes in response to high temperature in tetraploid genome, B. hybridum, probably induced by interaction of B. distachyon and B. stacei genomes ( Fig. 7C; Supplementary Fig. S2). Overall, our study revealed that altering the expression pattern of BdSBP genes may provide an important tool-box for the genetic improvement of the cereal crops.
Gene structure and phylogenetic analysis of BdSBP genes. The exon/intron structure of each BdSBP gene was predicted through gene structure display server program (https ://gsds.cbi.pku.edu.cn/index .php) by comparing their coding and genomic sequences. The TAIR (https ://www.arabi dopsi s.org/index .jsp) was used to obtain the Arabidopsis SBP sequences and rice genome annotation project database was used to obtain the rice SBP genes sequences. SBP sequences of wheat were obtained from a previous study 74  www.nature.com/scientificreports/ (https ://itol.embl.de). SBP domain sequences were aligned using MUSCLE tool followed by Gblocks curation utilities and maximum likelihood method was used to construct the phylogenetic tree using PhyML software (https ://www.phylo geny.fr).
Motif identification, miR156 target site prediction and alternative splicing event analysis. The MEME 4.11.0 tool (https ://meme-suite .org/tools /meme;) was used to search for conserved motifs within BdSBP proteins by using default settings, except that the maximum number of motifs to find was 10, the maximum width was 50 and the minimum width was 6. The sequence logo of the Brachypodium SBP domain was created with an online available WebLogo3 platform (https ://weblo go.three pluso ne.com/). The cDNA sequences of BdSBPs were subjected to psRNATarget tool (https ://plant grn.noble .org/psRNA Targe t/?funct ion) to predict the putative target sites of miR156. The Ensemble database (https ://plant s.ensem bl.org/Brach ypodi um_dista chyon /Info/Index ) was used to obtain the information on alternative splice events for each BdSBP gene (Supplementary Table S5).
Gene expression analysis of BdSBPs. The log2-transformed fragments per kilobase per million fragments measured (FPKM) values were used to study the expression of BdSBPs in nine tissues as described 13,74 .
A heat map of the expression of BdSBPs was generated by the average hierarchical clustering method using the MeV tool (https ://www.tm4.org/mev.html).
Co-expression, protein-protein interaction and gene duplication analysis. The co-expressed genes for BdSBP members were identified using online PlaNet (https ://arane t.mpimp -golm.mpg.de/index .html) tool 76 . PlaNet uses the Pearson correlation coefficient (PCC) and constructs a co-expression network, with PCC cut-off to 0.7 77 . Further, a highest reciprocal rank (HRR) co-expression network with standard edge cutoff of 30 was used. Additionally, a heuristic cluster chiseling algorithm (HCCA), which is optimized for HRRbased networks was used with standard parameters (stepsize = 3). The protein-protein interaction network was identified using STRING database (https ://strin gdb.org/cgi/input .pl?sessi onId=A92xE G08sQ Ek&input _page_show_searc h=on), which contained information from various datasets such as; gene coexistence, protein-protein interactions, gene fusion and co-expressed genes to calculate the semantic links between proteins 78  Genomic DNA and RNA isolation. Leaves from Brachypodium distachyon plants were collected and DNA isolation was performed using cetyl-trimethyl-ammonium bromide-based (CTAB) extraction method as described elsewhere 82 . PCR-based genotyping was performed using primers following recommendations from Joint Genome Institute 80 . The spectrum plant total RNA Kit (Sigma-Aldrich, St. Louis, MO, USA) was used for RNA isolation following the manufacturer's protocol. The RNA integrity and purity of all samples were verified on a Nanodrop ND-1000 (Nanodrop Technologies, Wilmington, DE, USA). Prior to synthesizing cDNA, RNA samples were treated with DNase I to remove genomic DNA contamination (Invitrogen, USA). The reaction mixtures were incubated at 23 °C for 15 min and after that 1 µl of 25 mM EDTA was added to each sample.
First strand cDNA synthesis and quantitative real-time PCR (qRT-PCR) analysis. First strand cDNA was synthesized from 2 µg total RNA sample using AffinityScript QPCR cDNA Synthesis Kit (Agilent technology, Canada). The qRT-PCR was run on Mx30005p qPCR system (Stratagene, USA) in a 20-μl volume containing 5 µM gene-specific primers, 1 μl diluted cDNA, and 10 μl Brilliant III Ultra-Fast SYBR Green QPCR Master Mix (Agilent, USA). Two biological and three technical replicates were used in all the experiments. The 2 −ΔΔCq method was used to quantify the relative level of gene expression (Livak and Schmittgen 2001). The genespecific primers for BdSBP genes used in semi-quantitative RT-PCR and qRT-PCR are listed in Supplementary  Table S7. PCR was performed in a 20 μl volume using GoTaq Green master mix (Promega, USA). BdUBC18 (Ubiquitin-conjugating enzyme 18) was used as a reference gene for different developmental stages and SamDC (S-adenosyl methionine decarboxylase) was used for heat stress 83 .

RNA ligase-mediated modified 5′ rapid amplification of cDNA ends (RLM-RACE). The miR156
mediated cleavage site in the BdSBP transcript was mapped by using First-choice RLM-RACE kit (Ambion, Austin, TX, USA). Total RNA was isolated from leaf and 7 days old spike (after heading). Without pre-treatment, 1 μg total RNA was ligated to the 5′ RACE RNA adapters. The M-MLV reverse transcriptase enzyme and 18mer oligo dT were used to reverse transcribe the adapter-ligated RNA. Primary and secondary nested PCR was carried out using 5′ RACE gene-specific outer and 5′-adapter outer primers and 5′ RACE gene-specific inner and 5′-adapter inner primers. The primer sequences used in nested PCR are listed in Supplementary Table S5. The 5′ RACE PCR amplified fragments were gel extracted and cloned into pGEMT easy vector. Further, clones were confirmed by EcoRI restriction analysis and Sanger sequencing.