Genome-Wide Identification of Key Components of RNA Silencing in Two Phaseolus vulgaris Genotypes of Contrasting Origin and Their Expression Analyses in Response to Fungal Infection

RNA silencing serves key roles in a multitude of cellular processes, including development, stress responses, metabolism, and maintenance of genome integrity. Dicer, Argonaute (AGO), double-stranded RNA binding (DRB) proteins, RNA-dependent RNA polymerase (RDR), and DNA-dependent RNA polymerases known as Pol IV and Pol V form core components to trigger RNA silencing. Common bean (Phaseolus vulgaris) is an important staple crop worldwide. In this study, we aimed to unravel the components of the RNA-guided silencing pathway in this non-model plant, taking advantage of the availability of two genome assemblies of Andean and Meso-American origin. We identified six PvDCLs, thirteen PvAGOs, 10 PvDRBs, 5 PvRDRs, in both genotypes, suggesting no recent gene amplification or deletion after the gene pool separation. In addition, we identified one PvNRPD1 and one PvNRPE1 encoding the largest subunits of Pol IV and Pol V, respectively. These genes were categorized into subgroups based on phylogenetic analyses. Comprehensive analyses of gene structure, genomic localization, and similarity among these genes were performed. Their expression patterns were investigated by means of expression models in different organs using online data and quantitative RT-PCR after pathogen infection. Several of the candidate genes were up-regulated after infection with the fungus Colletotrichum lindemuthianum.


Introduction
Small RNAs have regulatory roles in a multitude of biological processes, including stress responses, development, metabolism, and maintenance of genome integrity, in a sequence-specific manner [1]. Although heterogeneous in size, sequence, genomic distribution, biogenesis, and action, most of these small RNA molecules mediate repressive gene regulation through RNA silencing [2]. RNA silencing refers to a variety of mechanisms where a small RNA molecule interferes with a given nucleotide sequence. Plant RNA was estimated to have occurred ca. 110 000 to 165 000 years ago [29,30]. In that context, two genome assemblies of the common bean are available, one for genotype G19833 of Andean origin [30] and one for genotype BAT93 of Mesoamerican origin [31]. AGO and DCL genes have been analyzed in the Andean G19833 genotype leading to the identification of 15 PvAGO genes and 6 PvDCL genes [32]. Consequently, except for the report of de Sousa Cardoso et al. [32], our knowledge of the RNA silencing mechanism in common bean remains quite poor.
The aims of this study were to identify and characterize, by in silico analysis, the genes involved in RNA silencing, including AGO, DCL, RDR, DRB, NRPD1, NRPE1, and NRPD2/NRPE2 in common bean. Taking advantage of the availability of two genome assemblies of contrasting origins (Andean and Mesoamerican), we wanted to address the evolution of these genes on a short time scale. Their expression patterns were investigated in different organs using online data and after infection with the fungus C. lindemuthianum by quantitative RT-PCR. The identification of these core components to trigger RNA silencing in this non-model plant species of worldwide economic relevance pave the way for further investigation.
Candidate proteins were named based on their phylogenetic proximity to known members in A. thaliana, soybean, and/or M. truncatula. The prefix PvA or PvM was added for sequences originating from G19833 (Andean) or BAT93 (Meso-American), respectively.

