Analysis of TIR- and non-TIR-NBS-LRR disease resistance gene analogous in pepper: characterization, genetic variation, functional divergence and expression patterns

Identification of non-TIR- and TIR-NBS-LRR CaRGAs in pepper

Candidate non-TIR- and TIR-NBS-LRR CaRGAs were identified in pepper using two approaches, PCR amplification with degenerate primers and database mining. Two pairs of degenerate primers, previously designed based on conserved domains (P-loop and GLPL regions) among known NBS-LRR R genes from other plant species [37–39], were used. Two bands of the predicted size (~500 bp) were amplified using the genomic DNA of pepper (Additional file 1). The bands were then excised from agarose gels and cloned. A total of fifty clones were randomly selected for sequencing, twenty-four of which were highly homologous to NBS-LRR sequences or known R genes from other plant species. These sequences were designated as NBS-LRR CaRGAs. The remaining clones were homologous to either a putative polyprotein or a hypothetical LRR protein. This finding suggests that degenerate PCR amplification is a very effective method for isolating potential CaRGAs from pepper.

A total of fifty-four CaRGAs were identified using the key word ‘Capsicum resistance gene’ in a search of the National Center for Biotechnology Information (NCBI) non-redundant protein database (http://www.ncbi.nlm.nih.gov/). A total of seventy-eight CaRGAs were obtained using these two methods. Among the seventy-eight CaRGAs, fifty-five had uninterrupted open reading frames (seventeen from amplified products, thirty-eight from data mining) and twenty-three CaRGAs had stop codons in the reading frames. These CaRGAs may be non-functional genes. Hence, they were excluded from further analysis. Moreover, two R genes reported previously were selected for further analysis [32, 33]. In addition, more than one hundred and seventy and four hundred NBS-LRR R genes from Arabidopsis and rice [3, 26], respectively, were used to retrieve potential R genes or CaRGAs from pepper. However, no novel CaRGAs or R genes were found. The nucleotide and amino acid sequences of all these pepper CaRGAs are listed in Additional file 2.

Sequence analysis and phylogenetic relationship between non-TIR- and TIR-NBS-LRR RGAs

BLASTX searches revealed that the fifty-seven CaRGAs had a certain degree of identity with known R genes as well as some RGAs from other plant species. They contained the conserved NB-ARC (nucleotide-binding adaptor shared by Apaf-1, R proteins and Ced-4) domain of known R genes [40–42]. Further analysis indicated that these CaRGAs, except for CaRGA54 (AF513549), CaRGA55 (FJ605104), CaRGA56 (FJ605105) and CaRGA57 (FJ605107), also included typically conserved motifs such as P-loop, kinase-2 and GLPL.

To investigate the evolutionary relationships among the NBS-LRR CaRGAs of pepper, phylogenetic trees were constructed based on the region between the P-loop and GLPL motifs using Molecular Evolutionary Genetics Analysis (MEGA) 5.0 software [43]. The aforementioned four CaRGAs were excluded from subsequent analyses because they lacked conserved P-loop and GLPL motifs. Thus, a total of fifty-three CaRGAs of pepper were used for phylogenetic tree construction. Six known R genes (RPM1, Gpa2, L6, M, N and Prf), downloaded from the GenBank database, were also included in the analysis. The phylogenetic tree is shown in Figure 1. Two primary groups, non-TIR- and TIR-NBS-LRR, were obviously distinguished, which is consistent with previously described patterns [9, 15]. The non-TIR-NBS-LRR subfamily included thirty-eight members (CaRGA13–CaRGA43, CaRGA45–CaRGA48, CaRGA51, BS2 and CaMi), whereas only fifteen members (CaRGA01–CaRGA12, CaRGA49–CaRGA50 and CaRGA44) belonged to the TIR-NBS-LRR subfamily. These results support the view that both TIR-NBSLRR and non-TIR-NBS-LRR R genes occur in dicot species [15]. Within the non-TIR-NBS-LRR subfamily, CaRGAs could be separated into four subgroups, which were designated CaRGA I–IV. In the TIR-NBS-LRR subfamily, the CaRGAs could be divided into three subgroups, which were designated CaRGA V–VII. Each of these seven subgroups included different numbers of CaRGA members. The CaRGA III, CaRGA IV, CaRGA V and CaRGA VI subgroups contained two to four members, whereas the CaRGA I, CaRGA II and CaRGA VII subgroups consisted of eight to nineteen members. This result reflects a difference in the abundance of these RGA subgroups in the pepper genome (Figure 1). The characteristics of these non-TIR- and TIR-NBS CaRGA sequences from pepper are shown in Table 1.

Multiple sequence alignments of non-TIR- and TIR-NBS-LRR RGAs in pepper

The CaRGAs identified from pepper were obviously separated into non-TIR- and TIR-NBS-LRR subfamilies (Figure 1). Consequently, multiple sequence alignments of non-TIR- and TIR-NBS-LRR RGAs were performed separately. The sequence alignment between non-TIR-NBS-LRR RGAs and known disease R genes, including RPM1, Gpa2 and Prf, revealed six conserved motifs (i.e. P-loop, RNBS-A-non-TIR, kinase-2, RNBS-B, RNBS-C and GLPL) (Figure 2). Similarly, alignment analysis with known R genes N, M and L6 revealed that the TIR-NBS-LRR RGAs also contain six conserved motifs (i.e. P-loop, RNBS-A-TIR, kinase-2, RNBS-B, RNBS-C and GLPL) (Figure 3). The presence of these consensus motifs within the non-TIR- and TIR-NBS-LRR CaRGAs provides further evidence that the cloned sequences are NBS-LRR gene family members.

Analysis of functional divergence between the non-TIR- and TIR-NBS-LRR CaRGA subfamilies

Type I and II functional divergence between the non-TIR- and TIR-NBS-LRR subfamilies of pepper was assessed by posterior analysis using DIVERGE 2.0 software, which evaluates the (site-specific) shifted evolutionary rate after gene duplication or speciation [44, 45]. Posterior analysis results in a site-specific profile for predicting important amino acid residues responsible for functional divergence. The estimation was based on multiple sequence alignments and a neighbor-joining (NJ) tree of the NBS domain of the pepper RGAs, with clear separation of the two different subfamilies, non-TIR- and TIR-NBS-LRR (Figure 1). The coefficient of type-I functional divergence (θ_I) between the non-TIR- and TIR-NBS-LRR subfamilies was significantly greater than 0 (θ_I = 0.533 ± 0.156, P < 0.05). This result suggests that the altered functional constraint between the subfamilies is statistically significant and that some amino acid sites are subjected to different site-specific shifts in evolutionary rate that can lead to a subfamily-specific functional evolution after diversification. Compared with the findings for type-I functional divergence, the coefficient of type-II functional divergence (θ_II) between the non-TIR- and TIR-NBS-LRR subfamilies was less than 0. Therefore, type-I functional divergence was the primary pattern for the evolution of the non-TIR- and TIR-NBS-LRR subfamilies in pepper.

We further estimated the critical amino acid residues responsible for the functional divergence by calculating the site-specific profile based on a posterior probability (Q _k ) analysis of the non-TIR- and TIR-NBS-LRR subfamilies. Among all of the aligned sites, the Q _k values of most sites were <0.5 (Additional file 3). To reduce false positives, Q _k  > 0.70 was used as a cut-off to identify critical amino acid residues associated with type-I functional divergence between the non-TIR- and TIR-NBS-LRR subfamilies. A total of thirteen sites (positions 21, 22, 23, 69, 84, 104, 113, 116, 118, 125, 146, 155 and 165) were predicted (Additional file 3). Among these sites, the Q _k value of site 22 was 0.878, which was predicted to be highly related to functional divergence, whereas the degree of relation to functional divergence was lowest at site 155 (Q _k  = 0.729).

We also found that the degree of conservation of critical amino acid residues differed between non-TIR- and TIR-NBS-LRR RGAs. For example, the amino acid residue at site 22 was valine (V) and leucine (L) in the non-TIR- and TIR-NBS-LRR RGAs, respectively. Nevertheless, some amino acid residues were highly conserved between the subfamilies, such as those at sites 84 (E), 133 (F) and 155 (W) from the non-TIR-NBS-LRR subfamily and sites 23 (S), 116 (W), 118 (G) and 146 (V) from the TIR-NBS-LRR subfamily (Figure 4). These results suggest that these critical amino acid residue sites are subjected to strong functional constraints within the non-TIR- and TIR-NBS-LRR subfamilies.

Comparing evolutionary rates among NBS-LRR RGAs in pepper

Sliding window analysis

A sliding window analysis of the Ka/Ks ratios was carried out alongside the determination of the nucleotide sequences of the RGAs for pairwise comparison of any two groups using CRANN software to explore important regions that may have contributed to the functional diversification of NBS-LRR RGAs in pepper [48]. Within the non-TIR-NBS-LRR subfamily, we observed two groups of pairwise comparisons that showed similar tendencies. One was CaRGA I/CaRGA IV and CaRGA II/CaRGA IV, the other was CaRGA II/CaRGA III and CaRGA III/CaRGA IV. The remaining two groups showed different patterns. Overall, the sliding window profile revealed at least three regions of high peaks that were under positive selection, given that all Ka/Ks values were >1 (Figure 5A). Within the TIR-NBS-LRR subfamily, we observed 2–3 regions of high peaks with high Ka/Ks values, consistent with positive selection in three pairwise comparisons (Figure 5B). These regions may be exposed to strong functional constraints. Two of these regions may be located in the RNBS-B and RNBS-C domains, respectively. Therefore, we inferred that the observed difference in the regions of functional constraint may reflect the functional specificities of NBS-LRR RGAs in pepper.

Determination of functional specificity positions among pepper NBS-LRR RGAs

Expression analysis of CaRGAs in different organs and in response to defence signaling molecules

To analyze the expression levels of CaRGAs in the different organs of pepper, two representatives of each class were randomly used for expression analysis using reverse transcription (RT)-PCR. CaRGA14 and CaRGA15 were the only two members of Class CaRGAIV. Given that the sequence similarity among these RGAs was high, only one pair of specific primers was designed using Primer 5.0 software. A total of thirteen pairs of specific primers were obtained. As shown in Figure 6A, thirteen CaRGAs were expressed in different plant organs, but their expression levels were different. Among these, CaRGA04, CaRGA13 and CaRGA38 were expressed in the leaves, stems and roots at low levels, whereas CaRGA01, CaRGA03, CaRGA14 and CaRGA18 were expressed at relatively high levels. The remaining CaRGAs were expressed at intermediate levels. No expression of CaRGA44 and CaRGA49 genes was detected in the roots. The CaRGAV and II subgroups showed similar expression patterns.

After treatment with abscisic acid (ABA) and salicylic acid (SA), the expression levels of most of these CaRGA genes changed. We found that the addition of ABA increased the transcription levels of several of the analysed CaRGA genes, namely CaRGA18, CaRGA51, CaRGA23, CaRGA14 and CaRGA49. By contrast, ABA addition decreased the expression levels of CaRGA13, CaRGA01 and CaRGA05. However, the expression of CaRGA03, CaRGA36 and CaRGA38 were unchanged. SA is known to play a vital role in plant defence against pathogens [50]. SA also induces the expression of a range of pathogen defence genes in plants [51, 52]. Among the CaRGAs tested, CaRGA01, CaRGA05, CaRGA03, CaRGA49, CaRGA14, CaRGA51, CaRGA36 and CaRGA38 displayed the most marked responses to SA treatment (Figure 6A). The expression levels of CaRGA18 and CaRGA23 decreased. The remaining genes showed no response to the SA treatment conditions. Real time-PCR experiments confirmed the expression levels of the selected CaRGAs (Figure 6B).

We found that two members from a subgroup may have different expression patterns, such as CaRGA01 and CaRGA04. However, similar expression patterns were also observed (CaRGA44 and CaRGA49). We also found that the CaRGAIV and VI subgroups had similar expression patterns. However, whether these genes with similar expression patterns have similar functions remains unclear. In addition, some earlier studies have reported that signaling molecules not only function as a critical signal for downstream resistance events but also upregulate the expression of R genes [51–56]. Some CaRGA genes were activated by SA and ABA (Figures 6A and 6B). This suggests that these stimuli induce the expression of the CaRGA genes and that they may play a potential role in mediating cross-talk between signaling pathways.

In summary, this paper provided detailed characterization and data on the functional divergence of non-TIR- and TIR-NBS-LRR CaRGAs in pepper. The mode of selection (positive selection, purifying selection and neutral selection) among the non-TIR- and TIR-NBS-LRR CaRGA subfamilies was identified by Ka/Ks analysis. However, the kind of evolutionary mechanisms responsible for the evolution of R genes in pepper cannot be inferred with certainty without the complete set of NBS-LRR genes from the pepper genome. Future studies must focus on verifying and elucidating the biological function of these CaRGA genes using supplementary experimental approaches, particularly with virus- or Agrobacterium-mediated transient assays [57] or by performing loss-of-function experiments, such as virus-induced gene silencing [58].

Plant material

The sweet pepper breeding line PBC631B was selected to isolate potential NBS-type disease R genes. PBC631B seeds were germinated, and the seedlings were grown in growth chambers at 25°C for 12 h (day) and 18°C for 12 h (night). Relative humidity was maintained at 65–75%. Young leaves were harvested from 4 week-old plants, immediately frozen in liquid nitrogen and then stored at −80°C for nucleic acid extraction. For hormone treatments, the seedlings were cultured in Hoagland’s solution containing 100 μM SA and 100 μM ABA for 6 h. The treated samples were then harvested for testing. Genomic DNA was isolated using a commercial plant DNA extraction kit (Bioteke, Beijing, China) according to the manufacturer’s instructions.

Degenerate primers and PCR amplification

Two pairs of degenerate primers were selected for isolating the RGA sequences in pepper (Additional file 5). The positions of the degenerate primers were identified in the two conserved domains (P-loop and GLPL, respectively) in the plant R genes (Figure 7). PCR amplifications were performed using a PTC-100 thermal cycler. PCR amplification reactions were conducted in a total volume of 25 μL containing 20 ng of template DNA, 2 μL of 10× PCR Buffer, 1.5 μL of 25 mmol MgCl₂, 2 μL of dNTPs (2 mmol/L), 1 μL of Primer-F (10 μmol/L), 1 μL of Primer-R (10 μmol/L), 1 unit of ExTaq DNA (5 U/μL) and 15.3 μL of ddH₂O. The PCR reaction cycle was as follows: denaturation at 94°C for 4 min, followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 60 s, followed by a final extension at 72°C for 5 min.

Cloning and sequencing of PCR products

Full volumes of the PCR products were run on a 1.0% agarose gel. Bands of the expected size (~500 bp) were excised from the agarose gel and purified using a DNA gel purification kit (Sangon, Shanghai, China). The obtained DNA was cloned into a pGEM T-Easy vector (Promega, Madison, WI, USA) and transformed into competent Escherichia coli JM 10⁹ cells according to the manufacturer’s instructions. The cloned DNA fragment was sequenced by Bio-Asia Company (China).

Collection of other NBS-LRR RGAs through database mining

Other NBS-LRR RGAs from pepper were identified in the NCBI non-redundant protein database (http://www.ncbi.nlm.nih.gov/) using the key word ‘Capsicum resistance genes’. A total of 54 NBS-LRR-type sequences were identified. Six known disease R genes, including RPM1 (X87851), Gpa2 (AF195939), L6 (U27081), M (U73916), N (U15605) and Prf (U65391) were downloaded from the GenBank database, and their phylogenetic relationship with pepper NBS-LRR disease R genes was determined.

Sequence analysis and phylogenetic tree construction

Each of the acquired DNA sequences was trimmed of vector sequence contamination using VecScreen at NCBI. Identity and similarity searches of nucleotide and amino acid sequences were performed using BLAST at the NCBI GenBank database (http://www.ncbi.nlm.nih.gov/BLAST/). Sequence alignments were carried out using Clustal W (BioEdit software) [59]. The phylogenetic tree was constructed by the NJ method using MEGA 5.0 software [43]. The reliability of the interior nodes was assessed using 1000 bootstrap replicates. Human apoptosis activating factor-1 (Apaf-1), which contains homologous motifs with the NBS region in plant disease resistance genes, was included in the phylogenetic analysis as an outgroup sequence [60].

Analysis of functional divergence

The phylogenetic tree of the pepper NBS-LRR RGAs was broadly grouped into two clusters, namely the TIR- and non-TIR-NBS-LRR subfamilies. DIVERGE 2 software [61] was used to evaluate the potential functional divergence and to predict the important amino acid residues in these two subfamilies. The coefficients of type I and II functional divergence (θ_I and θ_II) between these two groups of pepper NBS-LRR RGAs were estimated through posterior analysis. A θ_I or θ_II value significantly >1 indicates altered selective constraints or a radical shift in amino acid physiochemical properties after gene duplication and/or speciation [44, 45]. A site-specific posterior analysis (Q _k ) was also used to predict amino acid residues important for functional divergence.

Calculation of Ka/Ks ratios

We detected the mode of selection (positive selection, purifying selection or neutral selection) among the non-TIR- and TIR-NBS-LRR RGA subfamilies. The Ka/Ks ratios were calculated according to Nei and Gojobori [62] using K-Estimator 6.0 software [63, 64]. We identified the NBS-LRR RGAs subject to different selection pressures. DnaSP 5.0 software [65], which can calculate pairwise distance as part of a sliding window analysis, was applied. Ka/Ks values were plotted using Microsoft Excel to produce a graph of Ka and Ks values. The resultant Ka and Ks values were the sum of every possible pairwise comparison between every subgroup of R gene candidates selected for that particular window.

Analysis of SDPs

SDPfox was used to predict the SDPs that may determine the functional specificity of homologous proteins [49]. SDP presents the statistical significance of the predictions in the form of Z-scores (the number of standard deviations away from the expected value) and displays the most significant positions in a multiple sequence alignment. Positions with high Z-scores are predicted to determine functional specificity.

RNA isolation, DNase l treatment, cDNA synthesis and semi-quantitative RT-PCR analysis

Total RNA from all samples was isolated using TRIZOL reagent according to the manufacturer’s protocol (Invitrogen). RNA integrity, RNA concentration, RNA quality, DNase l treatment and cDNA synthesis were performed as previously described [66]. Two representatives of each class of NBS-LRR R genes were selected for expression analysis. Gene-specific primer pairs were designed using Primer5.0 software. A total of 13 pairs of CaRGA-specific primers were obtained (Additional file 5). Nine CaRGAs (i.e. CaRGA23, CaRGA38, CaRGA18, CaRGA36, CaRGA51, CaRGA13, CaRGA14, CaRGA49 and CaRGA44) were selected from NCBI. The remaining four CaRGAs (i.e. CaRGA03, CaRGA05, CaRGA04 and CaRGA01) were selected from PCR amplification. Subsequently, we analyzed the visualization of amplicon fragments to verify whether these primers were specific. Primers that exhibited the electrophoresis pattern of a single amplicon with the correct predicted size were considered CaRGA-specific primers. RT-PCR reactions were carried out using an Eppendorf PCR system 5331 cycler. The cycling program was as follows: 10 min at 94°C, 30 cycles of 45 s at 94°C, 45 s at 55°C and 1 min at 72°C and a 7 min extension at 72°C. The glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was used as a reference [66].

Quantitative RT-PCR and data analyses

Primer specificity for quantitative RT-PCR was further confirmed by analyzing melting curves. Primers corresponding to the melting curves that yielded single sharp peaks were used for quantitative RT-PCR analysis. Real-time PCR reactions were carried out in a total volume of 25 μL containing 12.5 μL of 2× SYBRGreen PCR MasterMix (Applied Biosystems), 1 μL of each primer, 1 μL of template (10× diluted cDNA from samples) and 9.5 μL of sterile distilled water. The thermal conditions were as follows: 95°C for 10 min, followed by 40 cycles at 95°C for 15 s and a final step at 60°C for 1 min. Quantification analysis was performed through the comparative CT method. All reactions were performed in triplicate in 96-well reaction plates using an iQ5 machine (Bio-Rad). Two independent replicates were performed. GAPDH was used as a reference gene for the expression analysis of the pepper CaRGA genes [66].