Genome-Wide Identification of RNA Silencing-Related Genes and Their Expressional Analysis in Response to Heat Stress in Barley (Hordeum vulgare L.)

Barley (Hordeum vulgare L.) is an economically important crop cultivated in temperate climates all over the world. Adverse environmental factors negatively affect its survival and productivity. RNA silencing is a conserved pathway involved in the regulation of growth, development and stress responses. The key components of RNA silencing are the Dicer-like proteins (DCLs), Argonautes (AGOs) and RNA-dependent RNA polymerases (RDRs). Despite its economic importance, there is no available comprehensive report on barley RNA silencing machinery and its regulation. In this study, we in silico identified five DCL (HvDCL), eleven AGO (HvAGO) and seven RDR (HvRDR) genes in the barley genome. Genomic localization, phylogenetic analysis, domain organization and functional/catalytic motif identification were also performed. To understand the regulation of RNA silencing, we experimentally analysed the transcriptional changes in response to moderate, persistent or gradient heat stress treatments: transcriptional accumulation of siRNA- but not miRNA-based silencing factor was consistently detected. These results suggest that RNA silencing is dynamically regulated and may be involved in the coordination of development and environmental adaptation in barley. In summary, our work provides information about barley RNA silencing components and will be a ground for the selection of candidate factors and in-depth functional/mechanistic analyses.