Protein Sequence Alignment and Phylogenetic Tree Building
The complete sequence of each putative AGO, DCL, RDR, and DRB proteins were aligned using Muscle [40], and the resulting alignments were manually optimized using SeaView [41]. For a given gene, when more than one isoform was identified, the longest was selected for the alignment. Aligned sequences were then analyzed using ProtTest3 [42] Genes 2022, 13, 64 4 of 19 to estimate the best phylogenetic model. Maximum-likelihood trees were generated with PhyML [43]. Bootstrap values were computed with the consensus of 1000 trees generated with PhyML. The resulting phylogenetic trees were displayed using MEGA version 7 [44]. For phylogenetic analysis, the common bean sequences were completed with AGO sequences from soybean [17], DCL sequences from soybean, and M. truncatula [14], and RDR and DRB1 [also known as HYPONASTIC LEAVES 1 (HYL1)] sequences from soybean [15]. The location of each PvA_AGO, PvA_DCL, PvA_RDR, PvA_DRB, PvA_NRPD, PvA_NRPE gene on G19833 chromosomes was determined by tBLASTn searching against the G19833 genome. Molecular weights (Mol. Wt.) and isoelectric points (pI) were determined using the Pepstats program from EMBOSS [45] analysis package. The number of isoforms in G19833 (v1.0) and BAT93 (v10) were obtained from corresponding official annotations in the Phytozome (V9.0) and BAT93 genome data repositories, respectively. Protein similarity and identity percentage were calculated with needleglobal pairwise alignment [45]. The number of introns in the CDS was obtained from manual reannotation performed in the Artemis platform [35].
2.6. Plant Materials, Infection with C. lindemuthianum, RNA Extraction, and RT-qPCR Analysis Infections of the common bean Andean landrace JaloEEP558 with the incompatible strain 100 of C. lindemuthianum were carried out as previously described in Richard et al. [49]. A time-course gene expression analysis was conducted at 6, 24, 48, 72, and 96 hpi in JaloEEP558 seedlings infected with strain 100. For each time, one of the two cotyledonary leaves from three different inoculated plants and control plants were sampled and flashfrozen in liquid nitrogen for RNA isolation and RT-qPCR analysis.
Total RNA extraction and quantitative RT-PCR (RT-qPCR) experiments were performed as described in Richard et al. [49]. The expression analyses of the genes PvAGO1, PvAGO2a, PvDCL2a, PvDCL2b, PvAGO4a, PvAGO4b, and PvAGO4c were performed using the gene-specific primers listed in Supplementary Table S1. Gene expression was normalized with four reference genes (PvUkn1, PvUkn2, PvIDE, and PvAct11) [50] (Supplementary  Table S1). For each gene, gene expression in mock condition was used to calibrate gene expression in infected plants at each time point. Relative gene expression in inoculated leaves compared to mock leaves was calculated using the method 2 −∆∆Ct on three technical replicates and two biological replicates [51]. Statistical comparisons were carried out using unpaired t-tests between each mean value (at t = 6, 24, 48, 72, 96 hpi) and the corresponding mean value at t = 0 hpi.

Six Putative DCL Genes Are Present in P. vulgaris Genome
The search for homologous DCL sequences in the P. vulgaris genome generated six full-length DCL genes recovered from both G19833 and BAT93 genomes (Table 1). These genes were named using the prefix PvA_ or PvM_ to indicate genotype G19833 (Andean) or BAT93 (Meso-American), respectively, or PvA/M to indicate a gene present in both genotypes. PvA/M prefix was then followed by an identifier for their Arabidopsis homologs determined by phylogenetic analysis (e.g., PvA_DCL1 corresponds to the AtDCL1 gene). For paralogs, a letter (a, b, c . . . ) was used as the suffix. The same nomenclature was used for all genes involved in RNA silencing described in this study. Dicer-like 1-4 occurred as monophyletic groups containing DCLs from P. vulgaris, G. max, M. truncatula, and A. thaliana. Our manual annotation allowed us to identify PvM_DCL2c that was not present in the automatic annotation of BAT93 assembly. In P. vulgaris, for both BAT93 and G19833, DCL1, DCL3, and DCL4 occurred as single-copy genes, while DCL2 had three paralogs (PvA/M_DCL2a, PvA/M_DCL2b, PvA/M_DCL2c) ( Figure 1, Table 1). The six DCL genes in the common bean presented high levels of protein identity between BAT93 and G19833 (> 97% protein identity). PvA_DCL2a and PvA_DCL2b were separated by 2.5 kb on chromosome 6, while PvA_DCL2c was located on chromosome 8 ( Figure 2). Despite their tight physical linkage, DCL2a and DCL2b were phylogenetically separated (Figure 1), such that PvA/M_DCL2b and PvA/M_DCL2c grouped with GmDCL2b, while PvA/M_DCL2a grouped with GmDCL2a ( Figure 1). Manual inspection of flanking genes in the P. vulgaris and G. max genomes showed that both copies of DCL2 (a and b) are located in a syntenic region (Supplementary Figure S1). Indeed, in both species, the duplicated DCL2 genes were flanked by genes encoding a histidinol dehydrogenase and a protein male sterile 5 on one side and by genes encoding a stress up-regulated Nod 19 and 3-hydroxyisobutyrate dehydrogenase on the other side (Supplementary Figure S1). Amplification of DCL2 genes has also been observed in M. truncatula, which has three copies [14]; however, these DCL2s formed a separate clade ( Figure 1). The PvA/M_DCL proteins ranged in length from 1388 to 1975 amino acids (aa) ( Table 1). As observed for other legume species, the smaller DCL proteins occur within the DCL2 clade [52].

