Introduction 

Western white pine (Pinus monticola Dougl.) is an ecologically dominant, and highly important for forestry, species in the western North American ecosystems ([20]). In 1980, two P. monticola provenances from Jackson, Oregon were introduced in Greece in comparative species field trials of the Aristole University of Thessaloniki, Laboratory of Forest Genetics Experimental Centers located in Pertouli and Taxiarchis ([26], [27]). The aim of this P. monticola introduction was twofold: (a) the use of the species in plantation forestry for wood production, and (b) the evaluation of resistance to the white pine blister rust (Cronartium ribicola J.C. Fisch.). Trial results indicate large variation within provenance, growth comparable to that of local Abies borisii regis Mattf., but reduced survival ([28], [27]).

Cronartium ribicola is among the most destructive white pine fungal diseases which has caused severe mortality in North American white pines and especially in P. monticola, which is currently considered a species under threat ([21]). P. monticola species exhibits genetic variation in disease resistance to C. ribicola involving mainly a dominant disease resistance gene (Cr2 - [12], [19]). Disease resistance and maintenance of genetic diversity in production populations are very important aspects for P. monticola breeding ([4]). The development through breeding of trees resistant to fungi and the maintenance of a variety of genetic resources are of pivotal significance towards the substantial establishment of P. monticola reforestation and plantation forestry programs ([18]).

Classical breeding approaches, such as interspecific hybridization with the indigenous to Greece and resistant to C. ribicola Balkan white pine Pinus peuce Grisb., are time-consuming without warranties of specific results in a definite timeframe. Engagement of innovative molecular breeding approaches allows breeders to establish sustainable resistance more accurately and at a faster pace ([14]). PCR approaches implementing mainly degenerate primers were extensively applied in order to isolate members of the NBS family (nucleotide-binding site) / LRRs (Leucine-Rich Repeats) - containing genes, the most highly expanded group of genes linked directly to resistance (R) genes mediated functions in plants ([5]). This protein family is further classified into two major subfamilies based on structural characteristics of their putative amino-terminal signaling domain ([13]). One subfamily, TIR-NBS-LRR, contains an amino-terminal signaling domain that depicts homology with Toll and interleukin- 1 receptor (TIR), whilst the other subfamily, CC-NBS-LRR, instead of the TIR domain, carries a putative coiled-coil structure (CC) one ([13]). Both these subfamilies have been reported to mediate resistance to a wide range of phytopathenic fungi ([34]).

Resistance gene analogs (RGAs) are a large class of potential R-genes, being an important resource for the development of molecular markers, hence enhancing numerous disease resistance breeding efforts ([1], [30]). Moreover, these RGAs commonly contain tandem repeats such as LRRs which give rise to novel ligand-binding specificities under positive selection ([11]) via a number of different modes, such as point mutations, variations in repeat numbers and gene duplications ([8], [24], [37], [39]). Positive selection is defined as the rapid fixation of beneficial non-synonymous mutations and is an important evolutionary driving force for a number of paralogous genes towards their functional diversification, which might be selectively favored under biotic pressures in plants ([6]). Furthermore, the location of positively selected amino acid residues is an important aspect for obtaining novel gene functions ([25]). Previous studies have reported that solvent-exposed regions of the LRRs repeats are subjected to positive selection; such evidence has been interpreted as a sign of the involvement of these regions in pathogen recognition ([29], [39], [11]). Overall, the evolution patterns of plant NBS-LRR R genes and of their RGAs counterparts is likely a complex process ([40], [31]).

The TIR-NBS-LRR RGAs from P. monticola were previously investigated for their diversity and expression, in order to gain insights regarding their potential usage in disease resistance breeding programs ([17]). In this study, we address the question whether signatures of positive selection are systematically employed as an evolutionary driven force acting on TIR-NBS-LRR-encoding RGAs genes at P. monticola. Thus, we gain insights into the evolutionary profiles of the P. monticola TIR-NBS-LRR-containing RGAs genes, hypothesizing that those successive episodes of positive selection might contribute in the acquisition of novel pathogens recognition repertoires. Overall, our findings enrich and support the hypothesis that these P. monticola RGAs could be employed as biomarkers for both, monitoring the ongoing and, more effectively applying future P. monticola disease resistance breeding efforts.

  Materials and methods 

P. monticola TIR-NBS-LRR RGAs mining and phylogenetic analysis

