Two different approaches were used for genome-wide identification of the DCL, AGO, and RDR proteins: (i) In the first approach, the previously identified sRNA biogenesis protein sequences from Arabidopsis, rice, and soybean (Liu et al., 2014) were searched against the predicted gene models of chickpea, pigeonpea (http://ceg.icrisat.org/genomesequences.html) and groundnut (http://peanutbase.org/) using blastp program at an E-value threshold of 10⁻⁵; (ii) In parallel, the Hidden Markov Model (HMM) profiles of the domains with respect to DCL, AGO and RDR were downloaded from Pfam database and scanned against the predicted gene models of the legumes under study using Hmmer v 2.1.1 (Eddy, 2011).

A non-redundant set of putative sRNA biogenesis proteins identified from the two approaches were further confirmed for the presence of domains specific to each family using SMART and Pfam search. The identified genes were designated based on their phylogenetic relationship with their orthologs in Arabidopsis and soybean. The physio-chemical properties and subcellular localization of the identified proteins were calculated using ProtParam tool (web.expasy.org/protparam/) and ProtComp (www.softberry.com), respectively.

The genomic coordinates of the above identified DCL, AGO, and RDR genes were used to map their locations on a physical map of the respective legume crops. The exon-intron structure of the genes were determined using the genome annotations which were further represented using Gene Structure Display Server (GSDS; Hu et al., 2015). The members were screened for conserved motifs using Multiple Em for Motif Elicitation (MEME) from MEME suite (Bailey et al., 2009) with a motif width of 10–200 and maximum number of motifs 20. The identified motifs were annotated using Pfam database (E-value cut-off 1.0). The 1,500 bp sequence upstream of the translation initiation codon of all DCLs, AGOs, and RDRs of all three legumes were scanned for the presence of cis-regulatory elements using PlantCARE database (Lescot et al., 2002). Only the cis-elements whose matrix score was more than five were considered. In parallel, the complete set of plant miRNAs reported in miRBase (Kozomara and Griffiths-Jones, 2014) were searched against the transcript sequences of DCL, AGO, and RDR of chickpea, pigeonpea, and groundnut using psRNATarget server (Dai and Zhao, 2011) with default parameters.

The deduced protein sequences of the identified DCL, AGO, and RDR genes of chickpea, pigeonpea, and groundnut along with their counterparts from Arabidopsis and soybean were subjected to phylogenetic analysis. Initially, the multiple sequence alignments of the protein sequences from each family were carried out using ClustalW. The alignments were manually corrected and then used to calculate a p distance matrix after pairwise deletions of gaps. Finally, the phylogenetic tree was constructed based on p distance matrix using Neighbor-Joining method implemented in MEGA6 (Tamura et al., 2013) with 5,000 bootstrap replicates.

The duplicated genes (paralogs) in each legume were identified using blastp search with an E-value threshold of 10⁻⁵ and 80% sequence identity. For the identification of orthologs, the DCL, AGO, and RDR proteins of chickpea, pigeonpea, and groundnut were searched against the complete protein set of Medicago and soybean using blastp program with an E-value threshold of 10⁻¹⁰ and sequence identity of more than 80%. The orthologous relationships were further represented using circos (Krzywinski et al., 2009). The DCL, AGO, and RDR protein sequences of chickpea, pigeonpea, and groundnut were annotated using Blast2GO (Conesa et al., 2005). Further, the protein-protein interaction studies were carried out using web-based STRING database (Szklarczyk et al., 2011) which assesses and integrates both physical and functional protein-protein associations. The protein sequences of DCL, AGO, and RDR of chickpea, pigeonpea, and groundnut were searched against Arabidopsis proteins and corresponding best hits were used for network association studies. For association studies, STRING extracts curated data from various databases such as GO, KEGG, and tries to establish possible interactions among the proteins.

In chickpea, four AB responsive genotypes, two moderately resistant (ILC 3279 and ICCV 05530) and two susceptible (C 214 and Pb 7) were used for expression analysis at 7th and 11th day post-inoculation (dpi). Seedling raising and inoculum preparation were done as described earlier (Pande et al., 2011). Ten-day-old seedlings of both resistant and susceptible genotypes, acclimatized for 24 h to 20 ± 1°C temperature and 12 h of photoperiod in control environment facility, were sprayed with the A. rabiei spore suspension of 5 × 10⁴ conidia/ml until runoff. The leaf tissues were harvested on 7th (initial stress) and 11th (severe stress) dpi.

In the case of pigeonpea, three genotypes including one resistant (ICPL 20096) and two susceptible (ICPL 332-susceptible, ICPL 8863-highly susceptible) to SMD, were grown under glasshouse conditions (25 ± 2°C temperature, 16 h photoperiod). Seedlings emerging after 10 days of sowing were stapled with SMD-infected pigeonpea leaves with at least five live mites (Nene and Reddy, 1976). The first set of leaf samples was collected at 7th dpi as initial stage of stress, followed by the second sampling at 14th dpi as severe stress.

In the case of groundnut, two resistant (GPBD 4, ICGV 13208) and one susceptible (TAG 24) genotypes for rust and LLS were used in the study. ICGV 13208 was an introgression line with resistance contributed by genomic region from GPBD 4 (Varshney et al., 2014). The three genotypes were infected with both rust (Puccinia arachidis) and LLS (Phaeoisariopsis personata) pathogens. The 35-days-old plant leaves were uniformly sprayed with rust and LLS pathogens with a spore suspension of 5 × 10⁴ spores/ml using an atomizer under maintained high humid conditions for 24 h at 25°C. The leaf samples were collected at 21st, 35th, and 50th dpi. A control set was maintained without any inoculations for all the three crops. The harvested leaf samples were stored at −80°C until RNA isolation.

Total RNA from chickpea and groundnut leaf tissues was isolated using “NucleoSpin® RNA Plant” kit (Macherey-Nagel, Germany); and for pigeonpea, Plant RNA Miniprep kit, XcelGen (XG661-01, Xcelris, India) was used according to the manufacturer's instructions. The qualitative and quantitative assessment of these total RNA samples were conducted using Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) and Nano Drop 8000 Spectrophotometer (Thermo Scientific, USA). The RNA samples with RNA integrity value of more than seven and 260/280 ratio of 1.8–2.1 were used to synthesize cDNA samples from the three biological replicates. The expression of sRNA biogenesis (DCL, AGO, and RDR) genes were studied using quantitative real time PCR (qRT-PCR). cDNA was prepared using SuperScript® III First Strand Synthesis System for RT-PCR (Invitrogen, CA, USA) according to the manufacturer's instructions. The gene specific primers were designed using PrimerQuest (Integrated DNA Technologies, http://www.idtdna.com) with default parameters for qRT-PCR (Supplementary Table 1). The qRT-PCR reactions were performed using SYBR green master-mix in 96 well-plates with three biological replicates and two technical replicates using Actin in pigeonpea, GAPDH in chickpea, and ADH3 in groundnut as the endogenous control. The PCR conditions used were as follows: 2 min at 50°C, 10 min at 95°C, and 40 cycles of 15 s at 95°C, and 1 min at 60°C. The relative transcriptional level in terms of fold-change was calculated using the 2^−ΔΔCT method and student's t-test was used to calculate significance (Livak and Schmittgen, 2001).

Using similarity and HMM searches, a total of four DCLs (DCL1, DCL2, DCL3, and DCL4) were found both in chickpea and pigeonpea similar to Arabidopsis (Xie et al., 2004). In the case of groundnut, three and four DCLs were found in A. duranensis and A. ipaensis respectively (Table 1). A higher number of DCLs in A. ipaensis can be attributed to a higher number of duplicated events in A. ipaensis in comparison to A. duranensis (Bertioli et al., 2016). Varied number of DCLs were reported in other crops viz., rice (8), maize (5), soybean (7), sorghum (3), and foxtail millet (8) (Kapoor et al., 2008; Qian et al., 2011; Liu et al., 2014; Yadav et al., 2015). The phylogenetic analysis grouped DCLs into four clades depicted as DCL I, DCL II, DCL III, and DCL IV (Figure 1A). The DCLs ranged in length from 1,410 to 1,864 amino acid (aa) in chickpea, 1,203 to 1,913 aa in pigeonpea, 1,376 to 1,949 aa in A. duranensis, and 1,398 to 1,947 aa in A. ipaensis. Interestingly, in all three legumes the shortest length was shared by clade II DCLs. All DCLs were found to be localized in nucleus supporting the findings that sRNA precursors undergo nuclear processing in plants (Papp et al., 2003).

A total of 13 AGO genes each in chickpea, pigeonpea and 22 (11 each from A. duranensis and A. ipaensis) in the case of groundnut were identified (Table 2). However, in an earlier study in Arabidopsis, a total of 10 AGOs were identified with AGO1 involved in miRNA and ta-siRNA biogenesis; AGO4 and AGO6 responsible for DNA methylation (Baumberger and Baulcombe, 2005). The length of identified AGOs varied from 776 to 1,094 aa in chickpea, 821 to 1,054 aa in pigeonpea, 854 to 1,033 aa in A. duranensis, and 792 to 1,056 aa in A. ipaensis (Table 2). The phylogenetic tree classified the AGO proteins into three clades: AGO I, AGO II and AGO IV in all the three legumes, in accordance with the findings in Arabidopsis (Xie et al., 2004), and soybean (Liu et al., 2014). Chickpea and pigeonpea shared six, three, and four members each in clades, AGO I, AGO II, and AGO IV, respectively, in contrast to groundnut where A. duranensis and A. ipaensis had four, three, and four members each in AGO I, AGO II, and AGO IV clades, respectively (Figure 1B). Phylogenetic analysis clustered AGO1, AGO5, AGO10 in clade I, AGO3, AGO7 in clade II, and AGO4, AGO6, AGO9 in clade IV in all three legumes. We did not find AGO2 and AGO8 homologs in any of the three legumes considered in this study. Further, all the AGO3 genes were found to cluster with AGO2 and AGO3 genes of Arabidopsis. Interestingly, we found the presence of two homologs of AGO3 (AGO3a and AGO3b) in all three legumes, similar to soybean (Liu et al., 2014) but in contrast to Arabidopsis (Xie et al., 2004). It can be postulated that the two homologs of AGO3 compensates for the absence of AGO2 in legumes. Absence of AGO8 protein in the legumes is in agreement with the findings in Medicago and soybean that also lacked AGO8 homologs (Capitao et al., 2011; Liu et al., 2014). This could be attributed to the loss of AGO8 during the course of speciation. Interestingly, chickpea clade IV AGO was marked by the presence of two AGO6 homologs (AGO6a and AGO6b). This second AGO6 homolog was absent in pigeonpea and peanut. Instead, these species possess one AGO9 homolog that was absent in chickpea implying the legume specific diversification of sRNA biogenesis genes. All members of clade I AGOs of chickpea, pigeonpea, and groundnut were predicted to be localized in cytoplasm supporting the findings that AGOs are known to be present in high concentrations in cytoplasm and interaction of AGO1 with miRNA occurs in the cytoplasm to facilitate the RISC to find its target mRNA. However, most members of clade II and IV were predicted to be localized in the nucleus.

Chickpea, pigeonpea, A. duranensis, and A. ipaensis encodes for five RDRs each. A similar number of RDRs (6) were identified in Arabidopsis responsible for the formation of double stranded RNA (dsRNA) molecules (Wassenegger and Krczal, 2006). Length of RDRs ranged from 987 to 1,228 aa in chickpea, 957 to 1,135 aa in pigeonpea, 795 to 1,131 aa in A. duranensis, and 839 to 1,116 aa in A. ipaensis (Table 3). Similar to previous reports (Xie et al., 2004; Kapoor et al., 2008; Qian et al., 2011; Liu et al., 2014), the phylogenetic analysis grouped chickpea and pigeonpea RDRs into four clades—RDR I, RDR II, RDR III, and RDR VI. Each clade contains one RDR gene per species except clade I that contained two homologs per species. By contrast, RDRs from groundnut were clustered into three clades each one containing two members except clade RDR II (Figure 1C). The presence of two homologs of RDR III (RDR 3a and RDR 3b) and absence of RDR VI in groundnut could be explained by incorporation or loss of these genes in early diverging Dalbergioid species which diverged from other Papillionidoids some 55–56 million years ago (Rispail and Rubiales, 2016).

DCLs are large proteins with multiple domains, namely, DEAD, Helicase-C, Dicer-dimer, PAZ, two tandem Ribonuclease III, and double stranded RNA binding domain (dsrm). The characteristic DCL domain organization was found in all chickpea, pigeonpea, and peanut DCL homologs except CcDCL2, which lacks the DEAD domain. In the case of chickpea and groundnut, all DCLs contained two copies of dsrm domain except clade II DCLs. The clade II DCLs harbored a single copy of dsrm domain in accordance with clade II DCL proteins of Arabidopsis and soybean. Interestingly, in pigeonpea only clade I DCL was found to contain two copies of dsrm domain (Supplementary Figures 1A–C). The differences between sequences are shared by all species and correspond to DCL clade specific changes which might have a functional implication. Motif 1 (CY[QE]RLE[FY][LV]GD[AS]VLD) corresponding to RNase III domain was found among all the members of studied legumes. Another motif PK[AV]LGD[IL][VI]ES was seen in all members of chickpea, pigeonpea, and groundnut except for clade IV of groundnut. In the case of groundnut, motif 8 (WLPHDAPPSA), and motif 9 (EMQDFANEPENEEQPS) were seen to be specific for clade I and II, respectively, but were missing in both chickpea and pigeonpea (Supplementary Figures 1A–C). Exon-intron structure of DCL genes identified the maximum number of introns in clade IV in chickpea, clade III in pigeonpea, clade II in A. duranensis, and clade III in A. ipaensis. The number of introns ranged from 19 to 22, 17 to 26, 24 to 33, and 23 to 29 in DCL clade I, II, III, and IV, respectively (Supplementary Figures 1A–C).

AGOs are marked by the presence of four domains: N-terminal domain, PAZ, MID, and PIWI domain. MID domain is known to have a nucleotide specificity loop which makes it a recognition and binding center for sRNAs (Frank et al., 2012). In the current study, AGOs identified in the three legumes contained ArgoL1 (Argonaute linker 1 domain), ArgoL2 (Argonaute linker 2 domain), and Argomid (mid domain of Argonaute) domain beside the PIWI, PAZ, ArgoN (N-terminal domain of Argonaute) domain. An additional, Gly-rich_AGO1 domain was also found as a signature sequence of CcAGO1, CaAGO1, AdAGO1, and AiAGO1 proteins (Supplementary Figures 2A–C) in agreement with observations in Arabidopsis and soybean (Xie et al., 2004; Liu et al., 2014). In AGO family, motif 1 (PGTVVD[TS]KI[CT]HP) was found in all members of chickpea, pigeonpea, and groundnut. A conserved stretch of T[I/L]IFGADVTHP was seen in members of AGO I clade of all three species. This stretch was substituted by T[L/M]ILGMDVSHG in clade II and V[I/M] F[I/M]GADVTHP in clade III for both chickpea and pigeonpea. Motif 7, 16, and 20 in chickpea, motif 11, 12, and 13 in pigeonpea, and motif 6, 16, 18, and 19 in groundnut were shared by clade IV members only. Motifs 18 of chickpea, 17 of pigeonpea, and 8, 15 of groundnut were shared by members of clade II. Intronic organization identified 2–22, 1–21, and 1–26 introns in AGO genes of chickpea, pigeonpea, and groundnut, respectively. Members of AGO I and AGO IV clade were identified with higher number of introns as compared to AGO II members in all three legumes. These gene structures were similar to the ones found in rice (Kapoor et al., 2008). AGOs are endonucleatically active due to the presence of a conserved metal-chelating Asp-Asp-His (DDH) motif in their PIWI domain. Along with DDH (D760, D845, H986), one more conserved residue H798 of AtAGO1 is also known to play an important role in AGO slicer activity. However, studies on Arabidopsis, rice, maize and soybean reveals that DDH/H is not a quintessential motif. We also observed that in AGO proteins the DDH/H motif was substituted by patterns like DDD/H, DDH/P, DDH/S. Multiple sequence alignment of chickpea, pigeonpea, and groundnut AGOs with Arabidopsis AGOs unveiled the presence of conserved DDH/H motif in seven chickpea, eight pigeonpea, and nine groundnut (five A. duranensis, four A. ipaensis) AGOs (Table 4). It was seen that generally the members of the same clade exhibited the same patterns (Supplementary Figures 2A–C).

All RDR genes are known to have a conserved RdRP domain. However, clade I and clade II in the legumes studied here have an additional conserved RNA Recognition Motif (RRM) as also seen in clades 1 and 2 of Arabidopsis and soybean. In RDRs, among all the motifs analyzed, motif 2 in chickpea and pigeonpea, and motif 6 in groundnut was seen to be conserved among all the RDR members. A short sub-sequence of this motif, DLDGD was seen to be conserved in all RDR clades except clade III where it got converted to DFDGD. This motif is seen conserved in other species as well, and corresponds to the catalytic β′ subunit of DNA-dependent RNA polymerases (Zong et al., 2009). Another motif, PCLHPGD[V/I]R was seen to be conserved in all RDR clades except RDR III, where it got converted to PGLHFGDIH in all three legumes. The gene structure analysis identified introns ranging from 1 to 16 in chickpea, 2 to 17 in pigeonpea, 3to 17 in A. duranensis, and 2 to 18 in A. ipaensis. Interestingly, clade III RDRs were seen to have the maximum number of introns in all three legumes (Table 3, Supplementary Figures 3A–C). The motif analysis suggests a strong conservation of DCL, AGO, and RDR sequences among the studied legumes. However, the differences among the clades of a same protein family indicate the non-redundant structure and function of members of different clades of a family.

Chickpea DCLs were distributed on pseudomolecules Ca1, Ca3, Ca5, and Ca7 whereas in pigeonpea, only CcDCL1 was found to be present on pseudomolecule CcLG03 while rest of the DCLs were scattered on unassembled scaffolds (Table 1). All AGO genes except CaAGO1 in chickpea were found to be uniformly distributed on all eight pseudomolecules. Similarly, in the case of pigeonpea, most of the AGO genes (10) were identified on pseudomolecules with a maximum number of genes on pseudomolecule CcLG02 (Table 2). In chickpea, two RDRs—CaRDR1a and CaRDR1b were located on pseudomolecule Ca4. The other RDRs, CaRDR2, CaRDR3, and CaRDR6, were positioned on pseudomolecule Ca5, Ca7, and Ca8, respectively. In pigeonpea, only CcRDR2 was placed on pseudomolecule CcLG11; all other RDRs were scattered on scaffolds. In groundnut, all the DCLs, AGOs, and RDRs were distributed on pseudomolecules (Tables 1–3, Figures 2A–D).

Gene duplication underlines the phenomenon of evolution as it functions as a template for subsequent modifications acted upon by natural selection. Gene families are nothing but true paralogs resulting from gene duplication events. In this study, we identified the duplicated events of sRNA biogenesis genes in the genomes of three SAT legumes. A paralog for an identified gene was considered only if it retained at least one of all the essential domains responsible for defining that gene to be a part of a specific family. Considering the above criteria, no duplicated DCL gene was found in chickpea, whereas two were found in pigeonpea. Likewise, one duplicated DCL gene was found each in A. duranensis and A. ipaensis. In chickpea and pigeonpea, three and four pairs of duplicated AGO genes were found, respectively. In groundnut, A. duranensis harbored one, whereas A. ipaensis shared two events of duplication of the same AGO gene. No duplication was observed for RDR genes in both chickpea and pigeonpea. By contrast, in groundnut, two pairs of duplicated RDR genes were identified for each of the two groundnut sub genomes (Supplementary Table 2).

The orthologs of DCL, AGO, and RDR proteins of chickpea, pigeonpea, and groundnut were identified in two closely related legumes, Medicago, and soybean (Figures 2A–D). Chickpea was seen to share the maximum orthologs with Medicago and pigeonpea with soybean. In the case of groundnut, DCL and AGO orthologs were present in Medicago and soybean, whereas no RDR orthologs were seen in any of the two legumes (Supplementary Table 3). These findings can be correlated with the divergence of these legume species during the course of evolution. For instance, chickpea and Medicago belong to the same clade, i.e., Galegoid, pigeonpea and soybean belong to the Phaseoloid clade while groundnut belongs to the Dalbergioid clade.

Promoter analysis identified a number of cis-acting elements in all members of DCL, AGO, and RDR of all three legumes. Light responsive elements were the most frequent followed by those related to stress responses and hormones (Supplementary Tables 4–6). Prevalence of a high number of light responsive elements prompts the role of these genes in regulation of photoperiodic control of flowering; however, this would require further validation for arriving at a conclusion. Regulatory elements corresponding to endosperm expression, and circadian control were also seen. Cis-acting elements relating to drought, heat, and low-temperature stress suggest the possible role of these genes in regulating the expression in response to these stresses. The presence of hormones, such as ABA, ethylene, and gibberellin related elements prompt the involvement of these genes in signaling pathways as well.

All the DCL genes in chickpea, pigeonpea, and groundnut were recognized as miRNA targets. The DCL1 gene irrespective of the legume species studied in this study was found to be targeted by miR162. One of the miRNAs, miR1515 involved in hyper nodulation in soybean is known to target DCL2 (Li et al., 2010). The same miRNA was found to be targeting CcDCL2, AdDCL2, and AiDCL2 in this study (Supplementary Table 7). A total of 10, 11, 8, and 8 AGO genes were identified as targets for miRNAs in chickpea, pigeonpea, A. duranensis, and A. ipaensis, respectively. It has been reported that in Arabidopsis, the co-regulation of AGO1 and miR168 is indispensable for normal miRNA functioning and plant development (Vaucheret et al., 2004). AGO1 of chickpea and groundnut (both A. duranensis and A. ipaensis) was also seen to be targeted by miR168 substantiating the indispensability of miR168 in plant growth and regulation. However, in pigeonpea this pattern was not observed (Supplementary Table 8). In case of RDRs, a total of 2, 4, 4, and 3 RDR genes were identified as the targets for miRNAs in chickpea, pigeonpea, A. duranensis, and A. ipaensis, respectively. miR156, one of the stress responsive miRNAs (Stief et al., 2014), was seen to target a wide range of genes which include AGO genes in groundnut, and RDR genes in chickpea and pigeonpea suggesting the involvement of various mechanisms in sRNA regulation (Supplementary Table 9).

DCL, AGO, and RDR are broadly known to facilitate the sRNA biogenesis. However, each of them has a precise role in the molecular and biological modus operandi of sRNA biogenesis and regulation of gene expression. Annotation of the identified respective family members revealed that clade IV DCLs were involved in a biological process of mitotic cell cycle and virus induced gene silencing (VIGS; Supplementary Table 10). Among AGOs, clade IV AGO6 was found to be involved in siRNA binding. Clade II and VI RDRs together were involved in VIGS; however, clade II RDR was identified to play a role in DNA methylation and response to fungus (Supplementary Tables 11, 12). Only clade II RDRs showed antifungal response in all three legumes (Supplementary Table 12). In-silico protein interactions among the members of DCL, AGO, and RDR showed very strong associations among themselves suggesting that all three proteins interact with each other to accomplish the task of RNA interference (Supplementary Figures 4A,C,E). Along with this, these proteins also showed strong interactions among themselves in response to stress (Supplementary Figures 4B,D,F). It was seen that most of the interacting partners, such as HYL1, NRPD, SE of these families were common in all three species. These interactions indicate the role of these proteins in processing the pri-miRNA into miRNA in the light of the fact that DCL1, SE, TOUGH, and HYL1 in Arabidopsis are well-known to form a complex to process pri-miRNA into miRNA duplex in the nucleus (Fang and Spector, 2007; Song et al., 2007).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.