Introduction
Barley (Hordeum vulgare L.) is the fourth most important crop in the world in terms of cultivated area (50 MHa) and grain yield produced (140 Mt) [1]. Barley grain has a high nutritional component content [2] and is used as a human food and animal feed. Grain germination speed and consistency of endosperm cell wall breakdown into fermentable sugars makes it a primal raw material for the production of alcoholic beverages or biofuels [1, 3,4]. Barley straw contains a high level of lignocellulose, that may also be used as a form of renewable energy [4]. Beyond the obvious economic importance, studying barley is useful for understanding and increasing crop resilience.
Barley domestication over ten thousand years ago offers strong evidence of its capacity to acclimate to a wide range of environments; indeed, barley was proposed as a good model for adaptation studies [5]. Although relatively tolerant to abiotic stresses among cereal crops, heat stress during reproductive phase negatively affects barley grain yield and quality [6,7]. Barley genetic variation may be employed for the development of efficient strategies to enhance its productivity under diverse climatic conditions. There are over 400,000 barley accessions available [8]. Germplasm availability and expansion may be critical for the sustained high yield of crops under climate change and global warming conditions. RNA silencing is an evolutionary conserved sequence-specific gene-inactivation pathway. Originally evolved as an immunity system [9,10], RNA silencing acquired several cellular roles including developmental regulation, stress responses, or chromatin organisation. RNA silencing may act both at transcriptional (Transcriptional Gene Silencing, TGS) or post-transcriptional (Post-Transcriptional Gene Silencing, PTGS) levels [11,12]. First discovered in plants, it was later described in many other eukaryotic organisms [13][14][15].
AGOs, alongside the incorporated sRNA (and other associated proteins), form the RNAinduced Silencing Complex (RISC). AGOs are the executor components, while sequence-specificity of RISC activity is provided by the sRNA. Domain structures of the previously characterised AGOs revealed the conservation of four units: N-terminal, PAZ, MID and PIWI domains. N-terminal domain plays a critical role in target cleavage and dissociation of the cleaved RNA strand, PAZ is essential for sRNA's 2-nt 3′ overhang binding, MID anchors sRNA's 5′ phosphate, while PIWI of certain AGOs has RNase activity [33,34]. The four metal-coordinating catalytic residues aspartic acid, glutamic acid, aspartic acid, aspartic acid/histidine (D760, E788, D845 and D/H986, DED[D/H]) tetrad and a further conserved histidine residue (H798 in AtAGO1) within the PIWI domain are crucial for endonuclease/effector capacity of AGOs in vitro and in vivo [35][36][37]. The presence of the catalytic core, in some specific cases, is not sufficient for target cleavage ability [38,39]. RISCs find their target RNAs via sequence complementarity, then cleave or repress their translation during PTGS, or block their transcription due to DNA methylation and heterochromatin formation during TGS [9,40].
Environmental factors, such as extreme hot temperature, cause severe damage to the plant cells and organism by reducing fitness, endangering survival and propagation. Heat stress during cereal grain filling period causes a decrease in the synthesis of storage proteins and starch and leads to grain abortion [50]. There are a number of signal transduction pathways that initiate cellular responses in order to repair or prevent further damage [51,52]. sRNAs have been proposed to actively take part in fine-tuning the balance between development and stress responses during stress and in the post-stress period [53][54][55][56][57].
Accumulation of certain sRNAs and expression of RNA silencing factors is dynamically modulated during development and stress, suggesting that RNA silencing regulation is required for efficient response to environmental cues. Specific barley miRNAs are transcriptionally and posttranscriptionally regulated by heat [58]. A. thaliana AtAGO proteins have varying levels and patterning during different developmental stages and/or organs [59]. Brachypodium distachyon BdAGOs and BdDCLs are expressed in an organ-and developmental stage-specific manner [27,60,61]. Similarly, microarray-based profiling of rice (Oryza sativa) OsDCLs, OsAGOs and OsRDRs during vegetative and reproductive stage or in response to abiotic stresses revealed their specific and flexible regulation [25]. Maize (Zea mays) ZmAGO proteins are differentially produced under a variety of abiotic stresses [62]. These findings strongly suggest a perpetual and effective regulation of RNA silencing machinery at the transcriptional level.
DCL, AGO, and RDR gene families have several identified members in different crop species e.g., rice, maize, tomato, cucumber, grapevine, foxtail millet, and pepper [26,[63][64][65][66][67]. Although barley is among the first cultivated plants with great economic importance and with a fully sequenced genome, strikingly, there is still only scattered information available on its RNA silencing machinery and its transcriptional regulation. In the present work, we identified the members of barley DCL (HvDCL), AGO (HvAGO), and RDR (HvRDR) gene families, analysed their phylogenetic relationship to model and crop plants, investigated their domain architecture and core catalytic motifs/regions. Transcriptional changes of RNA silencing trans factors suggests that siRNA-based silencing may play a role in adaptation to heat stress conditions.

Phylogenetic Analysis, Chromosomal Localisation and Identification of Conserved Motifs
To compare the predicted HvDCL, HvAGO and HvRDR homologues with other monocotyledonous models and the dicotyledonous A. thaliana, amino acid sequences of Z. mays and O. sativa DCL, AGO, and RDR (Zm-and OsAGO, -DCL and -RDR) proteins were downloaded from UniProt database (www.uniprot.org) [68]. MEGA X software v10.1.8 [73] was employed to perform multiple alignments using the Clustal W algorithm [74] and generate neighbour-joining [75] phylogenetic trees with 1000 bootstrap replicates [76]. Barley's predicted RNA silencing components were named according to their phylogenetic relationship with the previously identified members of the same protein family. The chromosomal localisation of HvDCL, HvAGO and HvRDR genes were obtained from Ensembl Plants database [70]. SMART web server (http://smart.embl-heidelberg.de) was used to search for known protein domains. MEME web server (http://meme-suite.org/tools/meme) [77] was employed for conserved motif prediction. Search parameters were the following: optimum motif width 6 ≤ n ≤ 200, maximum number of motifs: 20. Those motifs, which did not belong to structural domains of the analysed protein family were rejected. The identified motifs were annotated using Pfam [71]. The percentage of similarity and identity of respective protein sequences were calculated with Ident and Sim online Sequence manipulation tool (http://www.bioinformatics.org/sms2/ident_sim.html) [78].

Plant Material and Abiotic Stress Treatments
Barley (H. vulgare L. cv. Golden promise) plants were grown in growth chambers at 18 °C, 16 h daylength. For expression analysis of RNA silencing-related genes, roots and leaves of two-weeksold plants were collected and frozen in liquid nitrogen. For the heat stress (HS) treatments, 16-dayold barley plants were subjected to persistent heat (40 °C for 24 h), or to mimic natural hot temperature conditions, gradient heat stress (gHS) was employed, i.e., temperature was elevated from 21 °C to 37 °C in a course of 4 h as described before [55]. Leaf samples were collected from control and treated plants immediately following the treatments, frozen in liquid nitrogen and stored at −70°C until use.

RNA Isolation, RT-PCR, Semiquantitative RT-PCR and RT-qPCR Analysis
Frozen leaf samples were homogenised in sterile mortars. RNA was isolated with Trizolate reagent (UD GenoMed Ltd., Debrecen, Hungary) according to the manufacturer's instructions. The RNA concentration and 260/280 ratio were determined with NanoDrop spectrophotometer (Thermo Fisher Scientific Inc. Waltham, MA, USA).
First-strand cDNA was synthesised from 4 μg of DNase I (Thermo Fisher Scientific Inc., Waltham, MA, USA ) treated total RNA using RevertAid kit (Thermo Fisher Scientific Inc., Waltham, MA, USA) according to the manufacturer's protocol. Every putative HvDCL, HvAGO and HvRDR genes were subjected to RT-PCR and semiquantitative PCR with manually designed genespecific primers (Table S1). HvACTIN (GeneBank ID: AY145451) was used as an internal control. For the PCRs, Phire II hot-start polymerase (Thermo Fisher Scientific Inc., Waltham, MA, USA) was used following the manufacturer's instructions. The reaction profile was the following: initial denaturation at 98 °C for 3 min, then 35 cycles of three steps consisting of denaturation at 98 °C for 5 sec, annealing of primers at a primer specific annealing temperature for 30 sec, elongation at 72 °C for 8-20 sec, and final elongation at 72 °C for 5 min.
For semiquantitative RT-PCR, the number of cycles was adjusted to 30.5 μL of the reactions were separated on 1.2% agarose gel (Lonza Inc. Rockland, ME, USA) and analysed with a ChemiDoc gel imaging system (Bio-Rad, Hercules, CA, USA) using ImageLab software (Bio-Rad, Hercules, CA, USA).
For RT-qPCR assays, 1 μg of DNase-treated total RNA and random primer was used for the first-strand complementary DNA reaction according to the manufacturer's instructions (New England Biolabs, Ipswich, MA, USA). RT-qPCRs were done using qPCR Master Mix (NEB, M3003S, www.neb.com) according to the manufacturer's instructions. RT-qPCR reactions were run in a LightCycler ® 96 Real-Time PCR machine (Roche, Basel, Switzerland). The reaction profile was the following: preincubation: 1 min 95 °C; amplification: 45 cycles of 15 sec 95 °C and 30 sec 60 °C; melting: 10 sec 95 °C, 60 sec 65 °C and continuous heating to 97 °C; cooling: 30 sec 37 °C. At least three biological and three technical replicates were assessed in each experiment and standard error bars are shown. p-values were calculated using unpaired two-tailed Student t-test to assess the significance of differences between the means of the treated and untreated samples. For a list of primers please see Table S2.

RNA-Seq Analysis
Sequence data associated with the study by Pacak et al., [79] were downloaded from the NCBI Gene Expression Omnibus database (accession number: GSE82134). The sequences were aligned to the H. vulgare IBSC v2 transcriptome and the normalised expression values (transcript per million, TPM) were calculated with kallisto v0.44.0, normalised between samples with sleuth v0.29.0 [80] and represented on heat maps that were prepared with the R package pheatmap [81]. Differential expression analysis was performed with sleuth v0.29.0 [80]. To test whether there is a significant difference between the means of the heat-stressed (HS) and the not-treated (NT) samples, a Waldtest was applied. The calculated p-values were corrected for multiple testing using the Benjamini-Hochberg method [82]. A gene was considered significantly differentially expressed if the Q-value was lower than 0.05 (5% false discovery rate).

Identification and Structural Analysis of DCL, AGO and RDR Genes in Barley
To identify barley's RNA silencing associated genes, the A. thaliana DCL, AGO and RDR protein's amino acid sequences were used as queries to build a Hidden Markov-model (HMM). According to HMM profile analysis, a total number of five HvDCLs, eleven HvAGOs and seven HvRDR genes were predicted in the barley genome. H. vulgare silencing components were named based on the closest relative in A. thaliana, or monocot species O. sativa and Z. mays ( Figure 1). The evolutionary history was inferred using the Neighbor-Joining method [75]. The bootstrap consensus tree inferred from 1000 replicates [76] is taken to represent the evolutionary history of the taxa analysed. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1000 replicates) are shown next to the branches. The evolutionary distances were computed using the Poisson correction method and are in the units of the number of amino acid substitutions per site. All ambiguous positions were removed for each sequence pair (pairwise deletion option). There were a total of 2145 (A), 1302 (B), and 1393 (C) positions in the final dataset. Evolutionary analyses were conducted in MEGA X [73].
The putative RNA silencing gene chromosomal location, structure and encoded protein predicted features are listed in Table 1. Exon number of HvDCLs was approximately the same as their A. thaliana homologues (Tables 1 and S3). Genes encoding DCL homologues in barley possess ORFs ranging from 4428 to 5994 bp. HvAGO ORF lengths are shorter, ranging from 2451 to 3654 bp. HvAGO2 and HvAGO7 possess only three exons while the other HvAGOs have a significantly higher number of exons (21)(22)(23); this is similar to A. thaliana orthologues (Tables 1 and S3). The exon number of HvRDR genes showed great variety between paralogues (e.g., HvRDR6a has two exons, HvRDR3 and 4 has 19), but are similar in the corresponding orthologues (Tables 1 and S3). HvRDRs have rather uniform ORF lengths, ranged between 3348 and 3684 bp. Isoelectric point and molecular weight of putative DCL, AGO, and RDR proteins were also estimated (Table 1). To explore the evolution of main RNA silencing gene families in barley, we analysed the genomic distribution by localising the genes on chromosomes. The exact chromosomal sites of HvDCL, HvAGO and HvRDR genes were obtained from Ensembl Plants barley (IBSC_v2) database (Table 1). The RNA silencing associated genes are distributed unevenly on the seven haploid chromosomes of barley, which is similar to RNA silencing genes in multiple monocot genomes [25,26].
The five HvDCL homologues are distributed on five chromosomes. In monocots, DCL3 gene has been duplicated: DCL3 paralogous gene duplication was suggested to occur before the common ancestor of barley and rice (ca. 60 mya) since DCL3 and DCL5 orthologues could be detected in rice, maize, and Aegilops tauschii, a progenitor of wheat [28]. HvAGO homologues are present on chromosomes 2, 3, 4, 6, 7; Besides this, HvAGO6 and HvAGO7 are located on the uncharacterised chromosome, their exact location is still uncertain, according to the latest barley genome assembly. To note, HvAGO1d is localised in close proximity of chromosome 7 telomere region. HvRDR homologous genes of barley were detected on chromosomes 2, 3, and 6, respectively. Chromosome 3 hosts four putative barley RDRs (HvRDR3, HvRDR4, HvRDR6a, and HvRDR6b). Interestingly, HvRDR1a and HvRDR1b are located in close proximity to each other, only three predicted genes separate them on chromosome 6. As they share about 81% sequence identity, they can be the result of gene duplication (Supplementary Material 1). It was suggested previously that duplication of HvRDR1 and HvRDR6 possibly reflects the divergence in disease resistance between monocots and dicots [83].