Thirteen AGO Genes in Common Bean Genome
The search for homologous AGO sequences in the P. vulgaris genome generated 13 fulllength AGO genes recovered from both G19833 and BAT93 genomes (Table 1). Our manual annotation allowed us to correct PvM_AGO2a by fusing two distinct genes from BAT93 automatic annotation leading to a 971 aa long PvM_AGO2a protein, sharing 99.3% of protein identity with the G19833 homolog ( Table 1). The length of the identified AGOs varied from 886 to 1063 aa. The Pv AGO genes were spread on 8 out of 11 common bean chromosomes, with two genes (PvA_AGO4a and PvA_AGO4b) organized in a tandem array on chromosome 8 and separated by ∼20 kb (Figure 2). The phylogenetic tree classified the AGOs proteins into three clades: AGO 1/5/10, AGO 4/6/8/9, and AGO 2/3/7 ( Figure 1). For each 13 Pv AGO genes, a clear orthology relationship was identified between G19833 (PvA_AGO) and BAT93 (PvM_AGO), testifying to the absence of recent gene duplication or deletion for this AGO gene family ( Figure 1). In particular, the gene duplication leading to PvA/M_AGO4a and PvA/M_AGO4b occurred prior to the Andean/Mesoamerican gene pool divergence.

Five RDR Genes in the Common Bean Genome
Common bean G19833 and BAT93 genomes contain five RDR genes each (Table 1), located on chromosomes 3, 4, and 9 ( Figure 2). Our manual annotation allowed us to correct PvM_RDR1a by fusing two distinct genes from BAT93 automatic annotation leading to an 1139 aa long PvM_RDR1a protein, sharing 99.3% of protein identity with its G19833 homolog ( Table 1). The length of RDRs ranged from 980 aa to 1222 aa ( Table 1). As previously observed [12,53], phylogenetic analysis grouped RDR into four clades (RDR1, RDR2, RDR3, RDR6) with clade RDR3 containing three Arabidopsis members (AtRDR3, AtRDR4, AtRDR5) out of the 6 AtRDR ( Figure 1). Concerning P. vulgaris, each clade contained a single Pv RDR gene, except for clade 1, which contained two RDR1 paralogs (PvA/M_RDR1a and PvA/M_RDR1b) closely linked on chromosome 3 and separated by 10 kb (Figures 1 and 2). Similarly, two RDR1 paralogs were also identified in chickpea and pigeonpea genomes [50], suggesting that they could correspond to an ancient gene duplication.

Ten DRB Genes in Common Bean Genome
Ten DRB genes were identified in both G19833 and BAT93 genomes (Table 1) with a clear orthology relationship, suggesting no recent duplication/deletion for this gene family in common bean (Figure 1). Our manual annotation led us to correct PvM_DRB4b by fusing two distinct genes from BAT93 automatic annotation leading to a 248 aa long PvM_DRB4b protein, sharing 99.6% of protein identity with its G19833 homolog (Table 1). Compared to Arabidopsis, the common bean genome experienced an amplification of the DRB1 gene family (five members) as well as the DRB4 gene family (two members). A clear ortholog of AtHYL1, a key interactor of DCL1 in miRNA biogenesis [54], referred to as PvA/M_HYL1, was identified on common bean chromosome 9 (Figures 1 and 2). The 10 common bean DRB genes were spread on seven chromosomes, with PvA/M_DRB1d and PvA/M_RDB1c genes tightly linked on chromosome 8.

