Identification and Expression Analysis of Key Gene Families in RNA Silencing and Their Responses to Heat Stress in Wheat (Triticum aestivum)

1. Introduction

RNA interference is the post-transcriptional mechanism that controls gene regulation in eukaryotes is mediated by the activity of short RNAs. In addition to being involved in cross-kingdom mutualistic relationships and interaction networks, gene silencing agents like miRNAs and siRNAs are also responsible for the endogenous control of gene expression [1]. Artificial miRNAs used in conjunction with RNAi technology have shown to be an efficient means of controlling a number of biotic hazards to plants [1,2]. A thorough knowledge of the parts of this machinery is required since the RNA silencing process is crucial for many regulatory elements of plants. Small RNA biogenesis is directly impacted by RNA dependent RNA polymerases (RDRs) and Dicer-like proteins (DCLs), whilst Argonaute (AGO) is a key part of the RNA induced silencing complex (RISC) (Graham et al., 2010). RDRs are in charge of creating dsRNAs from an RNA template, whereas DCLs are in function of cleaving those dsRNAs into functional short RNAs with a length of 21–24 nucleotides. To power the gene-silencing apparatus, these sRNAs—whether miRNAs or any other type of siRNAs—are integrated into the RISC [4]. The sRNAs bind to certain AGO proteins and then direct the RISC to their respective target genes via complementary base pairing between the target mRNA and the guide strand of the sRNA. These types of gene control may be mediated by target mRNA cleavage or translational inhibition (Liu et al., 2017). Based on the kinds of short RNAs they connect with, the AGO proteins of plants and animals may be divided into three categories. While the second category of AGO proteins, known as PIWI proteins, is only found in animals and interacts with PIWI-interacting RNAs, the first category of AGO proteins is known to interact mostly with miRNAs and siRNAs. Worms were shown to contain a third class of AGO proteins that bind to secondary siRNAs [5]. The N’- terminal domain (Argo-N), PAZ domain, MID domain, and PIWI domain are four characteristic domains that have been postulated to exist in AGO proteins by a number of studies [6]. The PAZ domain has a nucleotide-binding pocket that serves as an anchor for the two 3′ overhangs of short RNAs produced by DCLs' RNase III-like activity [7]. The PIWI domain has strong functional similarities to RNase H and is known to give the AGO proteins "slicer" action [8]. Small RNAs are known to be anchored to AGO proteins by the MID domain, which also binds their 5′ phosphates [9]. By disrupting the duplex structure, the Argo-N domain may make it easier to separate the tiny RNA: target duplex after slicing [10]. Along with these domains, there may also be two linker domains called Argo-L1 and Argo-L2 between the ArgoN-PAZ and PAZ-Piwi lobes, respectively. Different species of plants have varying amounts of AGOs in their genomes. For instance, 10 AGOs have been identified for Arabidopsis, 13 for Citrus, and 19 for Maize, and 21 for Rice, whereas 21 AGO genes have been found for Saccharum [11,12,13]. Arabidopsis has ten AGO proteins, several of which have been extensively studied. AGO1 interacts with the majority of miRNAs and certain siRNAs, such as trans-acting siRNAs, which regulate plant growth and stress responses via target mRNA cleavage and/or translation suppression. The DEAD-box helicase, Helicase C-terminal domain, a Dicer dimerization domain, PAZ, Ribonuclease-III, and dsRNA binding motif are among the six conserved domains discovered in the DCLs. However, even a functioning DCL protein may lack one or more of the aforementioned domains [14]. The sequence of RDRs only contains one distinct conserved domain, the RNA-dependent RNA polymerase (RdRP) [13].

Bread wheat (T. aestivum L.) is an important cereal crop, contributing 20% ​​of the daily caloric intake and is an important source of protein and carbohydrate in the human diet. Grain yield is the most important economic parameter and is influenced by a number of biotic and abiotic factors. The demand for wheat is forecast to increase by 60% by 2050 [15]. To cope with abiotic challenges such as salinity, desiccation and heat, plants employ a variety of molecular defense mechanisms. Globally, these environmental constraints lead to average production losses of more than 50% in main crops [16]. To meet the food demands of the growing human population, breeding stress-tolerant varieties resistant to various environmental challenges, as well as improved quality and yield [17] are future goals.

Studies have revealed that the sRNAs involved in the RNAi pathway are crucial in controlling gene expression in plants under biotic or abiotic stress. In this study a set of AGO, DCL, and RDR genes were characterised in wheat. We analysed the evolutionary relationship of RNA silencing families with those of Arabidopsis, rice and analysed their chromosomal locations, gene structure, duplication events, cis-element and expression profiling. The expression patterns obtained from in-silico analysis were then confirmed using qRT-PCR. The findings presented here provide new insights into the molecular mechanism of RNAi and will aid in the understanding of the mechanisms underlying stress resistance in wheat.

2. Materials and Methods

2.1. Identification of Putative Wheat AGO, DCL, and RDR Genes

