Conserved motifs in RALF proteins (Figure S1) were identified using MEME (Bailey and Elkan, 1994) using 37 previously identified RALFS as a query (Cao and Shi, 2012). These motifs were identified in other plant species by scanning with FIMO (Grant et al., 2011). Proteomes for the following species were retrieved from Phytozome v11.0 (Goodstein et al., 2012): Amaranthus hypochondriacus, Amborella trichopoda, Ananas comosus, Aquilegia coerulea, Arabidopsis halleri, Arabidopsis lyrata, Arabidopsis thaliana, Brachypodium distachyon, Brachypodium stacei, Brassica rapa, Capsella grandiflora, Capsella rubella, Carica papaya, Chlamydomonas reinhardtii, Citrus clementina, Citrus sinensis, Coccomyxa subellipsoidea C-169, Cucumis sativus, Eucalyptus grandis, Eutrema salsugineum, Fragaria vesca, Glycine max, Gossypium raimondii, Linum usitatissimum, Malus domestica, Medicago truncatula, Micromonas pusilla, Mimulus guttatus, Musa acuminate, Oryza sativa, Ostreococcus lucimarinus, Panicum hallii, Panicum virgatum, Phaseolus vulgaris, Physcomitrella patens, Populus trichocarpa, Prunus persica, Ricinus communis, Salix purpurea, Selaginella moellendorffii, Setaria italic, Setaria viridis, Solanum lycopersicum, Solanum tuberosum, Sorghum bicolor, Spirodela polyrhiza, Theobroma cacao, Vitis vinifera, Volvox carteri, Zea mays, and Zostera marina. Genomic data such as genome size and total gene number were extracted from the relevant publication for each species (see Table S1). A stringent q-value cut off of 0.05 was used to remove low-scoring hits to the six motifs identified by MEME.