Common Bean Pol IV and Pol V
In order to gain insight into the Pol IV and Pol V complex in the common bean genome, the largest and second-largest subunits of Pol IV and Pol V were searched by seeking AtNRPD1, AtNRPE1, and AtNRPD2/NRPE2 against common bean BAT93 and G19833 genomes with tBLASTn. Common bean encodes one NRPD1, one NRPE1, and one NRPD2/NRPE2, and hence they are named PvA/M_NRPD1, PvA/M_NRPE1, and PvA/M_NRPD2/NRPE2 (Table 1). These three proteins present a high level of identity (> 99%) between BAT93 and G19833. They are located on chromosome 2 (PvA_NRPD1), 11 (PvA_NRPE1), and 9 (PvA_NRPD2/NRPE2) (Figure 2).
3.6. In Silico Expression Pattern of AGO, DCL, RDR, DRB, NRPD1, NRPE1, and NRPD2 Candidate Genes In order to analyze the transcript abundance of these 37 genes in different organs of common bean, we performed a comprehensive gene expression in silico analysis using online RNAseq data for genotype G19833. The results are shown in Figure 3. After moderated log-counts-per-million transformation, we applied hierarchical clustering (with Euclidean distance and Ward method) on the 37 genes. The genes can be organized into three clusters. Cluster 1 corresponds to genes presenting a low expression level. This cluster comprises several DRB genes (PvA_DRB4b, 1d, 1b, 1a), two AGO genes (PvA_AGO2a, 2b), one DCL gene (PvA_DCL2c), and one RDR gene (PvA_RDR1a). Cluster 3 corresponds to genes that are highly expressed and comprises four AGO genes (PvA_AGO1, 4c, 5, 4b), two DRB genes (PvA_DRB2a, 2b), one DCL gene (PvA_DCL1), as well as PvA_NRPE1 and PvA_NRPD2. In particular, PvA_AGO1 seems to be highly expressed in all tested organs. Finally, the remaining 20 genes correspond to genes presenting an intermediary expression level (cluster 2; Figure 3). For most genes of this cluster, the expression level seems relatively homogenous in the 11 analyzed organs, except PvA_RDR3, which seems up-regulated in the nodules.

Common Bean Pol IV and Pol V
In order to gain insight into the Pol IV and Pol V complex in the common bean genome, the largest and second-largest subunits of Pol IV and Pol V were searched by seeking AtNRPD1, AtNRPE1, and AtNRPD2/NRPE2 against common bean BAT93 and G19833 genomes with tBLASTn. Common bean encodes one NRPD1, one NRPE1, and one NRPD2/NRPE2, and hence they are named PvA/M_NRPD1, PvA/M_NRPE1, and PvA/M_NRPD2/NRPE2 (Table 1). These three proteins present a high level of identity (> 99%) between BAT93 and G19833. They are located on chromosome 2 (PvA_NRPD1), 11 (PvA_NRPE1), and 9 (PvA_NRPD2/NRPE2) (Figure 2).

In silico Expression Pattern of AGO, DCL, RDR, DRB, NRPD1, NRPE1, and NRPD2 Candidate Genes
In order to analyze the transcript abundance of these 37 genes in different organs of common bean, we performed a comprehensive gene expression in silico analysis using online RNAseq data for genotype G19833. The results are shown in Figure 3. After moderated log-counts-per-million transformation, we applied hierarchical clustering (with Euclidean distance and Ward method) on the 37 genes. The genes can be organized into three clusters. Cluster 1 corresponds to genes presenting a low expression level. This cluster comprises several DRB genes (PvA_DRB4b, 1d, 1b, 1a), two AGO genes (PvA_AGO2a, 2b), one DCL gene (PvA_DCL2c), and one RDR gene (PvA_RDR1a). Cluster 3 corresponds to genes that are highly expressed and comprises four AGO genes (PvA_AGO1, 4c, 5, 4b), two DRB genes (PvA_DRB2a, 2b), one DCL gene (PvA_DCL1), as well as PvA_NRPE1 and PvA_NRPD2. In particular, PvA_AGO1 seems to be highly expressed in all tested organs. Finally, the remaining 20 genes correspond to genes presenting an intermediary expression level (cluster 2; Figure 3). For most genes of this cluster, the expression level seems relatively homogenous in the 11 analyzed organs, except PvA_RDR3, which seems upregulated in the nodules.