Domain Analysis and Phylogenetic Relationship of Identified Barley DCL, AGO and RDR Proteins
To investigate their potential functionality, we examined the presence and order of characteristic domains within the barley silencing proteins. For this, we performed domain search with the SMART tool [84]. All important regions of DCL proteins, including DEAD, Helicase-C, Dicer dimer (DUF283), PAZ, RNase III and dsRNA-binding motif (DSRM) were found in all HvDCLs. The domain order was consistent with DCLs of A. thaliana or other monocot species ( Figure 2). Domain search with MEME tool using the amino acid sequences of HvDCL candidates predicted the same domains and unravelled small variations between putative paralogs/orthologs, however, the domains were only partially recovered (Tables S4 and S5). Similarly, all predicted HvAGO proteins contained the characteristic domains in the precise order, namely N-terminal, Argonaute linker 1 (L1), PAZ, L2, MID and PIWI ( Figure 3, Table S4), as the previously characterised AGO proteins [33,37,85]. Besides these, we mined for glycine and arginine (GR)-rich domains upstream to the AGO N-domains, based on information from A. thaliana and Chlamydomonas reinhardtii systems [86]. HvAGO1a, 1d, 5a, 18, 2 and 7 were found to contain GR-rich regions consistent with previous data [86]. GR-rich regions may facilitate interaction with ribosomes to enable AGO tethering for the effector phase of silencing [86]. The observation that species including O. sativa, Z. mays, T. aestivum or B. distachyon also encode at least one AGO protein having GR-rich region [86], suggests that targeting of specific AGOs to translating mRNA targets is a widespread phenomenon in monocots. MEME search revealed complementary results with slight differences within certain AGO homologs (Tables S4 and S6). RDR proteins are defined by the presence of the RdRP domain. Additionally, RDRα class proteins contain an RNA-recognition motif (RRM). We have found the RdRP domain in all predicted HvRDRs and the RRM domain in RDR1a, 1b, 2, 6a and 6b proteins by SMART search (Figure 4). The de novo MEME prediction recovered the RdRP but not the RRM motifs (Tables S4  and S7). To investigate evolutionary relationship between the dicotyledonous model A. thaliana and the monocotyledonous crops rice, maize and barley RNA silencing proteins, unrooted NJ phylogenetic trees were constructed for each protein family from multiple alignments of full-length amino acid sequences (Figures 2-4). Full alignment of all A. thaliana and barley DCL, AGO and RDR proteins can be found in supplementary materials ( Figures S1-S3, respectively).
According to the phylogenetic analysis, the examined DCL sequences were grouped into clades with well-supported bootstrap values (Figures 1 and 2). Besides the four clades described in dicots (DCL1-4 clades), monocots evolved a functionally different DCL, the DCL5 [22,30]. This subbranch of DCL3 family was therefore assigned in a fifth clade (DCL5 clade). The five candidate HvDCLs were divided into these clades and named as HvDCL1, HvDCL2, HvDCL3, HvDCL4, and HvDCL5 based on their high level of sequence similarities with the other members of the same taxonomic group (Figure 2, Supplementary Material 1).
According to the phylogenetic analysis, RDR proteins of A. thaliana, maize, rice and barley were divided into four clades (Figure 4). AtRDR1, ZmRDR1 and OsRDR1 were clustered together in the first subfamily with two barley proteins designated HvRDR1a and HvRDR1b based on sequence similarity with AtRDR1 (RDR1 clade). Next clade contains one barley member, HvRDR2 (RDR2 clade). In the third clade, a total of six A. thaliana, maize and rice proteins were clustered together with two barley RDRs.
As sequence identity to each other was 42%, similarity 56%, these were named HvRDR3 and HvRDR4 (RDR3/4 clade). The last clade contained RDRs similar to AtRDR6. Two barley members of this group were designated HvRDR6a and HvRDR6b, according to their high sequence similarity to AtRDR6 (63% and 64%) and to each other (70%) (RDR6 clade) (Figure 4, Supplementary Material 1).

