non-neutral, selected substitutions at the nucleotide level among a group of homologous sequences. However, the above methods do not identify functional groupings of residues that are nearby in the folded protein but dispersed in the primary amino acid sequence. Some studies of positive selection on LRRs have manually approx- imated structural proximity to improve their power (e.g., [48]). Other computational approaches, such as Consurf [49], optimal docking area (ODA) [50], and conserved functional group (CFG) analysis [45], do use protein structural information obtained experimentally (X-ray or NMR) or by homology modeling. CFG and Consurf search homologous input sequences for conserved groupings of amino acids on the surface of the folded protein. ODA models desolvation energy, i.e. it searches for continuous surface patches that undergo favorable energy change when buried during a modeled protein-protein association; a low ODA value indicates a location that is predicted to interact with a ligand or another protein. These latter programs offer valuable para- digms but are limited to proteins with a crystal structure or a reliable homology model, and CFG and Consurf have constraints related to handling of hydrophobic and/or repeat motif residues that limit their efficacy for analysis of LRR domains. In hundreds of important proteins, LRR domains contain the ligand specificity region or other functional sites, and there is significant interest in identification and manipulation of these sites. In the present study we developed the Repeat Conservation Mapping (RCM) program that predicts functional sites in LRR domains by identifying, among a group of homologous LRR- containing proteins, the patches of predicted spatially adjacent residues on the surface of the LRR that exhibit the greatest conservation or greatest divergence. The RCM method utilizes linear amino acid sequences as input, along with existing generalized LRR structure principles, to predict likely adjacency of residues in the LRR. It then identifies regional conservation scores for predicted surface residues based on conservation of that residue and its proximal surface amino acids. We validate the method using previously solved co-crystal structures for LRR with ligand, and through discovery and in vivo validation of previously unknown functional sites of the plant EFR and FLS2 MAMP receptors, for which structures are not currently available. The RCM program can be run using the publicly accessible web server at http://www.plantpath.wisc.edu/RCM, and the source code is openly available via a GNU general public license. Methods Repeat Conservation Mapping program RCM analyses can be conducted at a publicly accessible web interface served from a Linux-based virtual machine. The RCM web site runs a set of linked php files that draw upon PHP, Perl, HTML, Python and C scripts, including local implementations of functions from ClustalW2 [51,52], LRRScan [53] and MATLAB (The Mathworks, Inc., Natick, MA). The code is available under GNU general public license v3; local installation is not recommended for most end-users. Figure 1B and the following briefly describe the method. The linear amino acid sequences of two or more LRR domains are aligned and compared, and a conservation score is determined for each amino acid position. Characteristic properties of LRR do- mains are then used to generate a generic super-helical structural model of the LRR domain being queried, which places amino acids in their likely relative locations to each other in a folded protein. A sliding window analysis then determines a center- weighted conservation score for all possible groups of adjacent amino acids that are predicted to reside on the surface of the helix. Repeat Conservation Mapping Typically a 565 matrix of 25 amino acids, spanning 5 LRR repeats, is queried to derive each regional conservation score. A colored heat-map of the result directs investigator attention to the most conserved (or most divergent) regions. Alternative score- weighting systems for the sliding window analysis are available to users; adjustment of these values significantly impacts the resulting conservation map. The default weighting values achieve a useful balance between the excessive smoothing (loss of resolution) that occurs with less center-weighting, and the loss of site identification (similarity to the initially calculated individual residue conservation scores) that occurs with greater center-weighting. Readers are referred to a more detailed step-by-step description of conservation mapping in Text S1. Comparisons to other methods For positive selection analysis, sequence alignment and phy- logenetic trees were constructed using MEGA 4.0 [54]; informa- tion from these trees were used to calculate positive selection using the CODEML module in PAML [46]. Optimal docking analysis (ODA) was implemented as described in [55]. The current web module of Consurf was used [49]. Homology Modeling The FLS2 homology model was obtained using the structure of PGIP2 from Phaseolus vulgaris as a template (PDB ID 1OGQ) [56]. Both proteins (as well as EFR) belong to the plant extracellular LRR protein subfamily [1] characterized by the same consensus sequence in the LRR domain, i.e. xxLxLxxNxLt/ sGxIPxxLxxLxxL. Furthermore, the LRR domain is capped at the N-terminus and at the C-terminus by two small cysteine-rich domains, which are also evolutionarily conserved among PGIP and the RLK receptors FLS2 and EFR [57]. However, PGIP2 contains 10 LRRs matching the above consensus while FLS2 is characterized by 28 complete repeats, hence four separate alignments were manually prepared and used for homology modeling. The N-terminal and first 9 LRR repeats of PGIP2 were manually aligned to the N-terminal and first 9 LRR repeats of FLS2. The PGIP LRR domain was separately aligned to FLS2 repeats 7–15 and 13–21. A final alignment encompassed FLS2 repeats 20–28 plus the C-terminal flanking region. These align- ments were used, together with the relevant PGIP2 coordinates, to obtain four independent partial models using Modeller, version 9.1 [58]. Twenty models were obtained from each alignment and the lowest energy models were selected according to the Modeller objective function. The four partial models obtained by this strategy are partially overlapping with two or three LRRs in common between consecutive ones. This allowed superposition of the models in their overlapping regions and merging of coordinates to obtain a single full-length model using the SSM algorithm as implemented in the program COOT [59]. This model was further energy minimized within Modeller to obtain a final model that encompasses residues 1 to 744 of FLS2. The geometrical quality of the model was very good as judged with PROCHECK [60] with 98.1% of residues lying in allowed regions of the Ramachandran plot, 1.2% in generously allowed regions and only 0.6% (4 residues) in the non allowed regions. The EFR LRR domain was aligned to FLS2 using ClustalW and the alignment was manually adjusted, when necessary, to match the plant extracellular LRR consensus sequence. An EFR homology model, encompassing residues 1–576, was calculated and energy minimized within Modeller using the FLS2 homology model as a template. The final model has a very good geometry with 96.8% residues in allowed regions of the Ramachandran plot, 2.4% in PLoS ONE | www.plosone.org 3 July 2011 | Volume 6 | Issue 7 | e21614