Expression Pattern Analysis after Fungus Infection
In order to investigate the role of RNA silencing in pathogen defense in common bean, we studied the expression profile of seven genes, including AGO1, AGO2a, DCL2a, DCL2b, AGO4a, AGO4b, and AGO4c (indicated by the arrows in Figure 3). The expression of these genes in response to infection with the hemibiotrophic fungus C. lindemuthianum was studied using RT-qPCR at 6, 24, 48, 72, 96 hpi in a resistant genotype (incompatible interaction). Significantly, temporal gene expression analysis revealed that DCL2a and DCL2b are both ca. nine-fold up-regulated after infection compared to the mock control at 72 hpi. Similarly, a clear upregulation is observed at 72 hpi for both AGO4a and AGO2a, and also for AGO4b but with a lower extend (Figure 4). Conversely, the expression of AGO1 and AGO4c was not modified upon C. lindemuthianum infection (Figure 4). change values are shown on the left of the heat map. DAP = days after planting. Arrows indicate genes analyzed in RT-qPCR experiments after C. lindemuthianum infection.

Expression Pattern Analysis after Fungus Infection
In order to investigate the role of RNA silencing in pathogen defense in common bean, we studied the expression profile of seven genes, including AGO1, AGO2a, DCL2a, DCL2b, AGO4a, AGO4b, and AGO4c (indicated by the arrows in Figure 3). The expression of these genes in response to infection with the hemibiotrophic fungus C. lindemuthianum was studied using RT-qPCR at 6, 24, 48, 72, 96 hpi in a resistant genotype (incompatible interaction). Significantly, temporal gene expression analysis revealed that DCL2a and DCL2b are both ca. nine-fold up-regulated after infection compared to the mock control at 72 hpi. Similarly, a clear upregulation is observed at 72 hpi for both AGO4a and AGO2a, and also for AGO4b but with a lower extend (Figure 4). Conversely, the expression of AGO1 and AGO4c was not modified upon C. lindemuthianum infection (Figure 4).

