The 33 species analyzed represent a broad spectrum of land plants (embryophyta), with one bryophyta genome, Physcomitrella patens (PHYPA, moss) (Rensing et al., 2008), the Selaginella moellendorffii genome (SELML, spikemoss) (Banks et al., 2011), representative of the lycopodiopsida, and 31 species of magnoliophyta (angiosperms), divided into eight monocot species (Phoenix dactylifera (PHODC, date palm) (Al-Dous et al., 2011), Musa acuminata (MUSAC, banana) (D’Hont et al., 2012), two subspecies of Oryza sativa (rice), Oryza sativa ssp. japonica (ORYSJ) (Goff et al., 2002) and Oryza sativa ssp. indica (ORYSI) (Yu et al., 2002; International Rice Genome Sequencing Project, 2005), Brachypodium distachyon (BRADI, purple false brome) (The International Brachypodium Initiative, 2010), Zea mays (MAIZE, corn) (Schnable et al., 2009), Sorghum bicolor (SORBI) (Paterson et al., 2009) and Setaria italica (SETIT, Foxtail millet)) (Bennetzen et al., 2012; Zhang et al., 2012), and 23 dicot species (two Solanum species, Solanum tuberosum (SOLTU, potato) (The Potato Genome Sequencing Consortium, 2011) and Solanum lycopersicum (SOLLC, tomato) (Tomato Genome Consortium, 2012), Vitis vinifera (VITVI, Grape Vine) (Jaillon et al., 2007), Lotus japonicus (LOTJA) (Sato et al., 2008), Medicago truncatula (MEDTR, Barrel Medic) (Young et al., 2011), Glycine max (GLYMA, soybean) (Schmutz et al., 2010), Cajanus cajan (CAJCA, pigeon pea) (Varshney et al., 2012), Prunus persica (PRUPE, peach) (Ahmad et al., 2011), Malus x domestica (MALDO, apple) (Velasco et al., 2010; Jung et al., 2012), Ricinus communis (RICCO, castor oil plant) (Chan et al., 2010), Jatropha curcas (JATCU) (Sato et al., 2011), Manihot esculenta (MANES, Cassava) (Prochnik et al., 2012), Populus trichocarpa (POPTR, black cottonwood) (Tuskan et al., 2006), two Cucumis species, Cucumis sativus (CUCSA, cucumber) (Huang et al., 2009) and Cucumis melo (CUCME, Muskmelon) (Sebastian et al., 2010; Garcia-Mas et al., 2012), Schrenkiella parvula (SCHPA (formerly EUTPR)) (Dassanayake et al., 2011), Eutrema salsugineum (EUTSA (formerly THEHA)) (Oh et al., 2010), Brassica rapa (BRARA) (Wang et al., 2011), two Arabidopsis species, Arabidopsis lyrata (ARALY) and Arabidopsis thaliana (ARATH, thale cress) (The Arabidopsis Genome Initiative, 2000; Hu et al., 2011), Carica papaya (CARPA, papaya) (Ming et al., 2008), Gossypium raimondii (GOSRA, cotton) (Wang et al., 2012) and Theobroma cacao (THECC, cacao tree) (Argout et al., 2011). Note that the genus Thellungiella is now known as Eutrema, so the species formerly known as Thellungiella halophila (THEHA) is now known as Eutrema salsugineum (EUTSA). In addition, the species of Eutrema sequenced at the JGI has been determined to be salsugineum (a close relative of halophila). Therefore, this genome is actually classified as Eutrema salsugineum. Since we started our analysis before this change, all sequences related to the Eutrema salsugineum genome are annotated “THEHA” or “EUTSA” in our paper. Also, in the course of our study, the name Eutrema parvula (EUTPR) has been changed for Schrenkiella parvula (SCHPA). So, all Schrenkiella parvula (SCHPA) sequences are annotated “EUTPR” or “SCHPA”. See Supplementary Table S1 for bibliography and web links for download. Species tree representation (Supplementary Material) is based on several studies (Koch et al., 2000; Oh and Potter, 2005; Rensing et al., 2008; Forest and Chase, 2009; Magallcdn and Castillo, 2009; Wang et al., 2009; Oh et al., 2010; Sebastian et al., 2010; Arakaki et al., 2011; Reineke et al., 2011; Wang et al., 2011, 2012; Young et al., 2011; Tomato Genome Consortium, 2012; Varshney et al., 2012; Zhang et al., 2012).

On each of the 33 plant proteomes, the hmmsearch program was run to extract peptide sequences containing both LRR(s) and a KD (Eddy, 2009). Sequences containing both LRRs and KD were classified into SGs using a global phylogenetic analysis (Fischer et al., 2016). First, the KD of all these sequences was aligned using MAFFT with a progressive strategy (Katoh et al., 2002). Then the alignment was cleaned with TrimAl configured to remove every sites with more than 20% of gaps or with a similarity score lower than 0.001 (Capella-Gutierrez et al., 2009). A similarity matrix was computed using ProtDist with a JTT model, and then a global distance phylogeny was inferred using FastME configured with default settings and SPR movements to optimize the tree topology (Felsenstein, 1989; Desper and Gascuel, 2002). SGs were defined manually in the global phylogeny using the Arabidopsis genes as reference (Shiu et al., 2004; Lehti-Shiu et al., 2009; Fischer et al., 2016). To extend this dataset to receptor kinases devoid of LRRs in their ECD (sequences annotated “No_LRR”), the BLASTP algorithm (default parameters) has been run SG per SG, using each of the 7,767 KD sequences to search a database composed of the 33 proteomes as query (Altschul et al., 1997). Blast outputs were parsed to keep only homolog sequences sharing more than 90% identity with the query sequence. The new “No_LRR” sequences retrieved by blast were assigned to the same SG as the query sequence. Then, phylogenies were inferred for each of the 20 SGs. Each group of sequences was aligned using MAFFT with an iterative strategy (maximum of 100 iterations) (Katoh et al., 2002). Alignments were cleaned using TrimAl configured this time to remove sites with more than 80% of gaps (Capella-Gutierrez et al., 2009). Then maximum likelihood phylogenies were inferred using PhyML 3.0, configured with LG+gamma model, and the best of NNI and SPR topology optimization (Guindon and Gascuel, 2003). Statistical branch supports were computed using the aLRT/SH-like strategy (Guindon et al., 2010; Anisimova et al., 2011). Each of the 20 phylogenetic trees has been reconciled with the species tree using RAP-Green (Dufayard et al., 2005)¹. By comparing the gene tree with the species tree, this analysis allows us to root phylogenetic trees (Dufayard et al., 2005). We tested this approach of rooting (by minimizing the number of inferred duplications and losses) and compared it with rooting with outgroups (data not shown). The two methods provided very close root locations that did not change the overall conclusions.

To define the monocots dicots (MD) OGs, monocots/dicots bifurcations (branch support threshold >0.85) have been manually located in each of the 20 SG-specific trees. To be considered as MD OGs, the minimum number of monocots and dicots species represented was 3 and 4, respectively, to avoid keeping groups of misannotated sequences as MD OGs. Thus, 101 MD OGs have been defined. For each of them, the number of sequences in each of the 31 studied angiosperm species was recorded. If no sequence was discovered in one species, it was considered lost in this species.

Number of motifs and positions (LRRs (PF00560.24) and KDs (PF00069.16)) are outputs of the hmmsearch program (default parameters) (Eddy, 2009). Island domains were determined based on predicted LRR positions in sequences. For SP and TM domains, the TMHMM and TOPPRED softwares have been used (Claros and von Heijne, 1994; Krogh et al., 2001). For the malectin domains in SG_I and SG_VIII-2, sequences of the domains were extracted in SMART and aligned to build hmm motifs with the hmmbuild program (Eddy, 2009; Letunic et al., 2009). For Cys-pairs, hmm motifs were built based on subsets of sequences known to possess these motifs (Eddy, 2009).

Twelve sub-trees were considered: the OGs selected were those with a ‘simple’ organization, i.e. with a gene topology fitting approximately with the species tree. Three to four sequences among those that are the most closely related to the OG were selected from the whole SG tree as outgroup (Supplementary Figure S4). The sequences were re-aligned and the alignment was cleaned as described previously (Fischer et al., 2014, 2016). We ran codeml branch/site models implemented in the PAML4 software (Yang, 2007). For each OG, the following branch partition was defined: all branches but one were tagged as ‘background’ branches and the branch between the duplication node and the node corresponding to the split between monocot and dicot tagged as ‘foreground’ branch. Then two models were compared: the null model (A₀), in which sites on the fore- and background branches evolved under the same selective pressure (purifying or neutral), and a model including positive selection (model A) in which some sites on the foreground branch evolved under positive selection whereas sites on the background branches still evolved under purifying selection or neutrality. The most likely model was inferred by a likelihood ratio test (LRT). To take into account multiple testing, a Bonferroni correction was applied: the significance threshold of 0.05 was divided by the number of tested branches (24). Sites detected to be under positive selection at the codon level were manually curated for alignment quality and reliability. In branches identified to have evolved under positive selection, Bayes empirical Bayes was used to calculate the posterior probabilities at each codon and detect those under positive selection (i.e., those with a posterior probability of ω > 1 strictly above 95%).

We conducted a phylogenetic analysis of the LRR-RLK gene family in 33 fully sequenced plant genomes to classify them into SGs and to highlight and describe general characteristics of these LRR-RLK gene sets (Fischer et al., 2016). Briefly, besides 31 angiosperms genomes – represented by eight monocots (including six poaceae) and 23 dicots – one bryophyte genome of Physcomitrella patens (PHYPA, moss) and one lycopodiopsida genome of Selaginella moellendorffii (SELML, spikemoss) were included (Supplementary Table S1, see Section “Materials and Methods” for details and five-digit species code). As it has been done previously, we based our classification of LRR-RLK genes into SGs on the KD phylogeny (Shiu et al., 2004; Lehti-Shiu et al., 2009; Fischer et al., 2016). In our previous study (Fischer et al., 2016), the LRR-RLK dataset contained 7,767 sequences possessing at least one LRR in their ECD. Since we focus on structural features of ECDs and presence/absence of the genes in LRR-RLK OGs in this present study, we included LRR-RLK homologues for which LRRs were completely lost or were degenerated in this new dataset. Thus, 593 sequences (prefixed “No_LRR”) were added to the original set of 7,767 sequences which lead to a total of 8,360 sequences (Supplementary Table S2). Within each of the 20 SGs, KDs were aligned and SG-specific trees were obtained using a likelihood-based method (PHYML) (Figure 1A and Supplementary Material). Altogether, these LRR-RLK genes represent on average 0.71 and 0.66% of the monocot and dicot proteomes, respectively. Interestingly, in moss (PHYPA) and spikemoss (SELML), the proportions of LRR-RLKs per genome (0.41 and 0.36%, respectively) are approximately half the ratio observed in angiosperms. Likewise, the average number of LRR-RLK genes in angiosperms is 263.6, with 260.7 LRR-RLK proteins [±20.2 (SE)] in dicots and 268.8 [±18.2 (SE)] in monocots. In PHYPA and SELML, 134 and 81 LRR-RLK genes have been retrieved, respectively. There is no significant difference in the average number of LRR-RLK genes between monocots and dicots but the number almost doubled in most angiosperms compared to PHYPA and SELML. It has to be noted that in some genomes (e.g., CARPA and LOTJA), the number of LRR-RLK genes is particularly low compared to other angiosperm genomes, suggesting that retention rates vary among genomes, and that many losses may have occured in some genomes (Fischer et al., 2016). Nevertheless, our results highlight the fact that, after the first wave of expansion in early land plants (Embryophyta), a second large amplification occured in angiosperm genomes which shaped the current LRR-RLK family size of more than 200 gene copies on average per genome.

Among the 20 SGs, SG_III, SG_XI, and SG_XIIa are the largest as they contain ∼50% of the total number of LRR-RLK genes in the analyzed plant genomes (Figure 1A). The extensive expansions leading to their size do not follow the same amplification pattern (Figure 1B). This observation is highlighted in Figure 1B by the color code used for each species in the 20 SG-specific trees (with branches of monocots species in pink and red, and branches of dicots species in yellow, blue and green). Our results reveal that the high numbers of SG_XIIa genes is the consequence of many lineage-specific expansions (LSE) (See Fischer et al., 2016 for details). These LSEs are relatively recent as they can be observed in phyla as well as species-specific lineages. On the contrary, in SG_III and SG_XI, expansions occured mainly before the early divergence of angiosperm lineages, even though LSEs can also be observed at different levels of resolution in the trees. Therefore, these numerous and diverse modes of expansions lead to complex paralogous and orthologous relationships.

With the aim of transfering functional annotation from well studied genes from model species to orthologous genes in other genomes, we first analyzed in depth orthologous relationships between monocots and dicots LRR-RLK genes in each SG. This analysis led us to define what we named the “core set” of LRR-RLK genes in angiosperms: i.e., orthologous genes which have not been completely lost in either monocots or dicots throughout the angiosperms evolutionary history. To do so, the 20 SG-specific trees were scaned to locate monocots/dicots bifurcations (Figure 2). Based on this analysis, 101 OGs containing monocots and dicots sequences (named MD OGs) were characterized and defined as the “core” set of LRR-RLK genes in angiosperms.

The SG analysis revealed that these 101 MD OGs are present in 19 of the 20 SGs, with the majority of them in SG_III and SG_XI (Table 1). This highlights again the fact that most of the LRR-RLK genes which underwent expansions before the monocots/dicots split have been retained in these SGs. In order to go further in the description of orthologous relationships, we qualified OGs as either “simple” or “complex”. In “simple” OGs, the presence or absence of duplications within the monocot or dicot clades can clearly be inferred from the phylogenetic tree. On the other hand, if several duplications occurred disorderly with no obvious connection to the species tree we described these OGs as “complex” (OG_c). Interestingly, these OG_c are over-represented in SG_IV, VIII-1, VIII-2, and XI. A total of 6739 genes are contained in the 101 OGs, representing 82.7% of the entire LRR-RLK gene family. However, while 2956 genes are included in the 24 OG_c (average of 123.2 genes per OG), 3783 genes belong to the 77 non-complex OGs (average of 49.2 genes per OG), highlighting differences in expansion/retention rates between these OGs. Moreover, looking at the percentage of genes contained in OGs per SG reveals that, except for SG_VIII-2 and XIIb, more than 70% of the LRR-RLK genes belong to OGs, with genes mainly in complex OGs in SG_I, IV, VIII-1, VIII-2, XI, and XIIa (Figure 3A). In some SGs, some OGs contain a very large number of genes, such as 405 genes in one of SG_I OG (SG_I-3c), or 307 and 708 genes in two SG_XIIa OGs (Figure 3B and Supplementary Table S2 for details). In these large OGs, many species-specific duplications occured, and among the 10 OGs containing more than 100 genes, 8 are complex. In SG_I-3c for example, several genes have been studied in Arabidopsis, such as IMPAIRED OOMYCETE SUSCEPTIBILITY 1 (IOS1), FLG22-INDUCED RECEPTOR-LIKE KINASE 1 (FRK1), and light-repressible receptor protein kinase (LRRPK) (Deeken and Kaldenhoff, 1997; Asai et al., 2002). All have been reported to be involved in abiotic and biotic responses in dicots but no gene from the same OG have been described so far in monocots. However, the fact that these genes are classified into the complex mode of expansion suggests that in monocots too, these genes could be involved in stress response.

For each OG, we investigated whether some species were lacking members, focusing particularly on the Brassicales, for which six species are included in our analysis (Figure 4). Moreover, this clade contains the model plant Arabidopsis which is the reference for many studies on LRR-RLK functions (Wu et al., 2016). Interestingly, among the 101 OGs, 14 (13.8%) have been completely lost in the Brassicaceae, and 3 of them are even absent in all the Brassicales. This observation is an incentive to the extensive study of these receptors in other plants than Arabidopsis, adding an argument to the fact that functions or interactions are sometimes phylum-specific. For example, the SG_I-2 OG contains the SYMBIOSIS RECEPTOR LIKE KINASE (SYMRK, also known as NORK or DMI2) receptor which is involved in actinorhizal plants and legumes, respectively, in phosphate-acquiring arbuscular mycorrhiza and nitrogen-fixing root nodule symbiosis (Antolin-Llovera et al., 2014a). For this gene, we noticed that besides being absent in the Brassicaceae, which do not form mycorrhizal associations and root nodule symbiosis with rhizobia, some other characteristics of these receptors have been observed in monocots (see below).

The use of orthologous relationships to infer functional annotations relies on the fact that orthologs are expected to carry equivalent functions in different organisms. However, this can only be reliably inferred if, at least, structural characteristics and domain architecture are conserved between othologs. In our analysis, the phylogeny of the LRR-RLK proteins was computed on the well conserved KDs. However, LRR-RLK sequences are composed of several domains and, notably, of LRRs in their ECD. One could wonder whether the domains belonging to genes of the same OG are conserved. First, we took a detailed look at the predicted number of LRRs of all these receptors. The number of LRR motifs per protein is an important feature for homo- and hetero-complex formation between LRR-RLKs (Macho and Zipfel, 2014). Second, we investigated the presence of island domains in between LRRs. These domains have been described to be the binding site for the BR hormone in some receptors (Kinoshita et al., 2005; Hothorn et al., 2011; She et al., 2011). Third, we analyzed the presence of the MLD, a carbohydrate-binding domain, and the GDPC, a protein cleavage motif, which were shown to be located before the LRRs in some SG_I receptors. Fourth, we investigated the presence of Cys-pairs surrounding the LRRs in some SGs. The presence and organization of these domains is functionally important and has to be taken into account for transfering functional informations between orthologous genes. The description of structural features localized in the ECDs of these receptors allows a subclassification which, although reflecting the phylogeny of the KD of these receptors, also takes the structural differences of the ECD into account. Therefore, we subdivided the 20 SGs further according to these characteristics in their ECD (Figure 5 and Supplementary Table S3 for details).

Twelve pairs of MD OGs present in almost all monocot and dicot species, and harboring a gene topology fitting approximately the species tree, appear to be issued from ancestral duplications predating the monocot/dicot divergence (Table 2). As these OG pairs had a similar ECD structural organization and were kept in almost all species studied here, we searched for potential positive selection footprints in the divergence leading to their differentiation. Although these genes do not all have a known function, it is expected that they underwent amino-acid changes leading to their sub- and/or neo-functionalization, and one can wonder if and how positive selection could have driven these changes.

To answer this question, we tested whether some sites underwent positive selection on the two branches starting from the ancestral duplication and ending at the monocot/dicot divergence node of each OG. The detailed results of this analysis are presented in Table 2 and Supplementary Table S5. Two pairs showed no signal on either of the two branches (SG_Xa-1/2, SG_XIIa-2/3). For the pairs SG_II-3/4 and SG_XIIIa-1/2, a signal was detected for one branch only but the signal on SG_XIIIa-1/2 may be a false positive or the sign of a lack of power, since no sites appeared to be significant (see Materials and Methods for details). The eight other pairs showed a signal of positive selection on each branch. Although the model indicating positive selection performs significantly better than the null model, two pairs (SG_III-3/4 and SG_XI-19/20) have no significant sites for one of the two tested branches. This again indicates either a false positive or a lack of power. It is also possible that positive selection acted on a large number of sites which results in none of them exceeding the significance threshold. Finally, five pairs (SG_III-8/9, SG_XI-8/9, SG_XI-14/15, SG_XI-17/18, and SG_XIIIb-1/2) have a strong signal with up to 26 sites validated after manual curation. This result shows that in about half of the tested cases, several amino acid changes fixed in the divergence between these genes are compatible with a signal of positive selection.

A total of 141 sites were manually validated as having experienced an episode of positive selection during MD OG genes divergence. For the five pairs with a strong signal, the repartition of these sites across the different domains of the LRR-RLK protein showed that the LRRs and KDs are the most affected (Figure 7). More than half of the sites (78) fall in the ECD, among which 68 are in the LRR domain; 51 sites fall in the ICD, among which 42 are in the KD (Table 2). Considering that these LRR and KD are the largest domains, we normalized the number of positively selected sites by domain size. Kinase and LRR appeared then to be affected equally (Chi-square test, p = 0.26) by positive selection. This result is very different from positive selection signatures observed in the recent paralogs, for which LRR is the most strongly affected domain (Fischer et al., 2016). The number of sites laying in the LRR domain allowed us to look for any specific distribution across the 24 amino acids composing the motif. The repartition of the sites affected by positive selection within the LRR is not homogeneous and the majority of them (67) fall in the 13 non-canonical positions (Supplementary Table S5). However no notable pattern emerges from their distribution (data not shown). Again, this result contrasts with what is observed in lineage-specific expanded genes for which four positions are predominantly affected (Fischer et al., 2016). The remaining domains are affected by a number of sites varying from 1 to 9.

This approach revealed a prevalence of sites targeted by positive selection in the ECD for the three couples of genes belonging to SG_XI as well as SG_XIIIb. On the opposite, a tendency to target ICDs can be observed for SG_III. Two pairs of OGs for which positive selection footprints are detected in the ECD correspond to OGs whose duplication gave rise to the PXY and PXLs clades in SG_XI-17/18 (Fisher and Turner, 2007; Jung et al., 2015); and to the ERECTA (ER) and ERECTA-like (ERL) clades in SG_XIIIb-1/2 (Torii et al., 1996; Sanchez-Rodriguez et al., 2009). Other pairs of OGs concern differentiation of the PEPR1 and 2 clade (SG_XI-9), of the STERILITY-REGULATING KINASE MEMBER 2 (SKM2) gene (SG_XI-14) or of the DspA/E-interacting protein of Malus x domestica Borkh 1 and 3 (DIPM1 and 3) clade (SG_III-9), from their respective sister clades, SG_XI-8, SG_XI-15, and SG_III-8 (Meng et al., 2006; Krol et al., 2010; Kang and Hardtke, 2016). In these clades, no gene has been described yet. The strong signal of positive selection detected for these five groups of genes indicates that the divergence between ancestral copies may have procured a selective advantageous: in the domain involved in ligands or partners binding for ECDs, or in the domain affecting downstream signaling pathways for ICDs. Indeed, during the early expansion of LRR-RLK that took place before angiosperm split, some duplicated LRR-RLK differentiated by fixation of a higher number of non-synonymous than synonymous mutations at some amino acid sites, indicating the emergence of probably new advantageous functions.