All protein alignments were initially created using the MUSCLE algorithm (Edgar, 2004) within the AliView alignment editor (Larsson, 2014). This was followed by manual optimisation to improve the alignment in regions that had been clearly misaligned by MUSCLE. Approximately-maximum likelihood phylogenetic trees were created for full-length proteins and mature peptides using FastTree v2.1 (Price et al., 2010), with 4 rounds of minimum-evolution SPR moves (option: -spr 4) and exhaustive nearest-neighbor interchanges (options: -mlacc 2 -slownni) to improve accuracy. All other parameters were left as default, including the calculation of local-support values by the Shimodaira-Hasegawa test (Shimodaira and Hasegawa, 1999). The inferred phyloXML trees were viewed in Archaeopteryx v0.9916 (Han and Zmasek, 2009) to identify RALF clades and sub-clades. WebLogo3 (Crooks et al., 2004) (http://weblogo.threeplusone.com/) was used to provide a visual summary of conserved residues within the alignments. The inferred phylogenetic trees have been uploaded to the TreeBASE repository (Piel et al., 2002) and can be accessed at http://purl.org/phylo/treebase/phylows/study/TB2:S20366.

To perform pairwise BLAST (Altschul et al., 1990) hits between each individual protein, full-length unaligned sequences of all 795 RALFs were uploaded to the online CLANS server (Frickey and Lupas, 2004) that is a part of the MPI Bioinformatics Toolkit (Alva et al., 2016). The BLOSUM80 substitution matrix was used for scoring and the HSP option was enabled. The output from this server was used within the standalone CLANS software to create the 2D similarity plots. A p-value threshold of 1e⁻¹⁰ was used to remove low-scoring BLAST searches and a minimum attraction value of 50 was used to restrict the proteins to within a reasonable 2D space. As CLANS is non-deterministic, we ran the analysis many times, each from different initialization states, to confirm that the separation of the proteins was reproducible and consistent. We stopped each analysis when the proteins had become stationary within the 2D plot.

The online multi-cleverMachine tool (Klus et al., 2014) was used to compare the physico-chemical properties of RALF proteins. Unaligned sequences for each clade were uploaded to the server in FASTA format, with “clade IV” proteins being denoted as the positive set and clade I, II, and III proteins as negative sets. Ten scales were used to predict each property and detailed statistical comparisons were viewed using the boxplot and ROC curve functions found within the output screen.

The mRNA expression of A. thaliana and Z. mays RALF genes was analyzed across publically-available RNAseq datasets that are included within Genevestigator (Hruz et al., 2008). A total of 1031 samples across all datasets were analyzed. A clustered heat-map representing the log-2 absolute expression values throughout the anatomy of the plant was obtained using the “Hierarchical Clustering” tool within Genevestigator, with both the genes and conditions subjected to Euclidian-distance clustering.

Previous genome-wide identifications and phylogenetic analyses of the RALF family have focused on either six (Cao and Shi, 2012) or four (Sharma et al., 2016) plant species. We sought to broaden this range by taking advantage of the Phytozome v11.0 genomic resource (Goodstein et al., 2012), which provides comprehensive sequence data and accompanying annotation for a variety of green plant species. A full list of the 51 species included in this study can be found in Table 1. This diverse range of species includes lineages that have been excluded from previous RALF phylogenetic analyses, such as the Rosaceae, and generally allows for a more informative analysis with greater resolution. Additionally, though a small number of RALFs have been previously found in plants such as potato (Germain et al., 2005) and M. truncatula (Pearce et al., 2001; Combier et al., 2008), these studies were performed before a full genome was available, and were instead reliant upon expressed sequence tags (ESTs). Therefore, revisiting such species in light of the numerous full plant genomes now available should allow for a more accurate representation of the RALF family.

To find RALF proteins across the 51 included species, MEME (Bailey and Elkan, 1994) was first used to detect up to six conserved amino acid motifs within the previously identified 37 A. thaliana RALFs (Cao and Shi, 2012; Figure S1). This was followed by the identification of regions matching to these motifs across the proteome of all 51 species using FIMO (Grant et al., 2011). We discovered a total of 795 RALFs, with a breakdown of the number per species shown in Table 1. None of the analyzed proteomes were found to contain more than the 37 RALFs identified in A. thaliana, though there are other eudicot species with large numbers of RALFs, such as the 33 found in the apple (M. domestica) proteome. In general, the eudicots contain more RALF proteins on average (~20) than the monocots (~14), however, the angiosperm with the fewest number of RALFs is the early-diverging rosid domesticated grape (V. vinifera), which only has four RALF proteins.

The difference in the average number of RALF genes between the monocots and eudicots could mean that the genetic mechanisms underlying the evolution of the RALF family was distinct in each group. On the other hand, such differences may simply be a consequence of the eudicots analyzed having larger genomes than the monocots. To investigate this, for each species we compared the number of RALF genes to the total number of genes in the genome (Figure 1), genome size (Mbp; Figure S2), and gene density (genes/Mbp). We found a strong positive correlation (r = 0.66) between the number of RALFs and genome size for monocots, but a weak negative correlation (r = −0.13) for eudicots. Furthermore, the number of RALF genes correlates very strongly (r = 0.93) with the total number of genes across the monocots, but there is only a weak correlation (r = 0.32) for the eudicots. This data suggests that in monocots, RALF diversification has occurred at a very consistent rate that is proportional to overall changes in genome size, such as those caused by genome duplications. In eudicots however, some species appear to have far more RALFs than can be explained by expansion of the genome alone. These have been circled on Figure 1 and Figure S2. Interestingly, we noted that five of these species belong to the Brassicaceae. When the Brassicaceae are omitted from this analysis, the number of RALF genes correlates more strongly with the genome size (r = 0.31) and gene number (r = 0.69) for eudicots, coefficients that are much higher than before but still below those of the monocots.

In order to understand the evolution of the RALF family in more detail we aligned the 795 identified RALF protein sequences and inferred a phylogenetic tree The tree separates into four major clades (Figure 2) that consistently shows very high support values (>0.8; Figure 2A) Clade III is the largest of the four major clades with 320 members, followed by clade IV with 264, clade II with 151 and clade I with only 49 proteins. All four clades contain RALFs from a variety of species, including both monocots and eudicots (Table 1), suggesting that all clades evolved before the divergence of these two angiosperm lineages. However, 90% of the genes within clades I and II are from eudicot species, with the Poaceae (grasses) being entirely absent from these clades. Conversely, there is a notable overrepresentation of monocots within clade III, as 70% of the monocot RALFs are found here.

Nine RALFs were identified within the A. trichopoda proteome, the most basal angiosperm (Zuccolo et al., 2011). Only 4 of the 43 angiosperm species studied here contain fewer than nine RALFs, meaning that there has been a general diversification, rather than a contraction, of the RALF family within almost all lineages since the early beginnings of the angiosperms. In support of this, A. trichopoda contains RALF proteins belonging to clades II, III and IV, explaining why these clades are represented across both the monocots and eudicots. The early-diverging species P. patens and S. moellendorffii have fewer RALF genes than A. trichopoda and these are instead found within clade III. We could identify no RALF genes within the five chlorophyte species (Table 1).

Clades II, III, and IV can be further split into distinct sub-clades, with strong local support (Figure 2). Each of the nine sub-clades contains a range of species from across the angiosperms (Table 2), suggesting that these sub-clades evolved within the ancient angiosperms. This is evidenced by the spread of the 9 A. trichopoda RALFs across 6 of the 10 sub-clades. However, clade III(D) is almost entirely absent from the eudicots, but in monocots has expanded to become the most prevalent clade. How each of these subclades differs in terms of amino acid sequence will be considered below.

It has been established that using a large number of aligned sequences does not necessarily lead to an accurate phylogenetic tree (Philippe et al., 2011). In order to assess the accuracy of our tree and to further investigate the relationship between the various clades, we sought an alternate method of visualizing protein relationships. CLANS performs all-against-all BLAST searches upon an unaligned dataset and returns a non-deterministic 2D or 3D map in which protein similarity can be visualized (Frickey and Lupas, 2004). Figure 3 shows a typical CLANS output of the 795 full-length RALF sequences, in which each protein is represented by a colored dot. The color of each dot denotes the sub-clade that the protein fell into within the phylogenetic tree (Figure 2B) and the distance between the proteins signifies their similarity. Of note is the high level of correspondence between CLANS and the phylogenetic tree, demonstrated by scarce intermixing of proteins from different major clades within the 2D space. Hence, CLANS independently finds that RALFs generally have more sequence similarity with other proteins from the same clade than they do with members of other clades. This confirms that the major splits of the phylogenetic tree (Figure 2B) represent detectable differences in the underlying amino acid sequences. Clade IV RALFs appear to be the most diverged, as they are mostly restricted to the periphery of the CLANS map. Clade I, II and III instead fall closely together, with the central mix of clade I and II proteins suggesting that some RALFs within these clades have very similar sequences.

Furthermore, RALFs consistently cluster within the CLANS plot (Figure 3) into their respective sub-clades, signifying that these sub-clades also represent distinguishable variations in sequence. Of interest is the separation of clades II(A) and II(B), with II(A) intermixing with clade I whilst II(B) clusters slightly further away. This is in contrast with the inferred phylogenetic tree and suggests that clades I and II(A) share more sequence similarity than either does with clade II(B).

We questioned whether the separation of the 795 identified RALFs into four clear clades was due to variations within the mature, functional peptide sequence or because of residues outside of this region that are perhaps subject to less selection pressure. We removed the N-terminus region from the aligned sequences, leaving only the YISY motif and all downstream residues. Although the beginning of the peptide is thought to be a few residues upstream of the YISY motif (Matos et al., 2008), variation within this region makes it difficult to accurately identify the start of the peptide and hence these residues were omitted for simplicity. A CLANS analysis of this peptide region (Figure 4A) produces a very similar distribution to the full-length preproprotein CLANS, suggesting that residues within the mature peptide have sufficiently diverged alongside those outside of this region. Of note, the RALFs found within clade IV on the full-length tree again fall to the periphery of the peptide-only CLANS plot, demonstrating substantial differences within their peptide sequence compared to the other clades. Conversely, there is a tight cluster of peptides from the full-length clades I, II, and III, highlighting the similarities of these peptides. In light of this, we aligned the mature peptide region from clades I, II, and III and created an additional approximate-maximum likelihood tree to see whether the peptides from these three clades could be distinguished. The tree (Figure 4B) reliably separated the clade III mature peptides from those in clades I and II, which appear to be mostly indistinguishable. The clear overlap between the phylogenetic analyses of the full-length preproprotein and mature peptide suggests that there have been divergences across the whole length of the protein, including the functional peptide region.

As both methods of assessing protein similarity broadly agreed that the RALF family has diverged into distinct groups, we carried out a more detailed analysis of the underlying protein sequences. We aligned individual clades in order to identify any distinguishing characteristics of each. Inspection of the underlying sequences reveals that the proteins of clades I and II(A) are indeed very similar. As shown in Figure 5A, the consensus mature peptide sequence of clade II(B) is much distinct from clades I and II(A), lacking a conserved YYNC motif whilst containing additional proline residues toward the C-terminus. Such differences are consistent with the relative placement of these three clade/sub-clades within the CLANS output (Figure 3).

The RALFs that occupy clade III are also very similar to those of clade I, in agreement with their close proximity in Figure 3. However, these proteins appear to have diversified somewhat from the remarkably well-conserved RALFs that belong to clades I and II(A). There are no clear characteristics that can distinguish clade III as a whole, demonstrated by the amount of variation at many residue positions of the mature peptide (Figure 5A). Instead, each of the four subclades has seemingly diversified differently, though there is noticeable diversification even within each subclade. The most distinguishing variations are found within the mature peptide region, as can be seen in Figure S3. Obvious examples include the insertion of a proline residue at position 29 within clade III(D), and the insertion of an alanine at position nine of clade III(C) RALFs. Additionally, the CRG motif that occupies the three terminal residues of almost all clade III(D) peptides is entirely unique to that sub-clade. Notably, clades III(A/B) commonly contain an additional di-basic site upstream of the YISY motif that is mostly absent from clade III(C) and entirely absent from clade III(D). It is possible that the presence of this di-basic site within close proximity to the di-basic RALF cleavage site (RRIL; Matos et al., 2008) could affect the processing of these preproproteins.

Whereas clades I, II, and III all possess the conserved YISY motif that is thought to be responsible for the binding of the peptide to its receptor (Pearce et al., 2010), remarkably this motif is rarely found within clade IV. Only 3/264 clade IV RALFs contain “YISY,” with the remainder showing a diverse range of substitutions, though many still contain an isoleucine and tyrosine at the second and fourth positions (XIXY). This lack of YISY conservation can be visualized in Figure 5A. Furthermore, many other typically-conserved residues are absent from this clade. The RRXL protease cleavage site, found upstream of YISY within the vast majority of clade I, II, and III RALFs, is almost entirely absent across the 264 clade IV proteins. This likely means that these RALFs are processed and cleaved through a different mechanism, if at all. Likewise, the acidic (glutamate/aspartate) region usually found between the signal peptide and the mature peptide is missing, with this probably impacting substantially upon the protein's structure and stability. In fact, only a minority of peptide residues are conserved within clade IV, with most residues being extremely variable (Figure 5A). The absence of these motifs and other residues results in the clade IV RALFs having a mean length of only 88 amino acids, in contrast to the other clades (Figure 5B), which contain RALF proteins with an average length of 125 amino acids. The missing regions likely explain the position of clade IV RALFs at the periphery of the CLANS 2D plot, away from the central zone occupied by the other three clades (Figure 3). The higher frequency of RALF proteins within the Brassicaceae is specifically due to an overrepresentation of clade IV RALFs, with this clade representing 56% of the RALF proteins within the Brassicaceae species analyzed, in comparison to the average representation of 34% across all 51 species.

We questioned whether the distinctive sequence patterns of the clade IV RALFs are likely to have any significant impact upon the physico-chemical properties of the translated protein. CleverMachine (Klus et al., 2014) allows for the detection of protein properties that differ between two datasets and has been used to distinguish P-bodies and stress granules from other globular proteins (Marchese et al., 2016) and to classify homo-repeat proteins (Yu Lobanov et al., 2016). A multi-cleverMachine property prediction and comparison of the RALFs reveals that the clade IV proteins differ in a variety of physico-chemical properties from the other clades (Figure S4), such as a reduced disorder propensity. The analysis revealed that most properties could distinguish the clade IV proteins with high accuracy, as demonstrated by the typical area under the ROC curves being >0.9, a score typically considered to indicate a highly accurate test (Greiner et al., 2000).

To identify whether genes within clades might have a common function and hence share similar expression profiles, we analyzed the expression of each A. thaliana RALF gene across various publically-available RNAseq datasets using Genevestigator (Hruz et al., 2008; Figure 6). We found that the expression of clade IV RALFs is almost exclusively restricted to inflorescence tissues, with only 1 of the 23 A. thaliana clade IV genes (AT4G14020) showing notable levels of expression within other anatomical regions such as the root. In contrast, genes from other clades exhibit a more widespread expression profile and expression is found within the root and shoot as well as flowers, with the specific exception of the subclade III(B) which has a very similar expression pattern to the clade IV genes. This data further suggests that there have been functional diversifications between clades and genes within a clade share expression patterns. We found similar expression patterns within Z. mays (Figure S5), with the clade IV genes again being restricted to the inflorescence tissues.

Whereas most residues within the mature peptide are highly conserved, there exists a great deal of variation at the C-terminus. Although the majority of RALFs terminate with an RCRR motif, many have additional residues downstream. The composition of these residues is highly variable and each is usually restricted to a few closely related species, suggesting that these are relatively recent additions to the protein. A selection of these is shown in Figure S6 to demonstrate their variability. The longest of these are found within L. usitatissimum, which contains two RALF proteins with lysine-rich C-terminal tails that are over 250 residues in length. We also identified a F. vesca gene (gene10567-v1.0-hybrid) which is a fusion of an N-terminus LOW PSII ACCUMULATION1 (LPA1) homolog to a C-terminus RALF protein. In Arabidopsis, LPA1 is known to be involved in the assembly of photosystem II (Peng et al., 2006). Finally, whereas some CLAVATA3/ESR-related (CLE) and C-TERMINALLY ENCODED PEPTIDE (CEP) genes containing multiple peptide motifs have been reported (Oelkers et al., 2008; Sawa et al., 2008; Roberts et al., 2013), we could find no evidence of RALF genes that contain more than one RALF peptide motif.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.