Discussion
Several studies have pointed out that genes involved in silencing evolve rapidly with a great variation of number even between closely related species. However, our comprehensive analysis of various gene families involved in RNA silencing in two common bean genomes of contrasting origins allowed us to identify the same number of AGO (13), DCL (6), DRB (10), RDR (5), NRPD1 (1), NRPE1 (1), and NRPD2/NRPE2 (1) genes in both G19833 (Andean) and BAT93 (Meso-American) genomes, suggesting that no recent gene duplication/deletion occurred after gene pool divergence. Indeed, for each of the 37 genes analyzed in the present study, orthologs presenting a high percentage of protein identity (> 94%) were unambiguously identified between BAT93 and G19833 (Table 1, Figure 1). This suggests that even if the genome assembly of BAT93 is of lower quality compared to G19833, and does not allow repeated sequence analysis [55], this quality is sufficient for gene analysis. However, we can not exclude that some genes were missing, and higher quality genome assembly based on long-read sequencing will be needed to unambiguously address this question. These 37 genes are distributed on all Pv chromosomes, except chromosome 10. Importantly, our manual annotation led us to correct misannotated genes, in particular in BAT93 (Table 1). Even if no recent gene dynamics were identified after gene pool divergence, an interesting pattern of evolution was identified for DCL and AGO gene families in the common bean.
DCL genes, and in particular DCL2 genes, present a complex pattern of evolution in legume species. Unlike the single-copy genes of DCL1, DCL3, and DCL4, in Pv and Mt, there were three copies of DCL2. Soybean contains seven DCL genes in its ancient polyploid (paleopolyploid) genome [56] (Figure 1). In soybean, the most recent genome doubling event occurred approximately 9-14 million years ago, and the soybean genome maintains at least one gene duplicate for ca. 75% of its genes, termed homoeologous gene pairs [57]. GmDCL4a/GmDCL4b and GmDCL1a/GmDCL1b correspond to such homoeologous gene pairs, while GmDCL3 is present as a single copy. By contrast, GmDCL2a and GmDCL2b are locally duplicated, separated by 5 kb, on chromosome 9. In soybean, the age of this GmDCL2a/GmDCL2b duplication was estimated to be 19.4 Mya [56], indicating that it predates the whole genome duplication of soybean at 9-14 Mya [56], and the split for common bean and soybean 19 Mya [58,59]. Consequently, this strongly suggests that in the putative Pv Gm common ancestor, a locally duplicated pair of DCL2 genes were present. In agreement with this, we found in common bean two DCL2 genes, Pv_DCL2a and Pv_DCL2b, organized in tandem array in the corresponding syntenic region with soybean (Figures 2 and S1). In the Pv genome, an additional paralog, DCL2c, present on chromosome 8, was putatively derived from Pv_DCL2b by a yet unknown mechanism that could involve transposable elements identified in the vicinity of DCL2c. Three copies of DCL2 were also identified in Mt [14]. However, this amplification appears to be independent of that observed in the common bean ( Figure 1). In contrast, only a single DCL2 has been identified in various other legume species, including chickpea (Cicer arietinum), pigeonpea (Cajanus cajan), and each of the two genomes composing the allotetraploid groundnut (Arachis duranensis and Arachis ipaensis) [52]. These independent amplifications of the DCL2 genes in specific legume species could lead to their functional diversification and probably reflect their functional importance. In Mt, a nodule-specific role for DCL2 has been proposed [14], while in soybean, DCL2 genes regulate traits such as seed color via the production of 22 nucleotide siRNA from long inverted repeats [60]. In another study, GmDCL2 paralogs exhibited a wide range of transcriptional changes in response to stress, suggesting DCL2s may play an important role in stress response [56]. Congruent with these findings, we found that Pv_DCL2a and Pv_DCL2b, mildly expressed in most organs (Figure 3), are up-regulated at 72 hpi in leaves infected by C. lindemuthianum (Figure 4).
We identified 13 AGO genes in Pv, while 15 AGO genes were reported in a previous analysis performed on the G19833 genome. Indeed, compared to de Sousa Cardoso et al. [32], our manual annotation led us to discard one AGO10 and one AGO2 gene. Similarly, 13 AGO genes were also identified in both chickpea and pigeonpea [52]. Within angiosperms, several AGO subgroups have expanded differently in monocots and eudicots, with lineage-specific gene duplications [61]. For example, the grasses exhibit an expanded AGO 1/5/10 clade [17]. More precisely, maize and rice harbor many AGO5 paralogs, and a grass-specific AGO18 family, a deep branch of the AGO1/5/10 clade, has been discovered and played important roles during plant reproduction and viral defense [17]. In common bean, expansion of the AGO1/5/10 clade was also observed, but it was the result of AGO10 gene amplification since four AGO10 genes were identified in the common bean genome (Table 1). Similar amplification of AGO10 was also observed in soybean, where eight paralogs were identified in its paleopolyploid genome [17]. Each Pv AGO10 gene clearly corresponds to two Gm orthologs (Figure 1), strongly suggesting that AGO10 amplification occurred prior to the soybean/common bean divergence. In soybean, the expansion of the AGO10 family presumably co-evolved with the expansion of the miR165/166 family since 21 copies of miR165/166 are annotated in the soybean genome [17]. Likewise, expansions of miR165/166 genes, with 10 copies, have also been identified in the Pv genome (Geffroy V. and Meyers B.C.; unpublished results). In addition to AGO10, expansion was also observed for AGO2 (two members) and AGO4 (three members) in the Pv genome. AtAGO4 primarily binds 24-nt, repeat and heterochromatin-associated siRNAs, and functions in RNA-directed DNA methylation [62], while AtAGO2 functions in antibacterial immunity [63]. In common bean, AGO2 and AGO4 genes have non-redundant expression profiles (Figure 3), suggesting that they may have acquired divergent functions. PvA_AGO1 seems to be highly expressed in all common bean-tested organs. In agreement with our results, AGO1 expression is detected in many organs, such as leaves, roots, and flowers, in Arabidopsis [64], rice [12], B. napus [65], and the emerging medicinal model plant Salvia miltiorrhiza [66].
Functional analysis of genes involved in RNA silencing revealed that most of them play multiple roles, not only in growth and development but also in immune defense against pathogens [1,[67][68][69]. The importance of RNA silencing in plant viral defense has been well documented for a long time [63]. In addition to viral defense, more and more evidence is accumulating, showing that RNA silencing also plays a role in plant interactions with bacterial pathogens [70]. More recently, the potential role of RNA silencing in plant defense has also been reported for several fungal pathogens, including Verticillium dahliae [71], Verticillium longisporum [72], Magnaporthe oryzae [73], and Botrytis cinerea [74]. The importance of RNA silencing in plant defense is illustrated by the fact that it has stimulated a counter-defense system from the pathogens to overcome it. Indeed, it is now well-known that pathogens of a different nature (viruses, bacteria, fungi, oomycetes, and phytoplasma) have evolved effectors that are able to target and suppress the host plant RNA silencing pathway [67,[75][76][77][78]. Suppressors of RNA silencing were first discovered in viruses (VSRs, for viral suppressors of RNA silencing) [4]. At present, there is no evidence of putative suppressors of silencing acting in C. lindemuthianum. However, there is growing evidence that this is a common mechanism exploited by fungal pathogens to promote their infection [79]. Consequently, such suppressors probably exist in C. lindemuthianum, although not yet identified.
To investigate the contribution of some of the genes involved in RNA silencing in the defense response in common bean, we performed quantitative RT-PCR-based expression analysis on leaves of resistant bean plants inoculated with the hemibiotrophic fungus C. lindemuthianum (in an incompatible context). Whereas expression levels of PvA_AGO2a, PvA_AGO4a, and PvA_DCL2 (a and b) is low to moderate without any biotic stress (Figure 3), a strong up-regulation of these genes was observed mainly at 72 hpi ( Figure 4). On the contrary, expression of PvA_AGO1 and PvA_AGO4c, which are ubiquitously and highly expressed in uninfected plants, was not significantly modified after infection. This suggests that after fungal infection, PvA_AGO2, PvA_AGO4, and Pv_DCL2 may play a prominent role in the sRNA-based regulation of defense gene expression in the common bean. Interestingly, the Argonaute proteins, AGO4 and AGO2, have both been linked to antibacterial defense. AGO4, a component of the RdDM pathway that directs DNA methylation at specific loci, mediates resistance to P. syringae, independently of the other components of the RdDM pathway [80]. AGO2 functions in antibacterial immunity by binding a specific miRNA to modulate the exocytosis of antimicrobial PR proteins [81]. In the literature, different pathosystems involving either hemibiotrophic pathogens or incompatible plant-microbe interactions present similar results. Notably, in susceptible wild tobacco plants challenged by the hemibiotrophic fungus Fusarium brachygibbosum as well as in resistant cowpea plants in response to CPSMV (Cowpea severe mosaic virus) infection, an increased expression of AGO4 has been reported, whereas no change in expression was observed for AGO1 [82,83]. Moreover, an up-regulation of AGO2 expression after infection was reported in Arabidopsis after infection by the biotrophic bacteria P. syringae, in the oil crop B. napus infected by the fungal necrotrophic Sclerotinia sclerotiorum, and in the cowpea Vigna unguiculata infected by CPSMV [20,82,84]. Concerning the role of Dicer proteins in plant defense, little is known about DCL2, except that it is involved in the processing of viral dsRNA. However, it has also been observed that the quantity of DCL2 transcripts increases at the local site of infection by the CLRDV (Cotton leafroll dwarf virus) in a resistant genotype cotton Gossypium hirsutum [85]. This is in agreement with our results, where up-regulation of both PvA_DCL2a and Pv_DCL2b was observed in incompatible interaction with C. lindemuthianum. This suggests that in the common bean, an increased expression of specific genes involved in RNA silencing, acting in both the miRNA and the siRNA pathways, could counteract the infectious process of C. lindemuthianum. However, how these genes could regulate resistance in the common bean requires further investigation.

Conclusions
This work will further provide a solid foundation for future functional analysis of AGO, DCL, RDR, DRB, NRPD1, NRPE1, and NRPD2 genes in the common bean genome. For example, taking advantage of the work presented here, we would like to utilize virusinduced gene silencing (VIGS) to silence NRPD1 and NRPE1, in order to gain insight into the mechanisms involved in the unusual methylation pattern observed for NLR genes in the common bean [86][87][88]. Furthermore, our work will also help to design specific primers for RT-qPCR experiments. Finally, given the genomic location of the 37 genes studied (Table 1, Figure 2), and considering that RNA silencing is involved in a large number of traits, our work may also provide candidate genes for QTL analysis.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes13010064/s1, Figure S1: Syntenic analysis between common bean and soybean of the genomic region containing DCL2a and DCL2b genes. Table S1: List of primer sequences used in this study.