All 33 P. monticola RGAs partial proteins previously identified by Liu & Ekramoddoullah ([17]) were downloaded from the GenBank by keyword searches based on their NB-ARC domain profile (PF00931.19/ CL0023). Using these amino acid sequence IDs, the respective ORF nucleotide sequences were also retrieved. All amino acid sequences had to fulfill the following criteria in older to be evaluated for evolutionary signatures acting upon them: (a) existence of valid structures by means of protein lengths up to 90 residues, positions of NB-ARC domain motifs, blast(p) hits against known plant resistance genes; and (b) verification of functional domains using the predictions against the Pfam database (⇒ http:/­/­pfam.­sanger.­ac.­uk/­search#tabview= tab1) and the InterProScan 5 software ([9]).

The phylogenetic relationships among these P. monticola RGAs amino acid sequences were revealed by aligning them using the MUSCLE software ([7]) and performing a tree reconstruction using the RAxML software ([32]) inferring a gamma model of rate heterogeneity. All these analyses were performed on the Geneious R7 platform ([10]).

P. monticola TIR-NBS-LRR RGAs evolutionary analysis

The YN00, CODEML and CODEMLSITES programs were employed from the PAML version 4.8 software package ([36]) in order to assess evidence of positive selection signatures acting in these P. monticola TIR-NBS-LRR RGAs. These tests compare the substitution rates of non-synonymous (dN) to synonymous (dS) mutations, and estimate the derived omega (ω) values, both among the nucleotide codon sequences and along the branch lineages. The rates of dN/dS nucleotide substitutions per site across all possible pairwise sequence comparisons were estimated using the counting approach ([35]) implemented in the YN00 program. The total dN, dS, and ω values for all pairwise sequence comparisons among these RGAs are shown in Tab. S1 (Supplementary material).

An NJ tree at amino acid level and the respective nucleotide alignment were both used in the CODEML and CODEMLSITES analyses according to Lynn et al. ([22], [23]); details of the exact methods for these analyses are found in Zambounis et al. ([39]). Briefly, the amino acid sequence alignment was conducted by employing the MUSCLE software ([7]), whilst the construction of the respective phylogenetic tree was performed using the MEGA 5 software ([33]). A Poisson correction was employed with a gamma distribution modeling among sites and rates variations ([39]). Finally, the Baseml software (PAML version 4.8 - [36]) was implemented for estimating the alpha shape parameter using the REV model.

CODEML was tested in P. monticola TIR-NBS-LRR RGAs dataset for estimating the variable selective pressures acting among lineages in the phylogeny, whilst CODEMLSITES was employed to test for site-specific codon substitution models allowing the identification of selective signatures across amino acid residues (Tab. 1). In both programs log-likelihood calculations were performed for each model and comparisons between each other were implemented by likelihood ratio tests (LRTs). The empirical Bayesian method further allowed the detection of codon residues subjected to positive selection by higher posterior probabilities. All data sets as well as the outputs of CODEML and CODEMLSITES programs are available upon request.

  Results 

P. monticola TIR-NBS-LRR RGAs structure and phylogenetic analysis

The 33 P. monticola TIR-NBS-LRR RGAs studied were found to share Pfam hits with the NB-ARC domain (PF00931.19) and the InterProScan ones with LRR-containing proteins (PTHR23155). The RGAs were corresponding to regions between conserved kinase-2 and P-loop domains. A RAxML phylogenetical analysis at the amino acid level was conducted in order to explore their expansion (Fig. 1). The backbone topology of this tree was properly resolved with 63 branches (Fig. 1), whilst the overall pairwise identity among the amino acid sequences was 59.6%, indicative of a high degree of divergence.

Fig. 1 - RAxML phylogenetic reconstruction of the amino acid sequences of P. monticola TIR-NBS-LRR RGAs. In various colors the six clades which were being assigned whilst, asterisks (*) indicate tree branches under positive selection (ω values higher than 1) with posterior probabilities using the CODEML method.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Positive selection acting across P. monticola TIR-NBS-LRR RGAs lineages

In order to verify whether accelerated evolution and signs of positive selection might have also contributed to the divergence of these TIR-NBS-LRR RGAs genes, our RGAs dataset was subjected to several successive tests of positive selection using the counting (YN00) and the ML (CODEML/ CODEMLSITES) methods. A robust and reliable nucleotide alignment resulted, based on the presence of zero or few gaps across the aligned sequences. The dN and dS mutation values across the entire RGAs sequences were calculated using the counting YN00 method (see Tab. S1 in Supplementary material). Any substitution saturation effects were checked by calculating dS rates between the aligned nucleotide sequences. The average numbers of synonymous (dS) mutations per synonymous sites over all sequence pairs were below the value of 2, the crucial threshold value above which the corresponding sequences would have to be excluded from further evolutionary analyses in order to bypass any saturation effects based on nucleotide substitutions ([35]).

Finally, the CODEML approach was employed in order to identify branches subjected to positive selection. Statistically significant evidence of positive selection (ω > 1.00) was detected in numerous branches (Fig. 1) with posterior probabilities P < 0.0001. These signs were found across clades I, IV and V acting with a non-general rule, as positive episodes were evidenced both in the ancient and terminal branches (reflecting probably a series of recently successive bursts of gene duplications).

Positive selective pressures acting among amino acid residues across the P. monticola TIR-NBS-LRR RGAs

Possible positive selection at the individual amino acid residues level was investigated by the CODEMLSITES method ([36]) which compared various models (M0 against M3, M1 against M2, and finally M7 against the stringent M8) in our RGAs dataset. Two models (M3 and M8) revealed positive selective pressures acting among amino acid residues under significant posterior probabilities (P < 0.0001 and P < 0.05, respectively - Tab. 1). These extensive signs of positive selection were found acting rather widely upon amino acid residues and across the respective tree branches. The accurate locations of these positively selected sites (with ω values up to 1.48) are depicted in Tab. 1.

  Discussion 

Conifer genomes are in general large with DNA content reaching more than 26.500 Mbp per haploid genome ([2] and references therein). Taking also into consideration their long generation times, classical breeding, especially for disease resistance is a very slow process. On the contrary, molecular breeding approaches may provide promising and effective strategies for both R genes and RGAs counterpart identification based on their highly conserved motifs ([20]). Previously, it was postulated that there is a highly diverse group of TIR-NBS-LRR genes in P. monticola ([17]) along with the existence of numerous CC-NBS-LRR-containing proteins ([20]). Similar abundance and diversity of these genes was also observed in other plant genomes ([40], [3], [31]). Hitherto, TIR-NBS-LRR RGAs in P. monticola exhibited high variation at amino acid levels, whilst almost half of them were found to be pseudogenes representing a rather complex TIR-NBS-LRR subfamily ([17]).

In our study, the analysis of 33 P. monticola TIR-NBS-LRR RGAs showed a high level of diversity that may contribute to their evolutionary fitness, and the putatively functional character of their transcription activity ([17]). These RGAs are a potential source for the development of molecular markers, especially regarding selection of resistance to white pine blister rust. Furthermore, RGAs were found to share conserved structures rendering them ideal candidates for applying evolutionary bioinformatics approaches.

The CODEML analysis, showed recent episodes of positive selection acting extensively in clades I, IV and V, which in turn were overlapping similar more ancient events. Similar positive selective pressure acting on plant RGAs has also been reported in other surveys ([15], [34], [30], [11]). In principle, this reservoir of genetic variation could drive the evolution of new R-gene specificities ([30], [31]). To our knowledge this is the first systematically report of evolutionary pressures acting in P. monticola RGAs. Positively selected RGAs branches were not only at terminal positions, implying rather recent positive selective signatures, but also at ancestral positions on the phylogeny. This finding is in accordance with previous reports which showed that numerous selective episodes occur even from the birth of a paralogous locus ([37], [39]).

Signs of positively selected residues were found to be distributed rather evenly across the RGAs amino acid sequences. Positive selection on NBS-containing gene families in plants and especially among their LRR domains, were also reported in various plant hosts ([3], [16], [39], [38], [30], [11]). These findings support the view that selection for durable disease resistance may be a crucial component of plant breeding programs. Additionally, comparative analyses of Rosaceae trees regarding their R genes revealed that solvent-exposed residues of the LRRs domains are hyper-variable with intensive positive selective pressures acting on them ([30]). Moreover, such evidence is consistent with host-pathogen co-evolution processes leading to the acquisition of novel resistance functional specificities ([39], [30]). These findings would facilitate the development of informative RGA-derived molecular markers for identifying fungal resistant P. monticola and P. monticola × P. peuce genotypes in breeding programs using segregating progenies. They may also prove to be of importance for P. peuce itself, as despite the fact that Balkan white pine is a resistant species, future selection pressures on C. ribicola in light of climatic change may alter the virulence of the fungus towards its host.

  Conclusions 

The P. monticola RGAs gene homologues encoding for R gene members of the TIR-NBS-LRR subfamily were found to be subjected to intensive positive selection, acting widely in three clades across their phylogeny, both on terminal and ancestral lineages. The majority of positively selected residues sites were found to be localized widely across these RGAs sequences, putatively affecting the structures of their ligand-binding domains. Practical applications in the development of disease resistant plants through breeding applications may also be considered, facilitating for instance, the development of informative RGAs-derived sequence characterized amplified region (SCAR) markers, based on RGAs sequences with extensive positive selection acting upon them. Such markers might be particularly useful for identifying fungal resistant genotypes in routine marker-assisted selection (MAS) breeding programs.

  List of Abbreviations 

The following abbreviations have been used throughout the paper:

1. CC: N-terminal coiled-coil domain;
2. LRR: Leucine Rich Repeat;
3. NBS: nucleotide-binding site domain;
4. PAML: Phylogenetic Analysis by Maximum Likelihood;
5. RaxML: Randomized Axelerated Maximum Likelihood;
6. RGA: Resistance Gene Analog;
7. TIR: Toll/interleukin-1 receptor domain;
8. WWP: Western white pine.

  Acknowledgments 

The co-authors have not a conflict of interest to declare.

  References

  Supplementary Material

Tab. S1 - YN00 method results of the P. monticola TIR-NBS-LRR RGAs.