Functionally Conserved Amino Acid Residues in Domains of the Identified RNA Silencing Factors
To further corroborate the potential functionality and get an insight into the possible mechanistic activities of the identified barley silencing proteins, we analysed the conservation of regions involved in RNA-binding, enzyme catalysis or other essential features [31,[88][89][90].
In-depth structural data on plant DCLs are very limited and comprises only domains outside the processing centre [31,89]. The catalytic core of AtDCL4, however, was computationally modelled, providing information about the amino acids involved in dsRNA recognition, binding, or cleavage [90]. Based on these data, we thoroughly characterised barley DCLs (Figures 5 and S4). Both the surrounding region and the catalytic residues directly involved in dsRNA cleavage of RNase IIIA (E1122, D1126, N1159, K1233, D1237 and E1240) and IIIB (E1330, D1334, N1367, K1418, D1422 and E1425) were found to be well conserved between all A. thaliana and barley DCLs ( Figures 5A,B) [88,90]. Connector helix core L/IPSI/L/MM(X)11LK/R was also conserved, except in HvDCL3 ( Figure S4). The 3′ RNA binding pocket of PAZ and 5′ RNA binding pocket have high variability among the studied structures; accordingly, we could not locate these unambiguously. Nterminal part of PAZ-loop (NLL motif) responsible for dsRNA binding, was conserved, except in HvDCL2 ( Figure S4). The AtDCL4 RNase IIIA TEKCHER motif and RNase IIIB analogous HPSYN loop were predicted to interact with dsRNA; only limited conservation of these could be observed ( Figures 5A,B). T1150, R1151 (RNase IIIA) and T1358, R1415 (RNase IIIB) residues were predicted to contribute to dsRNA positioning at cleavage sites; at T1150, a conserved threonine or serine was found, implying that a hydroxylic side chain is required at this point ( Figure 5A, open arrow); T1358 residue (RNase IIIB) was also conserved ( Figure 5B, open arrow). In contrast, residues corresponding to R1151 (RNase IIIA) and R1415 (RNase IIIB) were highly variable, implying that these are not critical for dsRNA positioning. In summary, data suggest that, although binding/positioning of dsRNA is potentially different, the cleavage catalysis driven by HvDCLs mechanistically is very likely similar to the RNase III type enzymes described before [90]. Besides the regions discussed above, we have found other highly conserved residues and stretches ( Figure  S1), implying their structural and/or functional relevance in a clade-independent manner. The biological relevance of these remains to be validated experimentally. Next, we took a closer look into the conservation of functional motifs in HvAGOs. PAZ domain recognise and bind the 3′ end of sRNAs [91,92]. As data on plants are not available, we searched for conserved motifs based on D. melanogaster and human AGO data [92,93]. Only limited conservation was found in the amino acids involved in the sRNA 3′ end binding between DmAGO1 aromatic cluster α3 (L, Y, F and Y, red squares), human HsAGO-eIF2C1 (H, Y, F, Y and L, blue circles) and A. thaliana or barley AGOs ( Figure S5). These residues potentially may take part in sRNA 3′-end binding. Besides these, other highly conserved residues/motifs with potential functional importance are present within this area ( Figures S2 and S5).
We also investigated if HvAGOs retained the DED[D/H] (DEDH or DEDD) tetrad and the additional histidine (AtAGO1-H798 analogue) residue within PIWI, directly involved in enzyme catalysis [35][36][37] (Figure 6B). The DED[D/H] tetrad was fully conserved in all HvAGO proteins. Five out of eleven HvAGOs (HvAGO1a, 1d, 5a, 7, and 18) also possessed the additional histidine analogous to AtAGO1-H798. HvAGO2 contained the DEDD tetrad and the additional H, signature typical to AGO2/3 proteins [104]. However, HvAGO4a, 4b, 5b, 6, and 10 had amino acid substitutions in the AtAGO1-H798 analogue position. The highly conserved glutaminephenylalanine-valine (QF-V) motif essential in sRNA duplex recognition and sorting [105] was fully conserved in all HvAGOs ( Figure 6B). The clade-specific surroundings of D760, H798, E801 residues and QF-V motif hints to roles in the functional diversification of AGO clades. We also analysed the AGO-hook motifs involved in GW protein partner and sRNA 5′-end binding [106][107][108]: several essential residues of HsAGO2 and hydrophobic stretches within pockets 1 and 2 show conservation between A. thaliana and barley AGOs ( Figure S6) pointing to biological relevance. Interestingly, some of these regions possess clade-specific amino acid swaps. Residues directly involved in GW protein AGO hooking and sRNA 5′-end binding need to be defined experimentally.
Finally, we checked the presence of the catalytic core and its surrounding region within barley RDRs . HvRDR1a, 1b, 2, 6a and 6b, belonging to RDRα clade all possessed the canonical DLDGD, while RDRγ clade proteins HvRDR3 and 4 contained the atypical DFDGD catalytic core (Figure 7). The enclosing of the catalytic centre also carried differences between RDRα and RDRγ groups: the PHX2EC/AS upstream stretch or the WD dipeptide downstream to the catalytic core is RDRαspecific. These observations further underpin the motif and phylogenetic analysis on HvRDRs (Figure 4). Retention of the active core within all RDRα members suggests that these factors are potentially required at some stage during barley lifecycle. The high similarity between barley and corresponding genes/proteins of other species, the presence of characteristic domains in the rigorous order (Figures 2-4) and functionally important residues/motifs ( Figures 5-7) suggest that the predicted barley genes/proteins are potentially the orthologues of the previously characterised DCLs, AGOs and RDRs of A. thaliana or rice and that these genes are probably functional and required during development and environmental adaptation.

Autoregulation of RNA Silencing
RNA silencing factors including AGO and DCL transcripts were shown to be autoregulated in many plant species by a negative feedback loop involving sRNAs [109,110]. Ath-miR168 and the secondary siRNAs arising from miR168-guided cleavage negatively regulate AtAGO1 in a complex but also robust manner [39]. This autoregulatory loop seems to be widely conserved among vascular plants, and potentially also occurs in monocots. MiR168 was found in rice, wheat, maize, B. distachyon and barley as well [111][112][113][114]. For these reasons, we analysed the potential of AGO1 being negatively regulated through miR168 in barley. MiR168 target region was found in both HvAGO1a and HvAGO1d (Figure 8). HvAGO1a and HvAGO1d target-miR168 guide pairing topology and locations of mismatches are very similar to the A. thaliana system, suggesting that both HvAGO1a and 1b are probably subjects of an autoregulatory loop. This hypothesis is also supported by a degradome data of barley HvAGO1a miR168-mediated cleavage [115]. MiR403 is the negative regulator of AGO2 and AGO3. MiR403 is present only in specific dicot groups, having very similar sequence and few isoforms suggesting a recent evolutionary origin. MiR403 is probably absent in monocots [113,116]. In accordance with this, we were unable to find any potential target region of miR403 in HvAGO2. Ath-miR162 has been shown to negatively regulate AtDCL1 in A. thaliana [117]. Although found in rice, miR162 was not detected in barley, wheat or B. distachyon [112]. In line with this, we could not find or predict miR162 target region within HvDCL1 gene.

Heat Stress Significantly Alters the Expression of RNA Silencing Genes
High-temperature stress (heat stress, HS) is considered to be one of the major abiotic stresses affecting both composition of natural habitats and distribution and productivity of agriculturally important plants worldwide [118][119][120]. Plants alter their developmental pathways to re-allocate resources and ensure versatile stress management. To preliminary assess the transcriptional alterations of RNA silencing components in response to HS, RNA-seq data of barley cv. Rolap shoots provided by Pacak et al. [79] was analysed. A heat-map was generated using the normalised expression values of every putative HvDCL, HvAGO and HvRDR gene ( Figure 9A). According to the RNA-seq data, HvDCL3, HvDCL5, HvAGO2, HvAGO6, and HvRDR2 were significantly induced, while HvAGO1a was slightly downregulated in heat-stressed shoots (Figure 9). To verify the results of the in silico analysis, expression of selected silencing-related genes, including the ones that changed significantly upon heat stress in the RNA-seq experiment by Pacak et al. [79], were validated by RT-qPCR ( Figure S7). Additionally, semi-quantitative PCR consistently measured elevated HvDCL3, HvAGO6, HvRDR2 and HvRDR6a but unchanged levels of HsDCL1 and HvAGO1 following persistent, long term heat stress (40 °C/24 h) treatments ( Figure S8). Finally, to mimic natural conditions, we performed the experiment by employing gradient heat shock treatment (gHS, 21 °C to 37 °C in the course of 4 h, see Materials and Methods). Significant accumulation of HvDCL3, HvAGO2, HvAGO6, HvRDR2 and HvRDR6a transcripts was again certified. The two main factors of miRNA pathway HvDCL1 and HvAGO1a however, were not altered by this treatment either ( Figure 10). In summary, moderate direct HS (35.5 °C/48 h), prolonged heat (40 °C/24 h) and gradient heat (21-37 °C/4 h) all lead to transcriptional accumulation of siRNA-but not miRNA-based pathway components (Table S8).

Discussion
Barley is a multi-purpose crop plant: besides its economic importance, is a model species for cereal research (almost 20,000 papers were published on barley based on a PubMed keyword search). As a close relative to wheat (separated ca. ten million years), it may be a good starting point for studies difficult to conduct in wheat. Fertile hybrids may also be produced between barley and wheat, enlarging the genetic utensils available. Importantly, barley is a more resilient plant and adapts better to regions/habitats where wheat cannot be cultivated. Due to its diploid genome, it is suitable for easy obtaining of mutants through classical breeding or molecular (e.g., CRISPR mutagenesis) tools. Therefore, understanding cellular pathways and their specific features in barley is more convenient and could give a good model for crop science. Since RNA silencing plays a crucial role in the life of plants, there is a growing need to identify and characterise its central protein components in crop plants like barley.
DCL proteins are key factors of sRNA biogenesis. Unlike mammals, which have only one DCL enzyme, plants possess at least a set of four DCLs of monophyletic origin. During evolution, more complex plants tend to evolve more DCLs as a result of gene duplication events [28]. The close phylogenetic relationship within DCL1, 2, 3 and 4 clade orthologues suggests functional conservation. Mild differences, however, exist within the HvDCL and their orthologues: (i) The domain organisation of HvDCL1 is very similar to the AtDCL1, while in rice and maize DCL1 PAZ and Dicer-dimer domains are slightly altered. This change may have occurred after the divergence of the common ancestor of wheat and barley from the ancestor of rice and maize (ca. 60 million years ago) [28]; (ii) PAZ domain within monocot DCL2s (including barley) has a potentially modified structure, suggesting a distinct folding and RNA binding; (iii) when analysed, we found a highly conserved DCL catalytic core within all DCLs ( Figure 5). However, dsRNA binding may be possibly different between DCL orthologs/paralogs; (iv) in contrast to dicots, monocots evolved a fifth DCL, DCL5 [28]. DCL5 is required for the production of specific 24-nt-long pha-siRNAs in reproductive tissues of rice [30]. In maize, the absence of DCL5 causes temperature-sensitive male sterility [121]. The presence and domain conservation of HvDCL5 suggests a conserved and similar function in barley. Based on the roles of 24-nt siRNAs in A. thaliana, we speculate that perhaps DCL5 (and the 24-nt pha-siRNAs) may be involved in TGS during reproduction in barley. To reveal the functionally important differences in DCLs between dicots and monocots or within monocot lineages, further in silico, biochemical and genetic analyses are required.
AGO proteins are present in all eukaryotic organisms and can be identified by the combined presence of PAZ and PIWI domains [85]. The A. thaliana, rice, maize, and B. distachyon genomes encode 9 [39], 19 [25], 18 [26], and 16 [27] AGO genes, respectively. In wheat, two AGO genes have been described so far [122]. We describe here eleven AGO candidates in barley. Interestingly, the expansion of AGO family characteristic to monocots is not observed in barley. Although barley has one of the largest diploid genome (ca. 5.3 Gb), it seems that the restricted number of HvAGOs are able to fulfil all tasks required. Existence of functionally distinct AGO clades seems to widely spread in plants. From these clades, the AGO1/5/10 contains the so-called binder and slicers, while the AGO4/6/8/9 the modifiers. It was suggested that the AGO2/3/7 clade members are potentially more flexible, and could provide both binder/slicer and modifier functions, depending on the circumstances [123].
In barley, we have found five members of AGO1/5/10 clade. HvAGO1 and HvAGO5 have been duplicated. Grasses in general exhibit an expanded AGO1/5/10 clade [87]. OsAGO1 family possess the catalytic residues and have slicer activity [100]. We have also found perfect conservation of DEDH + H motifs in HvAGO1a, HvAGO1d and HvAGO5a but not in HvAGO5b, where the histidine residue (analogous to AtAGO1-H798, polar, hydrophobic) was replaced by a proline (nonpolar, hydrophilic) (DEDH + P, Figure 6). It was shown experimentally that H798P change turns AtAGO1 cleavage-deficient [35]. At 798 position a proline (P) residue is present in AtAGO6, raising the possibility that HvAGO5b also acts as a chromatin modifier. The observation that OsAGO5c/MEL1 is required for H3 modification reprogramming during male meiosis in rice supports this idea [124]. AGO10 has a particular function in A. thaliana: it sequesters miR165/166 to block its activity and enhance its decay. AtAGO10 is required for stem cell maintenance in SAM [125,126]. Whether this specialised role of AGO10 is preserved in monocots is not known. The catalytic core (DEDH + Y) of HvAGO10 was unique amongst clade members (AtAGO10 has DEDH + H). Although tyrosine (Y) is similar to histidine (H), as both having aromatic and hydrophobic side chains, the absence of the charge (H has a positive charge at physiological pH) may alter the activity of HvAGO10. Overall, we postulate that within HvAGO1/5/10 clade, HvAGO1a and 1b may act canonically and redundantly, while HvAGO5a, HvAGO5b and HvAGO10 probably perform differentiated/ specialised functions.
Grasses have evolved a further subgroup of AGO1/5/10, the AGO18 clade (Figure 3), that seems to be required for specific tasks: AGO18 confers broad virus-resistance in rice [127]; on the other hand, it was proposed that it binds 24-nt pha-siRNAs to regulate male reproductive organ development [128]. The phylogenetic analysis of HvAGO18 placed it close to Z. mays and O. sativa AGO18 ortholog proteins, implying similar functions.
AtAGO2 of the AGO2/3/7 clade binds 21-nt sRNAs and was shown to have roles during pathogen defence and double-stranded DNA repair [33,116,129,130]. In spite of the fact that AtAGO2 and AtAGO3 are very similar, surprisingly, AtAGO3 binds 24-nt and targets transposable elements through TGS [131]. In barley, only the homolog of AGO2 has been detected. HvAGO2 also contained the distinctive DEDD + H motif specific to AtAGO2/3 proteins ( Figure 6). Which AGO may substitute for AGO3 activity in barley, is an exciting question. The function of the second member of the clade AGO7 is more conserved: it regulates organ development through miR390 binding and generation of ta-siRNAs from TAS3 precursor [100,132]. Domain organisation and catalytic core in HvAGO7 are similar to AtAGO7 ( Figure 3). As multiple TAS3 loci have been found in barley [133], we assume conserved functions for HvAGO7.
The modifier AGO clade (AGO4/6/8/9) contains three functional members in A. thaliana (AtAGO8 is a pseudogene), with partial redundancy and specificity [134][135][136]. The clade members bind primarily 24-nt hc-siRNAs with 5′-A but also 21-nt ta-siRNAs to direct RdDM [17,137,138]. Barley genome also encodes three members of this clade, HvAGO4a, 4b and 6. The DEDH + P catalytic core of HvAGO4a, 4b and 6 is identical to AtAGO6 ( Figure 6B). Domain topology and phylogenetic analysis of HvAGO shows high similarity to A. thaliana, rice and maize clade members suggesting chromatin regulatory roles. The observation that ZmAGO104 (ZmAGO9) is needed for non-CG methylation of centromere region and knot-repeat DNA backs this theory [139]. We could not define in silico the functional homolog of AGO9. At/ZmAGO9 proteins are predominantly expressed in ovules and regulate cell fate [136,139]. Which of barley AGO4/6/8/9 clade members fulfils this task remains a future question.
Initially, RDR proteins were studied for their antiviral roles but later it became evident that they are also required during gene expression control and chromatin regulation. Several barley HvRDR genes were identified previously [83], however, we expanded this gene family with systematic identification of multiple members ( Figure 1 and Table 1). We noted a slight expansion of RDRα (5 members) and a reduction of RDRγ (2 members) clades. Phylogenetic analysis, domain organisation and catalytic core type reinforce their classification (Figures 1 and 4). There are only scattered data on RDRα protein functions in monocots: OsRDR1 is induced and required for antiviral defence [140]; maize RDR2 ortholog MOP1 takes part in heritable chromatin silencing [141]; SHL2, 4 and SHO1 proteins (RDR6 orthologs) are involved in ta-siRNA pathway in rice [142]. The close phylogenetic relationship and available functional data suggest conserved functions of RDR1, 2 and 6 in barley. The existence of two HvRDR1 and HvRDR6 homologues suggests an evolutionary selection for specialisation either during development or under stress responses. Indeed, HvRDR1a/b and HvRDR6a/b transcription were selectively induced in different organs, in response to different pathogens or under elevated temperature [83]. On the other hand, the presence of the two RDRγ proteins (HvRDR3, 4) informs about a need for their activity. Absence of one or multiple RDR3, 4, and 5 proteins from different plants, however, suggests non-essential roles [143].
To function in an equilibrated and flexible manner, RNA silencing needs to be strictly regulated at multiple layers including miRNA-directed feedback loops. Several mechanisms have been described, including miR168/AGO1, miR403/AGO2/3, miR162/DCL1 feedback loops [109,110,116,117]. In monocots, multiple loci of the miR168 family are present suggesting that AGO1s self-regulation occurs in monocots as well [113]. In our study, miR168 target sites were in silico detected in both HvAGO1a and 1d transcripts, implying that miR168-mediated regulation potentially occur in barley (Figure 8) [115]. AGO1 regulation by miR168 was also proposed in maize [114]. The presence of miR168 and that of target sites within HvAGO1a and 1b, nonetheless, is not the ultimate proof [144], therefore miR168-mediated self-regulation in monocots lacks the final biochemical evidence. miR162/DCL1 and miR403/AGO2/3 feedback regulation [116,117] seems to be absent from barley. Other components of RNA silencing (e.g., RDRs) may also be regulated. The monocot-specific miR444 indirectly activates OsRDR1 to boost antiviral silencing response [140]. MiR444 is also present in barley, raising the possibility of a similar regulation [115].
Heat stress is a major threat to barley crop yield and quality. Previous studies demonstrated the heat-regulated changes of specific miRNAs, ta-siRNAs, hc-siRNAs despite the rather stable global level of sRNAs [53,54,58,[145][146][147][148]. Transcription of silencing trans factor is also modulated by heat. Solanum lycopersicum AGO10a and AGO10b (SlAGO10a and SlAGO10b) were both activated by heat [65]. AtAGO1 is required for heat-stress memory [56]. Prolonged elevated temperature releases transgene-induced PTGS, which was epigenetically inherited trans-generationally [145]. These data show an intimate connection between environmental temperature, sRNA biogenesis/activity, silencing trans factor regulation and epigenetic/chromatin reprogramming.
To unravel heat-induced transcriptional regulation of silencing trans factors in barley, we studied RNA-seq data available [79] (Figure 9). Based on these, we selected silencing factors and assessed their change during three different heat stress regimes, including moderate HS (35.5 °C/48 h), prolonged heat (40 °C/24 h) and gradient elevation of heat (21-37 °C/4 h). Significant accumulation of HvDCL3, HvAGO2, HvAGO6, HvRDR2 and HvRDR6a mRNA was confirmed (Figures 9, 10, S7 and S8). The highly conserved miR390/TAS3/ARF pathway exerts its function via the siRNA-based subgroup of silencing factors. We searched for TAS3-derived tasiRNA targets in barley and identified three potential ARF target transcripts. Two of these targets were significantly down-regulated upon heat stress ( Figure S9) suggesting that the activity of siRNA-based silencing is potentially elevated at higher temperature. Contrarily, the principal trans factors of miRNA pathway, HvDCL1 and HvAGO1a seem to be much stable under the investigated circumstances. In summary, data from RNA-seq, semi-quantitative and RT-qPCR measurements all converge and points towards the transcriptional accumulation of factors enrolled primarily in siRNA-based silencing, including hc-siRNA, pha-siRNA and RDR6-dependent sRNA pathways. Importantly, our gradient heat treatment mimics natural situations, e.g., temperature changes during a summer day, therefore may be relevant in field conditions. As DCL3, AGO2, AGO6, RDR2 and RDR6 factors were all involved in TGS or double-stranded DNA break repair [130,137,149], RNA silencing could have chromatin regulatory and protective roles in barley during heat stress acclimation.

Conclusions
We have identified members of gene families having key roles in RNA silencing of barley and provided basic data on their genomic location, clade, phylogenetic relations, domain and motif organisation, and catalytic core build-up. Our data firmly suggests that these players are potentially functional and likely required at some point during barley's lifecycle. Transcriptional accumulation of siRNA pathway factors hints to a probable role in environmental adaptation. This work will be a stepping-stone to ask further fundamental and exciting questions that remained pending in barley and monocot RNA silencing field.
Supplementary Materials: The following are available online at www.mdpi.com/2218-273X/10/6/929/s1, Figure S1: Full-length protein alignment of all A. thaliana and barley DCL proteins, Figure S2: Full-length protein alignment of all A. thaliana and barley AGO proteins, Figure S3: Full-length protein alignment of all A. thaliana and barley RDR proteins, Figure S4: Partial alignment of A. thaliana and barley DCLs' PAZ domains, Figure S5: Partial alignment of A. thaliana and barley AGOs' PAZ domains, Figure S6: Partial alignment of A. thaliana and barley AGOs' PIWI domains, Figure S7: RT-qPCR analysis of selected RNA silencing transcripts in response to heat treatment (40 °C/24 h; HS) and non-treated (NT) barley leaves, Figure S8: Semi-quantitative RT-PCR measurement of barley RNA silencing factor expressional changes in response to sustained heat stress (40 °C/24 h), Figure S9: Targets of the conserved TAS3-derived siRNA in barley, Supplemental Material 1: Identity and Similarity scores of A. thaliana, O. sativa, Z. mays and H. vulgare homologue proteins, Table S1: List of primers used for semiquantitative PCR, Table S2: List of primers used for RT-qPCR, Table S3: Number of exons in barley RNAi genes and their Arabidopsis homologues, Table S4: Domain organisation of barley RNAirelated proteins according to Pfam, Table S5: Conserved motifs of HvDCLs identified by MEME, Table S6: Conservative motifs of HvAGOs identified by MEME, Table S7: Conserved motifs of HvRDRs identified by MEME.