To identify the RNA silencing genes in wheat, the Phytozome database and the Rice Annotation Project (RAP) (https://rapdb.dna.afrc.go.jp/) were used to identify the AGO, DCL, RDR protein sequences of related plant species, viz A. thaliana and O. sativa,. The bread wheat proteome sequences (fp://fp.ensemblgenomes.org/pub/plants/release-51/fasta/triticum aestivum/pep/) were used as the database for BLASTP, while sequences from other species, such as A. thaliana, and O. sativa, were used as query sequences. These techniques were used to identify potential RNA silencing candidates. After eliminating duplicate results, the remaining sequences were scanned using HMMscan (https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan), the SMART database (http://smart.embl-heidelberg. de/), NCBI CDD and PFAM [18], among others [19]. AGO proteins were identified by sequences with N-terminal (pfam16486), PAZ (pfam02170), and PIWI (pfam02171) domains. However, Linker and Mid domains might or might not be found in every gene found. The presence of all the conserved domains of DCL proteins, including DEAD (pfam00270), Helicase-C (pfam00271), Dicer-dimer (pfam03368), PAZ (pfam02170), RNaseIII (pfam00636), and DSRM, was determined by batch CD search against the pfam and SMART databases using wheat peptides that showed the presence of RNase III domains (pfam00035). The peptides that showed the presence of the RdRP domain were taken into consideration as potential wheat RDRs for identifying RDRs. Bio sequencing analysis utilizing Hidden Markov Models in the HMMER database [20] also validated the locations and structural integrity of the discovered domains. However, the naming of identified genes was according to their positions in phylogenetic trees, including designated AGOs, DCLs, and RDRs of A. thaliana [21]. The physio-chemical properties were calculated by using the ProtParam tool in EaPASy web server including length of amino acids, isoelectric point (pI), aliphatic index, instability index, GRAVY value and other properties of the TaAGO, TaDCL, and TaRDR proteins.

2.3. Phylogenetic Relationship Analysis and Classification of RNA Silencing Genes in Wheat

The ClustalX 2.1 tool was used to construct multiple sequence alignments for the putative TaAGO, TaDCL, and TaRDR proteins using default parameters, and GeneDoc software v1.0 (https://genedoc.sofware.informer.com) was used to view the results (Larkin et al., 2007). The evolutionary and phylogenetic analyses were performed using MEGAX software (https://www.megasofware.net) [22]. By employing the Neighbour Joining technique on a matrix of pairwise distances calculated using the Jones-Taylor-Tshorton (JTT) matrix-based model, preliminary trees for random search were created (Kakar et al., 2018). Using bootstraps of 1000 repeats and the Maximum Likelihood technique based on the Jones-Taylor-Torton (JTT) model, were used for final phylogenetic trees construction.

2.4. Prediction of Motifs and Gene Structure

The MEME web server was used to identify the conserved motifs in the identified proteins, keeping the range of ideal motif width between 6 and 200 and the maximum number of distinct motifs at 10. (Bailey et al., 2015). Pfam and hmmscan were used to annotate the found motifs (Potter et al., 2018; Cui et al., 2020). GSDS 2.0 was used to obtain a depiction of the exon-intron organization of the AGO, DCL, and RDR family genes based on alignments of their coding sequences with the relevant genomic sequences (Hu et al., 2015).

2.5. Chromosomal Locations, Interaction Network and Orthologous Events Analysis

The physical locations of RNA silencing genes were determined from the ensembl database (https://plants.ensembl.org/Triticum_aestivum/Info/Index). The three gene families’ chromosomal positions were mapped using Mapchart v2.3 [23].

For the purpose of identifying protein-protein functional interactions, the STRING v1054 databases were employed [24]. RNA silencing protein sequences were entered into the STRING57 programme, and H. vulgare was used as the reference species to search in the database. The best-hit gene for each gene was chosen using Cytoscape58 to create a PPI network [25]. Using Cytoscape's cyto Hubba plugin, the top hub gene from the interaction network were finally determined. Orthologous genes of the TaAGO, TaDCL, and TaRDR were identified in A. thaliana by using BLASTN tool of NCBI with parameters like e value of 1^e-10 [26]. However, the paralogous genes were shown by circos using the TBtools [27].

2.6. Identification of Cis-Acting Regulatory Elements

The upstream region from the promoter sequences of 2 kb in size of each candidate gene in the TaAGO, TaDCL, and TaRDR families was used to search cis-elements by the PlantCARE program (https://bioinformatics.psb.ugent.be/webtools/plant care/html/).

2.7. Prediction of miRNA Targeting Genes

Prediction of miRNA target sites within the RNA silencing transcripts, cDNA sequences of identified transcripts were used as target gene as input to the psRNATarget server and analysed against the conserved plant miRNAs using the default parameters.

2.8. RNA-seq Derived Gene Expression Profiling

To determine in silico gene expression in different tissues, the transcripts per million (TPM) value for each TaAGO, TaDCL, and TaRDR genes from the expVIP database (http://www.wheat-expression.com/) is required. Heatmaps were created with Clustvis (https://biit.cs.ut.ee/clustvis/) [28].

2.9. Plant Material and Heat Stress Treatment

In the current study, seeds of the two contrasting wheat genotype (RAJ3765 and HS240) were obtained from the Germplasm Unit of the ICAR-Indian Institute of Wheat and Barley Research in Karnal, India. Under controlled conditions, the seeds were germinated in Petri dishes at 22°C after being sterilized with 1% sodium hypochlorite for 10 minutes and washed three times with distilled water. After five days of germination, seedlings were transferred to full strength Hoagland's solution phyta-jars and incubated for fourteen days in a BOD incubator with two sets of three biological replicates of each genotype. Two distinct wheat genotypes RAJ3765 (heat resistant) and HS240 (heat sensitive) were selected for heat stress. These plants were kept at basal temperature of 37°C for three hours, then kept at room temperature for three hours, and finally 42°C for three hours (at Acquired). Leaf samples from the baseline and acquired stress levels were taken at the time interval mentioned above. For total RNA isolation, all acquired samples were immediately wrapped in foil and frozen in liquid nitrogen at -80°C.

2.10. Expression Analysis of AGO, DCL, and RDR under Heat Stress

TRIzol Reagent was used to isolate the RNA in accordance with the manufacturer's instructions. The extracted RNA was subjected to DNase I treatment to get DNA-free RNA (NEB, USA). Using Superscript-III reverse transcriptase, 1gm of total RNA was converted into the first strand of cDNA (Invitrogen, USA). For real-time qRT-PCR analysis, the cDNA was diluted to a ratio of 1:2, and 1µl of the diluted cDNA was used as a template in a 10 µl reaction volume in accordance with the manufacturer's instructions. SYBR Green-based real-time quantitative RT-PCR analysis was performed using the BIO-RAD CFX96 (Bio-Rad). The expression data was standardised using the endogenous control wheat actin [29]. The expression was calculated as a relative fold change using the 2^△△-Ct technique [30]. The error bars display the standard deviation of the three biological replicate’s expression.

3. Results

3.1. Identification and In Silico Analysis of AGO, RDR, and DCL Genes in Wheat

To find the RNA silencing associated genes in wheat. Arabidopsis and rice protein sequencse of AGO, DCL, and RDR were used as queries species to get the HMM profiling (Supplementary file S1). After HMM profiles (Hidden Markov Model) of conserved domain, and found the respective genes from the wheat database (ensembl database) to confirm the domain of AGO, DCL, and RDR. we identified 82 encoding AGO proteins (TaAGOs), 32 genes encoding DCL proteins (TaDCLs), and 31 encoding RDR (TaRDRs) against the pfam database (Supplementary Table S2). However in wheat we found that AGO have four distinct domains- a variable N-terminal domain, conserved PAZ, MID and PIWI at C-terminal domain. Furthermore, DCLs identified conserved domains viz., RNaseIII, PAZ, and dsRNA binding motif and RdRP (for RDRs) followed by analysis against the pfam database. After evaluating the structural integrity of conserved domain we found the 81 TaAGO, 31 TaDCL, and 31 TaRDR in wheat genome. These RNA silencing genes were named based on the phylogenetic relationship exhibited by Arabidopsis RNA silencing genes sequences as previously reported (Figure 1a–c).

3.2. Physio-Chemical Properties of TaAGO, TaDCL, and TaRDR Gene Families

The analysis of physio-chemical properties of all identified genes revealed a significant variance in this study which includes molecular weights, amino acid length, isoelectric point (IP) and GRAVY of TaAGO TaDCL , and TaRDR proteins of T. aestivum as mentioned in Supplementary table S1. The majority of TaAGO proteins, in contrast to TaDCL and TaRDR proteins, have a relatively high isoelectric point (pI > 9) (Supplementary Table S1). The GRAVY (Grand Average of Hydropathicity) values of all the identified proteins are all negative, indicating that the genes in all three families are non-polar or hydrophilic by nature. In comparison to the TaDCL and TaRDR family members, TaAGOs are often more hydrophilic. (Supplementary Table S1).

3.3. Phylogenetic Analysis of DCL, AGO, and RDR Proteins

The potential orthologus gene from the RNA silencing family in wheat was found using a phylogenetic method. Arabidopsis and rice full-length amino acid sequences were used to create an ML phylogenetic tree for the studied families (Figure 1). The AGO gene family was categorized into ten subfamilies : AGO1, AGO2, AGO4, AGO5, AGO6, AGO7, AGO9, AGO10, AGO17, and AGO18 (Figure 1a). We discovered 81 AGO genes in T. aestivum, which we named as TaAGO (Figure 1a). Fifteen TaAGO proteins belong to the AGO1 clade because their amino acid identity with OsAGO1. Eight AGO members (TaAGO10, TaAGO11, TaAGO12, TaAGO16, TaAGO17, TaAGO18, TaAGO22, and TaAGO23) were identified on the AGO2 clade had sequence similarity with the AtAGO2 proteins. Likewise TaAGO64, TaAGO68, TaAGO69, TaAGO70, TaAGO77, and TaAGO78 proteins were grouped with AtAGO4, belongs to the AGO4 clade based on amino acid sequence similarities. Twenty wheat proteins, one Arabidopsis (AtAGO5) protein, and three rice proteins lead to AGO5 clade (OsAGO4, OsAGO5 and OsAGO9) (Figure 1a). Furthermore, because of high sequence similarity with Arabidopsis and rice, the Clades AGO6, AGO7, AGO9, and AGO10 had 3, 5, 10, and 6 genes, respectively. Clade AGO17 and AGO18 clade, on the other hand, each had only four genes.

Thirty one TaDCLs were distributed in in 5 clades: investigated DCL sequences were categorized into clades with well-supported phylogenetic relationships, according to the results of phylogenetic analysis. Monocots evolved a functionally distinct DCL, the DCL5, in addition to the four clades reported in dicots (DCL1–4). This fifth, seprate clade has the domain RNase3_2 and ds RBD. Apart from that, DCL1 have (TaDCL13, TaDCL15 and TaDCL17), DCL2 (TaDCL16, TaDCL18, TaDCL19) DCL3 (TaDCL1, TaDCL2, TaDCL3, TaDCL8, TaDCL9 and TaDCL10) and DCL4 have (TaDCL4, TaDCL5 and TaDCL6) amino acid sequences, respectively, and the remainder of the sequences belong to clade 5. However, TaDCL3 and TaDCL4 are two very similar proteins with 90 percent amino acid similarity that belong to the DCL3 and DCL4 clades, respectively (Figure 1b).

Four sets of RDR genes were discovered in the phylogenetic tree studies (Group I-III). However, wheat provided 26 RDR genes. Seven wheat genes (TaRDR2, TaRDR3, MaRDR4, TaRDR5, TaRDR6, TaRDR7, TaRDR15, TaRDR18, TaRDR19, TaRDR20) and one Arabidopsis and rice gene each make up Group I (Figure 1c). Furthermore, based on sequence similarity with the A. thaliana AGO protein AtRDR2 in group II, the TaRDR21 to TaRDR27 proteins belongs to the RDR2 subfamily. TaRDR8, TaRDR9, TaRDR10, TaRDR11, TaRDR13, TaRDR28, TaRDR31 protein is associated with AtRDR3, AtRDR4, and AtRDR5 a gene in Group III, but it is also tightly clustered with AtRDR5 and OsRDR3 and OsRDR5. Group IV has two genes (TaRDR12 and TaRDR14), one Arabidopsis gene (AtRDR6), and one rice gene (OsRDR5) (Figure 1c).

3.4. Conserved Domain and Motif Analysis of DCL, AGO and RDR Proteins in Wheat

The important information to further identify a protein's regulatory functions for gene expression may be found in the prediction of motifs from insight into a protein sequence. We identified and predicted typical well-distributed motifs in TaAGO, TaDCL, and TaRDR proteins. In TaAGO genes, Motifs 1, 2, 3, 5, 7, 8, 9, and 10 are associated with Piwi domain, whereas motifs 4 represent Argo-L1, rpoC2, and RNA polymerase beta. Motif 6 contained PAZ and argonaute domain (Figure 2a). The identification of eight motifs that were matched to known domains came from the study of TaDCL proteins using conserved motif analysis. Pfam annotation indicates that motifs 1, 2, 4, and 7 correspond to portions of the RIBOc and RNaseIII domains, respectively, whereas motifs 3 and 7 correspond to the C-terminal and DEAD box domains of the helicase, respectively (Figure 2b). Motifs 6 and 10 have secondary associations, which is yet unclear. The majority of the motifs found in TaRDRs are components of the RdRP domain since it is the sole conserved domain found in plant RDRs, with the exception of motifs 7 and 10, which have no specified annotations in pfam (Figure 2c). (Supplementary Figure S1) has the logos for the appropriate motifs.

3.5. Gene Structure

According to the gene structure analysis, the predicted RNA silencing genes had well-conserved gene structure and similarity with the reference Arabidopsis genes (Figure 3). Most of the AGO genes were more than thirty in number. The length of introns also varies among different TaAGOs. Out of 82 TaAGO genes, majority of TaAGO genes exhibited 22–29 intron, except for TaAGO9, TaAGO10, TaAGO11, TaAGO12, TaAGO16, TaAGO17, TaAGO18, TaAGO21, TaAGO22, TaAGO23, TaAGO26 and TaAGO78 that has only 2-7 introns (Figure 3a). The TaAGOs displayed the most introns (2–30), which were quite comparable to the AtAGOs' gene structure. Except TaDCL4, the exon-intron numbers of displayed TaDCLs, as compared to AtDCLs, showed variations (Figure 3b). The TaDCLs intron numbers (1-28) demonstrated similarity with AtDCLs. On the other hand, four genes TaRDR01, TaRDR16, TaRDR17, TaRDR29 were intronless. Whereas, three genes have single introns and other displayed 5–10 numbers of the intron in their gene structure (Figure 3c). Whereas, TaDCL, TaAGO, and TaRDR gene structures share a higher degree of similarity with their Arabidopsis orthologs, pointing to their closely related functional participation in the RNAi pathway.

3.6. Chromosomal Location of Wheat DCL, RDR and AGO Genes

Next, we analyzed the chromosomal distribution of genes encoding the AGO, DCL, and RDR in wheat using the wheat database, to get the chromosomal location (Figure 4). The corresponding genes from TaAGO, TaDCL, and, TaRDR, proteins were unevenly located in all chromosomes. However, it was observed that gene pairs e.g., TaAGO6/7 located on chromosome 2A, TaAGO42/44, TaAGO43/46 on chromosome 5B, TaAGO69/70 on chromosome 7B, and TaAGO80/82 present on chromosome 7D, exhibited the presence of these genes in close vicinity with each other whereas the rest of AGO genes were distributed on to separate chromosomes. Similarly, the presence of homologous genes in the same location was also observed in the case of TaDCL genes, such as the homologous pair of TaDCL22/23 and TaDCL29/31 are located on chromosomes 7A and 7D respectively (Figure 4). Other DCL genes were unevenly distributed on the wheat chromosome except chromosome 6. It was found that no homologous gene pairs from the same chromosome. Each wheat chromosome contains at least one RDR gene except chromosomes 1 and 5. This shows that tandem gene duplication, may have occurred as a result of the evolution of these genes. Tandem duplications are frequently regarded as one of the main driving forces behind the emergence of new biological functions.

3.7. Protein-Interaction Network Analysis

In order to study the protein-protein interaction among the 144 proteins associated with RNA silencing genes identified in wheat was constructed in-silico using the String database and using the H. vulgare as a query species (Figure 5). A conspicuous network can be seen which shows intensive cross-links among various genes. PPI network analysis showed that 25 of 144 identified proteins interacted with each other, including TaAGO9, TaAGO12, TaAGO68, TaDCL1, TaDCL4, TaDCL16, TaRDR12, TaRDR16, and TaRDR26 were considered hub-genes that interacted with multiple genes belonging to different families (Figure 5). Some of these proteins e.g., TaAGO2, TaAGO3, and TaAGO5 strongly interact only with TaDCL4. Whereas, TaAGO8 and TaAGO10 interact with TaDCL16 only. However, the other genes TaAGO24, TaAGO38, TaAGO55, TaDCL7, TaDCL8, and TaRDR13. These results indicated that various combinations of three core components of TaAGOs, TaDCLs, and TaRDRs may participate in different RNA silencing pathways in wheat.

3.8. Identification of Orthologous and Paralogous RNA Silencing Genes

We thoroughly analyzed the orthologs of T. aestivum and H. vulgare to evaluate the evolutionary connection between wheat and other plants in terms of proteins involved in small RNA machinery. There were discovered to be 62, 27, and 31 orthologous protein sequences for TaAGOs, TaDCLs, and TaRDRs, respectively. Contrarily, there were 70, 26, and 23 orthologous genes of TaAGO and TaDCL, and TaRDR of wheat and rice respectively (Supplementary Table S3). Additionally, it was found that the orthologous genes of wheat and Z. mays have 62, 21, and 20 genes of TaAGO, TaDCL, and TaRDR respectively. However, the paralogous gene pairs of TaAGO, TaDCL, and TaRDR were used to draw the circos (Figure 6).

3.9. cis-acting Regulatory Elements

The promoter regions of the identified genes (2 kb upstream of the translation start site) contained a variety of cis-acting regulatory elements. The majority of genes' upstream sequences have been shown to include the TATA box, which is one of the key regulatory elements and is located about 30 bases before the translation start point. There are also common cis-acting regulatory elements like the CAAT box and A-box in the promoter regions of many genes. According to their functional characteristics, the additional cis-acting elements found in TaAGOs, TADCLs, and TaRDRs may be divided into four groups: hormone-responsive elements, stress and defense response, plant growth and development, and light-responsive elements. The number of these elements found in each gene's promoter region is displayed in (Supplementary Table S3).

3.10. MicroRNA Targeting AGO, DCL and RDR Genes

miRNA control the various important biological processes e.g. signal transduction, growth and stress reponse. The RNA silencing genes are themselves participated in the biogenesis of miRNA and identifying its target sites of these gene families may help to understand the many processes like self-regulatory or its feedback mechanism involved in plant miRNA biogenesis. For the target analysis, using the following criteria e.g. e-value and threshold 2.0 and the conserved plant miRNA targeting the AGO, DCL, and RDR were collected from miRBase(Figure 7a). Whereas, eight and ten effective miRNAs targeted seventeen DCL genes (Figure 7b) and forty-eight RDR genes sites (Figure 7c), respectively. The conserved miRNAs are placed on the 3′ strands of the stem-loop hairpin precursors. The UPE (Unpaired energy) value of AGO genes varied from 4.2 (miR390a-5p) to 23.2 (miR172), and the DCL & RDR gene’s UPE value ranges from 14.3 (miR396-5p) to 24.1 (miR167g) and 13.5 (miR172a) to 24.5 (miR1120a) respectively (Supplementary Table S4). A lower value of the UPE indicates a higher likelihood of interaction between the miRNA and the target mRNA since it shows the relative energy needed to open the secondary structure of the miRNA around the target mRNA.

3.11. Expression Analysis of RNA Silencing Genes in T. aestivum

To analyze expression patterns of genes encoding RNA silencing in T. aestivum, transcriptional expression analysis was done in different tissues from seedlings as well as from mature plants using the RNA-seq. database (Figure 8). Overall, it was found that no clear temporal expression pattern was observed among the 144 identified genes, but expression patterns that tended to be specific to some tissues were observed. The majority of genes encoding for TaAGOs exhibited no or very weak expression levels in either leaf or shoot or other tissues. The TaAGO genes (TaAGO6, TaAGO12, TaAGO18, TaAGO25, TaAGO28, TaAGO31, TaAGO61, and TaAGO74) have higher expression in normal conditions and moderate expression under stressed conditions (Figure 8a). However, the majority of genes expressed very high in spike under normal conditions. For the TaDCL genes, seven genes (TaDCL10, TaDCL13, TaDCL15, TaDCL17, TaDCL23, TaDCL26, and TaDCL27) had high expression levels in leaf, shoot, root, grain, and spike under normal as well as stressed conditions (Figure 8b). Whereas, TaDCL6, TaDCL8, and TaDCL9 showed higher expression levels in spike tissue compared to other tissues. Out of thirty-one TaRDRs genes, it was found that only TaRDR1 is expressed higher in different tissues under normal as well as stressed conditions (Figure 8c). Whereas, TaRDR25 is expressed in leaf, shoot, root, and spike under normal conditions as compared to stressed conditions. Eleven TaRDR genes were expressed only in spike. However, TaRDR12 is expressed in roots and spike tissues relative to other tissues (Figure 8c).

3.12. qRT-PCR Expression Analysis of TaAGO, TaDCL, and TaRDR Genes under Heat Stress

The RNA silencing genes expression was determined in heat-tolerant (Raj3765) and heat-susceptible (HS277) varieties under normal growth and stressed conditions (Figure 9). However, there are thirteen genes showing variable expression in the leaf tissue under heat stress. The expression pattern of the following genes, TaAGO12, TaAGO18, TaAGO25, TaAGO28, TaAGO31, TaDCL13, TaDCL17, and TaDCL27 was found to be higher in both genotypes at 42 °C (Figure 9). However, there are some genes viz TaAGO74 (7.03 folds) which belongs to the AGO1 subfamily TaDCL23 (~5.01), and TaRDR1 (2.88 folds) that expressed higher at 42°C in Raj3765. Apart from this, the expression was downregulated in both genotypes at 37°C except TaDCL17 (~2.3 folds) (Figure 9).

4. Discussion

RNA interference is a flexible mechanism that employs machinery to silence particular genes based on their sequence in order to control the amount of gene transcript accumulation. However, investigating the role of genes of this silencing machinery's varied expression patterns under various conditions becomes critical. The AGO, DCL, and RDR genes work together to process several forms of short RNAs, which indirectly includes them in the control of a range of biological activities [31]. These three gene families' members play key roles in the synthesis of sRNAs and the efficient silencing of their targets. In previous report analysis [32,33] RNA silencing gene family characterization in wheat using computational analysis. However, in the current study we have extended the comprehensive knowledge of RNA silencing gene in terms of stress responsive experimental validation.

To date, at least 20 plant species have been used to investigate the genetic diversity of the RNA silencing gene families. RNA silencing control the process of transcription and post transcription of gene expression. AGO is one of the examples which contributed in this process. As previously reported, from three or fewer members in prokaryotic unicellular green algae like C. reinhardtii to many members in flowering plants like A. thaliana and O. sativa, the AGO family proteins have multiplied by duplication during the development of plants [34]. The current study, involved the systematic identification and thorough analysis of 82 TaAGOs which is higher in number as compared with earlier reports in barley, maize, wheat, tea and Banana [33,35,36,37] . According to the findings, whole-genome duplication or segmental duplication may have contributed to the expansion of wheat. The domain analysis of the AGO family in wheat revealed the presence of the N-terminal domain, PAZ, Piwi, Mid, and linker-1 and linker-2 domains, which are similar to the AGO domain in barley [36]. Nevertheless, it was found that only the Argo-N, PAZ, Piwi, and Argo-L2 domains were shared by all of the identified AGOs, whereas the TaAGO7 and TaAGO4 families lacked the Mid domain. TaAGO4 was the only one lacking the linker-1 domain. According to Bologna and Voinnet (2014), TaAGO1 has a glycine-rich, AGO1-specific region that is also present in the AGO1 proteins of various other plants, such as Arabidopsis and Coffea. The AGO gene family has several members, but AGO1 tends to be the one of the most frequently used in miRNA-mediated gene silencing. Other times, however, other AGO gene homologs also help to complete the specific miRNAs' silencing mechanism.

AtAGO2 of the AGO2/3/7 clade binds 21-nt sRNAs and has been implicated in double-stranded DNA repair and pathogen defense [38,39,40,41]. Only the AGO2 homolog has been found in wheat. The characteristic DEDD + H motif unique to AtAGO2/3 proteins was also present in TaAGO2 (Figure 2). It is yet unclear which subfamily may replace AGO3 functioning in wheat, and results were consistent with barley. The second member of the AGO7 clade serves a more universal purpose by controlling organ development by binding to miR390 and producing ta-siRNAs from TAS3 precursors [42,43]. TaAGO7's domain organization and catalytic core resemble those of AtAGO7 and HvAGO7 (Figure 2a).

However, A. thaliana has three functional members of the modifier AGO clade (AGO4/6/8/9), with some redundancy and specificity, whereas, AtAGO8 is reported as a pseudogene [44,45,46]. The clade members predominantly bind 24-nt hc-siRNAs with 50 -A but may also direct RdDM with 21-nt ta-siRNAs [47,48,49]. AGO4/6/9 were present in wheat, whereas AGO8 was not. Similar to AtAGO6, TaAGO4 and TaAGO6 subfamilies are encoded in the wheat genome.

Based on phylogenetic analysis and Domain topology of TaAGO have shown high similarity to A. thaliana and rice clade members suggesting chromatin regulatory roles. As previously reported, ZmAGO104 (ZmAGO9) is needed for non-CG methylation of centromere region and knot-repeat DNA backs this theory [50]. At/ZmAGO9 proteins are predominantly expressed in ovules and regulate cell fate [45,50].

The AGO18 clade, another subset of AGO1/5/10 that has developed in grasses and appears to be necessary for certain functions, is shown in Figure 1a. AGO18 gives widespread viral resistance in rice [51]; nonetheless, it was hypothesized that it binds 24-nt pha-siRNAs to control the development of the male reproductive system [52]. TaAGO18 subfamily was aligned with Z. mays and O. sativa. AGO18 ortholog proteins in the phylogenetic study, suggesting that both had comparable activities, similar results were reported in prior reports [32,33] .

DCL proteins are play an important role in sRNA biogenesis. According to earlier studies, there are a set of four DCLs of monophyletic origin in plants. The genes encoding DCLs varied among plant species from four AtDCLs in Arabidopsis thaliana [53] to eight OsDCLs in rice [54], B. napus [53], and millet [55]. Due to gene duplication events, more complex plants often evolve more DCLs during evolution [14]. The DCL1, 2, 3, and 4 clade orthologues' tight evolutionary connection supports functional conservation. However, there are slight variations between the TaDCL and their orthologues: (1) While rice and maize DCL1 PAZ and Dicer-dimer domains alter somewhat from those of TaDCL1, TaDCL1's domain organization is extremely similar to that of AtDCL1. This change may have taken place following the split between the rice and maize ancestor and the common ancestor of wheat and barley (roughly 60 million years ago) [14]; (2) the PAZ domain in monocot DCL2s (including wheat and barley) may have a modified structure, indicating a different folding and RNA binding [32,36]; and (3) upon analysis, we discovered a highly conserved DCL catalytic core within all DCLs (Figure 2b). Although dsRNA binding may change between DCL orthologs and paralogs, (4) unlike dicots, monocots have developed a fifth DCL, DCL5 [14]. In the reproductive tissues of rice, DCL5 is necessary for the synthesis of certain 24-nt-long pha-siRNAs [56]. However, we found a set of genes distinct from DCL1 to 4 and belonging to a different lineage. We can infer that this clade may belong to DCL5. The lack of DCL5 in maize results in temperature-sensitive male sterility [57]. In barley, it was reported that the conserved domain of HvDCL5 suggests a conservation and functional similarity [36]. In A. thaliana, it was observed the roles of 24-nt siRNAs, and speculated that perhaps DCL5 (and the 24-nt pha-siRNAs) may contribute to TGS during reproduction in wheat. 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.

RDR proteins were first investigated for their antiviral functions, but it has now become clear that they are also necessary for the regulation of chromatin and gene expression. Although some wheat TaRDR genes had already been discovered [58], our systematic discovery of several more members allowed us to enlarge this gene family (Figure 1c and Supplementary Table S1). We observed a small decline in RDRγ clades and an increase of RDRα. Their categorization is supported by phylogenetic analysis, domain organization, and catalytic core type (Figure 1 and Figure 2). Information on RDR protein activities in monocots is scarce: The maize RDR2 ortholog MOP1 participates in heritable chromatin silencing [59], the rice SHL2, 4 and SHO1 proteins (RDR6 orthologs) are implicated in the ta-siRNA pathway [60], and the osRDR1 is activated and necessary for antiviral defense [61]. The close evolutionary lineages and the functional data available imply that RDR1, 2 and 6 in wheat have preserved roles. One TaRDR1 homologue provide evidence of evolutionary specialization, either during development or in response to stress. Indeed, in response to various infections or at high temperatures, TaRDR1, belong to RDR1 and TaRDR12, TaRDR14 belongs to RDR6 subfamilies transcription was specifically activated in many tissues under control and stress conditions [58]. On the other hand, the two RDR proteins (TaRDR3, TaRDR4) indicate that their activity is required by their existence. However, the absence of one or more RDR3, RDR4, and RDR5 proteins from many plants indicates non-essential activities [62].

RNA silencing needs to be controlled at various steps, including miRNA-directed feedback loops, in order to work in an equilibrium and flexible way. MiR168/AGO1, miR403/AGO2/3, and miR162/DCL1 feedback loops are a few of the processes that have been discussed [39,63,64]. Numerous miR168 family locations are found in monocots, which suggests that AGO1 self-regulation also occurs in monocots [65]. In this study, it was found that wheat eleven transcripts from AGO1 that were the subject of our analysis, miR168 target sites were found in silico, suggesting that miR168-mediated regulation may be present in wheat (Figure 8) [66]. In maize and barley miR168 control the AGO1 has also been postulated [67]. The presence of miR168 and that of target sites within TaAGO8, TaAGO14, TaAGO20, TaAGO22, TaAGO54, TaAGO56, TaAGO57, TaAGO59, TaAGO60, TaAGO61, and TaAGO74, nonetheless, is not the ultimate proof [68], therefore miR168-mediated self-regulation in monocots lacks the final biochemical evidence. Like Barley, wheat also appears to lack miR403/AGO2/3 and miR162/DCL1 feedback regulation [39,69]. RDRs are one example of a component of RNA silencing that may also be controlled. Indirect activation of OsRDR1 by the monocot-specific miR444 enhances the antiviral silencing response [61]. Barley has the same MiR444 regulation as wheat, suggesting that there may be a comparable control [66]. The activity of the miR166-AGO5 family is essential for meristem growth in plants (Iki, 2017). TaAGO9 family expression was higher in tissues undergoing active meristematic development, such as roots and seeds.

Wheat crop production and quality are seriously threatened by heat stress. Previous research revealed that, particular ta-siRNAs, miRNAs, and hc-siRNAs underwent heat-regulated alterations, despite the relatively consistent global level of sRNAs (Yu et al., 2012; Zhonga et al., 2013; Kruszka et al., 2014; Li et al., 2014; Gyula et al., 2018; Szaker et al., 2020). Heat modulates the transcription of the trans factor silencing as well. Heat activated both subfamilies AGO10a and AGO10b in S. lycopersicum [70]. For heat-stress memory, AtAGO1 is necessary [71]. Transgene-induced PTGS, which was epigenetically inherited trans-generationally, is released by sustained high temperatures [72]. 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 wheat, we studied RNA-seq data available (Figure 8). Based on these, we selected silencing factors and assessed their change during different heat stress regimes, including acquired condition (37°C/3h) and basal condition (42°C/3h). Significant accumulation of TaAGO25, TaAGO28, and TaAGO31 mRNA was confirmed (Figure 8, Figure 9). Contrarily, the principal trans factors of miRNA pathway, TaDCL1, TaAGO9, and TaRDR1 seem to be much stable under the investigated circumstances. In summary, RNA-seq analysis, RT-qPCR measurements all converge and points towards the transcriptional accumulation of factors enrolled primarily in siRNA-based silencing, including pha-siRNA, hc-siRNA, and RDR6-dependent sRNA pathways. The fact that our heat treatment simulates actual natural occurrences, such temperature variations over a summer day, makes it potentially useful in environmental conditions. RNA silencing may have chromatin regulatory and protective roles in wheat during heat stress acclimation since AGO1, AGO9, DCL2, RDR1 and RDR6 factors were previously shown to be involved in TGS or double-stranded DNA break repair (Wei et al., 2012; Wu et al., 2012; Nuthikattu et al., 2013).

5. Conclusions

In eukaryotes, the functional connection between AGOs, DCLs, and RDRs is responsible for enhancing gene regulatory mechanisms such as RNA interference and RdDM. Evaluating the possible involvement of these major gene families in a commercially important crop like wheat undoubtedly aids in the engineering of wheat crops to improve crop yield and quality. In this study, 81 TaAGO, 31 TaDCL, and 31 TaRDR genes were identified in the wheat. We aim to gain a better knowledge of these gene families through using bioinformatics, including their phylogenetic relationships, domain components, genomic location, functional annotation, and subcellular localization. Despite the fact that the newly found genes have evolved in wheat, their patterns of expression in various tissues and under stress suggest that the paralogous members of the gene family have overlapping roles. These genes may have a role in wheat adaptation to stress due to the presence of stress hormone-related promoter elements in their upstream region. The genes discovered in this study may be used as crop improvement targets to create wheat cultivars that can withstand heat stress.