The evolutionary history of plant T2/S-type ribonucleases

T2/S-RNase sequences and alignment

We obtained known T2/S-RNase amino acid sequences from the protein database of GenBank release 202 (see Section S1 for query strings). Groups of sequences that shared 90% or higher sequence identity were identified using UCLUST clustering algorithm in USEARCH version 7.0.1090, and only the longest sequence from each group was retained (Edgar, 2010). The resulting set was then used to query core nucleotide (NT) and expressed sequence tags (EST) databases with the tblastn algorithm in BLAST version 2.2.29, with default settings and an expected value cut-off of 1 × 10⁻¹⁰ (Altschul et al., 1990). The limited taxonomic search included only data from the Viridiplantae (green algae and land plants), and was restricted to sequences between 300 to 10,000 bp (see Section S2 for exact query). The BLAST results from EST and NT databases contained 3,679 and 2,380 unique accessions, respectively, from 411 species.

With the exception of Petunia × hybrida and Solanum lycopersicum, sequences from domesticated and hybrid species where excluded (n = 26, see Section S3 for a list of excluded species). The number of allelic S-RNase sequences was deliberately reduced for efficiency, and only the longest high-quality representative sequences were kept for all available genera. Coding regions from nuclear DNA and mRNA sequences from NT and EST databases were aligned using MAFFT version 7.158b (Katoh & Standley, 2013). A maximum likelihood guide tree was constructed using RAxML version 8.0.26 (Stamatakis, 2014). With the exclusion of known S-alleles, groups of monophyletic congeneric sequences that shared at least 98% sequence identity were identified using UCLUST. Sequences within each group were aligned using MAFFT version 7.158b and visually inspected in Geneious version 7 (created by Biomatters, available from http://www.geneious.com). We removed any present polyadenylation tails or ambiguous characters at either end of each sequence. Overlapping groups (unigenes) were collapsed to a majority rule consensus sequence with ambiguities introduced as necessary. The resulting set of sequences was reviewed, and only the sequences longer than 350 bp that contained at least three out of five conserved sequence motifs found in this gene family were kept. Catalytically active histidine residue (CAS II) is known to be essential to ribonucleic activity of T2-type RNases (Jost et al., 1991; Taylor et al., 1993; Irie, 1999). We identified and kept the proteins that were apparently missing this residue, as well as proteins established as catalytically inactive in functional studies (MacIntosh et al., 2010; Ohkama-Ohtsu et al., 2004; Wei et al. , 2006; Gausing, 2000; Van Damme et al., 2000; Kim et al., 2004).

In the guide tree constructed during the review process, sequences in our dataset formed two distinct clades. One of these clades consisted of T2/S-RNases, the other was composed of Thioredoxin-domain-2-containing disulphide isomerases. These sequences, excluded from analyses, were likely detected by BLAST search because they contain a sequence motif similar to T2/S-RNase conserved region 3. The processed dataset sourced from GenBank searches contained 618 sequences (see Table S1 for the list of GenBank accessions).

Retrieval and processing of sequences from genomes

We also used the above-generated dataset to query 146 available sequenced plant genomes (as of September 2014), using the blastn algorithm (see Table S2 for the list of genomes used). The gene structures of the BLAST hit results including the adjoining upstream and downstream 3 kb segments were annotated with Exonerate 2.2.0 (Slater & Birney, 2005) using the translated GenBank dataset as templates. Annotated sequences containing premature stop codons or ambiguous amino acid residues were removed. To construct the final T2/S-RNase sequence dataset for alignment, we extracted the coding regions obtained from the genomic sequences (see Table S1 for a list of genomic sequences), concatenated them, and then pooled them with the GenBank dataset. Groups of sequences that shared 95% sequence identity (or higher) were identified using UCLUST, and then further subdivided by genus. The longest sequence from each such group was retained.

Final sequence dataset and alignment

After the inclusion of sequences obtained from the genomes and subsequent processing, the final set consisted of 715 sequences annotated as T2/S-RNases, with detectable T2/S-RNase features, and/or grouping with T2/S-RNases in preliminary analyses. The bulk of these (711) represented land plants (embryophytes), the rest came from distantly related chlorophyte algae to be used as potential outgroups (one from Chlamydomonas reinhardtii and Bryopsis maxima, and two from Volvox carteri). No sequences from more closely related streptophyte algae were available. These sequences were translation-aligned using MAFFT version 7.164b, the alignment was reviewed, adjusted manually and mapped back to nucleotide sequences.

Phylogenetic analyses

Separate gene trees were inferred using nucleotide and amino acid substitution models using MrBayes version 3.2.2 (Ronquist et al., 2012). Both analyses consisted of four independent runs with one cold and seven heated chains. We implemented general time-reversible models with Gamma-distributed among-site rate variation, and four Gamma rate categories. The runs were allowed to complete 780 and 200 million generations for nucleotide and amino acid models, respectively. The trees and parameters were sampled every 1,000 generations for both runs. The temperature parameter was periodically adjusted throughout the runs to ensure that the acceptance rates of attempted swaps between the cold and the heated chains fell within the target window of 20%–60%. The proposal probabilities for different moves were tuned so the acceptance rates fell within the target window of 20%–70%. Parameter convergence was assessed using R package RWTY (Warren, Geneva & Lanfear, 2017). For both analyses, all parameters have converged within the first fifty million generations. The tree topologies took much longer to converge, as judged with treespace plots. Based on this information, the burn-in was set to 580 million generations for nucleotide trees and 150 million generations for amino acid trees. Both posterior tree sets were used to generate maximum credibility trees as well as consensus trees with minimum clade frequency threshold of 0.75 using the program SumTrees version 3.3.1 (Sukumaran & Holder, 2010). The posterior sets from both analyses were resampled (every four million generations for nucleotide trees and every one million generations for amino acid trees) to obtain a total of 100 trees which were used as starting trees for maximum likelihood inference with RAxML version 8.1.17 (Stamatakis, 2014). Two sets of analyses were performed by fitting a general time-reversible models of nucleotide and amino acid substitution with the CAT model of rate heterogeneity. Support values for the highest scoring RAxML trees were calculated from the respective MrBayes posterior set of trees using SumTrees.

Analysis of intron positions

In order to investigate T2/S-RNase intron/exon structure, genomic sequences with gene structure annotations were first translation-aligned using MAFFT version 7.164b in Geneious and the alignment was manually adjusted. Introns were treated as homologous across sequences if their starting positions overlapped within a seven nucleotide window in the alignment. Next, introns were classified by their phase (position within a codon). Phase zero introns occur before the first base, phase one and two introns interrupt a codon triplet after the first or second base, respectively. Aside from plants, intron positions were also identified in the T2-type RNase loci from algae (Volvox carteri, Bryopsis maxima), animals (Amphimedon queenslandica, Strongylocentrotus purpuratus, Hydra vulgaris, Homo sapiens), and (fungi Saccharomyces cerevisiae, Aspergillus oryzae).

Isoelectric point (pI) value calculations

Isoelectric point (pI) is the pH at which a molecule, on average, carries no net electric charge. The pI value for each sequence included in the phylogenetic analyses was calculated using methods described in Bjellqvist et al., (1993) and Bjellqvist et al., (1994) implemented in the ProteinAnalysis tool in Biopython 1.64 (Cock et al., 2009). Signal peptide sequences were not included in pI calculations.

Identification of putative SFB genes located near T2/S-RNases

F-box motif-containing genes empirically linked with SI function are co-located with S-RNases and are approximately 1 kb long (Lai et al., 2002; Ushijima et al., 2003; Sijacic et al., 2004; Entani et al., 2003). Most lack introns, although a single intron has been reported in untranslated upstream region of Prunus avium SFBs (Yamane et al., 2003a; Vaughan et al., 2006). Open reading frames (ORFs) between 900 bp and 1.8 kb, each containing this motif, were identified within the upstream/downstream 2 Mb regions flanking genomic RNase loci. Each ORF was used to query the GenBank NT database with an expected cutoff value of 1 × 10⁻²⁰. Resulting hits containing the terms “f-box” or “fbox” in their descriptions were treated as potential SFB genes and were combined with known S-locus associated F-box sequences from Solanaceae, Plantaginaceae, Rubiaceae, and Rosaceae. These sequences were translation-aligned using MAFFT version 7.309 and a maximum likelihood tree was constructed using RAxML version 8.2.9, with a general time-reversible model of nucleotide substitution and CAT model of rate heterogeneity. Sequences belonging to a clade that included the known S-locus F-box sequences were extracted and realigned using MAFFT. A gene tree was constructed with RAxML using general time-reversible model of nucleotide substitution and CAT model of rate heterogeneity. Rapid bootstrap analysis was conducted using ‘-f a’ option. Two thousand bootstrap replicates were obtained.

Online service for the phylogenetic placement of T2/S-RNases

Alignment and phylogenetic reconstruction with highly divergent sequences, like these diverse members of the T2/S-RNases, is a tedious process. In order to facilitate phylogenetic placement and classification of new T2/S-type RNase sequences, we provide an online service available at http://t2.karol.is. The service takes one or more nucleotide or amino acid sequences and adds them to the sequence alignment used in this study using MAFFT option ‘-add’. This new alignment is then fed into RAxML, which adds user provided sequences to the maximum likelihood tree obtained in this study using RAxML evolutionary placement algorithm, option ‘-f v’. The results, provided for download to the user, include the alignment and the tree files, as well as calculated isoelectric point values (for amino acid sequences), and classification of user sequences as putative members of class I, II, or III.

Phylogenetic relationships among T2/S-RNases

We recover three distinct clades of T2/S-RNases (Fig. 1 and Fig. S1), which mirror the previously described three ‘classes’ (e.g., Igić & Kohn, 2001; Steinbachs & Holsinger, 2002). These clades are present in the consensus and maximum credibility trees derived from MrBayes posterior tree sets, as well as the best-scoring RAxML maximum likelihood trees, inferred using both nucleotide and amino acid alignments and substitution models. (Only the best-scoring maximum likelihood tree obtained using nucleotide alignment and substitution model is used in the figures. Trees obtained using other reconstruction methods are available as supplementary material.) The posterior support values for classes I, II, & III were 1.00, 1.00, and 0.99, respectively, in the analyses of nucleotide alignment and model of sequence evolution, and 0.96, 1.00, and 0.78 when amino acid alignment and substitution model was used. The three classes are defined somewhat arbitrarily, but are remarkably well-supported by other lines of evidence examined, and we enumerate these in turn, when they relate to our principal results. No angiosperm genome we examined contains fewer than four members of this superfamily.

The dataset contained 349 sequences inferred to belong within class I. The members of this class of sequences were represented in all major land plant lineages, including the so-called ‘early-diverging’ groups—marchantiophytes, bryophytes, lycophytes, and ferns. Several well-characterized T2/S-RNases belong to class I, such as Arabidopsis thaliana RNS1 and RNS3 (Bariola et al., 1994; Bariola, MacIntosh & Green, 1999; Hillwig et al., 2008; Hillwig et al., 2011; LeBrasseur et al., 2002; Nishimura et al., 2014), Nicotiana glutinosa RNase NW and RNase NT (Kariu et al., 1998; Hino, Kawano & Kimura, 2002; Kawano et al., 2006; Kurata et al., 2002), and Solanum lycopersicum RNase LE and RNase LX (Groß, Wasternack & Köck, 2004; Jost et al., 1991; Köck et al., 2004; Köck, Stenzel & Zimmer, 2006; Lers et al., 1998; Lers et al., 2006; Löffler et al., 1992; Nürnberger et al., 1990; Tanaka et al., 2000). Many genes included in this group are secreted active RNases, expressed during senescence and phosphate starvation (see Table S3 for a list of studies of T2/S-RNase expression patterns and functions).

Class II RNases are generally found as single-copy genes within the seed plants, with the exception of recent polyploids and few apparent instances of segregating paralogous copies. The dataset contained 125 sequences placed in this class, including the genes coding for Arabidopsis thaliana RNS2 (Taylor et al., 1993; Bariola, MacIntosh & Green, 1999; Hillwig et al., 2011), Nicotiana glutinosa NGR2 (Kurata et al., 2002), and Solanum lycopersicum RNase LER (Köthke & Köck, 2011). These genes appear often constitutively expressed, and their expression levels are not necessarily increased by wounding. We did not find any class II genes in the genomes of Selaginella moellendorffii, a lycophyte, or Physcomitrella patens, a moss (bryophyte).

On the other hand, the genome of Marchantia polymorpha, a liverwort, surprisingly does contain a single class II sequence. This is unexpected because liverworts are more distantly related to flowering plants, than are either lycophytes or mosses. The discordance could be due to a variety of errors (e.g., flawed genome assemblies), independent losses in lycophytes and mosses, or unusual evolutionary processes (e.g., horizontal transfer). The possibilities could be disentangled with a broader phylogenetic coverage, but no fully sequenced fern genome has been published to date, and all T2/S-RNase sequences from ferns deposited in GenBank cluster with our class I. Finally, we found no class II sequences in the draft genomes of Salvinia cucullata and Azolla filiculoides (F-W Li, pers. comm., 2017).

Class III RNases are comprised by S-RNases and an astonishing diversity of non-S-RNase sequences. The occurrence of this group is restricted to core eudicots (The Angiosperm Phylogeny Group, 2016), although they are not present in all families whose representatives have been sequenced to date. Most notably, class III members appear absent from the well-characterized genomes in Arabidopsis and Brassica, although a distant relative, Carica papaya (Caricaceae, Brassicales), does contain a putative class III gene. Class III sequences were similarly absent from the published genome assembly of Lactuca sativa (Asteraceae).

Our dataset contained 237 class III sequences. More than a half of these (122) were specifically included because of their reportedly known S-RNase identity and function. Without exception, all functional S-RNases belong to class III, but they form an apparently polyphyletic group, although this assignment is complicated (see ‘Discussion’; Igić & Kohn, 2001). Roles of the non-S-RNase genes in this class are poorly understood, if at all. Based solely on the high diversity of primary sequence features, and expression patterns, they appear potentially highly functionally disparate.

Estimates of ancestor-descendant relationships among the three classes are challenging without additional information, because of the high sequence divergence and uncertainty over the prior expectation for rooting lineage(s). In the analyses using nucleotide alignments and models of sequence evolution, plausible outgroup sequences from algae disrupt the monophyly of land plant RNases, very possibly as an artifact of deep divergence (over 450 My), at the limits of inference for a relatively short gene (ca. 600 bp). Two sequences from Volvox carteri, and one from Chlamydomonas reinhardtii form a clade that does not include the sequence from Bryopsis maxima. Rooting the trees using the Volvox/Chlamydomonas clade instead results in monophyly of classes II and III, while rooting the tree using the Bryopsis sequence results in monophyly of classes I and III. Representative genes from algae are monophyletic in the trees produced using amino acid model of sequence evolution (AA trees). Placing the root between the algal clade and land plant RNases, class II sequences are sister to classes I and III. None of the three possible arrangements among the classes receives significant support. The most apparent incongruence found in all RNase trees—independent of reconstruction method used—is the placement of class I lycophyte, fern, and gymnosperm sequences within the angiosperm RNase clades (Fig. S1d). However, the posterior support values for this arrangement vary wildly (0.01 and 0.84 for best-scoring maximum likelihood trees inferred using amino acid and nucleotide data, respectively). Additionally, in analyses with nucleotide data, gymnosperm class II genes are placed within monocot sequences with posterior support of 0.69 (Fig. S1a).

We identified a set of 91 sequences either corresponding to proteins established to be catalytically inactive in functional studies (MacIntosh et al., 2010; Ohkama-Ohtsu et al., 2004; Wei et al., 2006; Gausing, 2000; Van Damme et al., 2000; Kim et al., 2004), or lacking a conserved histidine residue that is essential for ribonuclease activity (Jost et al., 1991; Taylor et al., 1993; Irie, 1999).

The placement of these sequences is shown in Fig. 2. They are distributed across the tree, clustering in several small clades, often associated with changes in pI value (see below). In our dataset, class I contains 67 such inferred non-functional sequences, scattered across nine clades, class II contains three sequences in three clades, and class III contains 21 sequences in 12 clades (Fig. 2). The pattern is consistent with independent losses of ribonuclease function, perhaps following gene duplications or losses of SI.

Intron positions

Apart from several scattered and easily identifiable recent gains and losses, the patterns of intron presence and absence are well-conserved and remarkably concordant with our T2/S-RNase gene trees (Fig. 2). We identified 11 observed intron positions across all land plants (embryophytes) examined. The position, phase (reading frame), and numerical abundance of introns are summarized in Fig. 3. Introns at a given position, which we supposed to share ancestry based on extremely similar position in sequence alignment (within seven nucleotides in final alignment), did not exhibit any phase variation across land plants. Such conservation of exon-intron boundary positions within codon triplets reinforces our classification of intron occurrences. We also noted the phylogenetic distribution and number of genes matching each intron pattern in species with complete genomes (Fig. 3).

Intron patterns are largely concordant with the phylogenetic classes. Class I genes display four intron presence–absence patterns. Most sequences contain three introns at positions 2, 5, and 9 (pattern I-D). Absence of intron at position 9 defines the second most frequent pattern, I-C. Two additional single intron patterns, I-A and I-B, containing introns at positions 2 and 5, respectively, occur only in Poaceae (grass family). Class II sequences are remarkable, in that all 46 sequences examined to date contain eight introns, and exhibit no apparent variation in intron pattern. Most class III sequences contain a single intron at position 5 (pattern III-C). Almost all known S-RNases exhibit this intron pattern although Prunus S-RNases contain an additional intron at position 1 (pattern III-D). Pattern III-A (intron positions 5 and 9) is found in the sequences of several distantly related Rosid species (Ricinus communis, Carica papaya, Cajanus cajan, Glycine max, Theobroma cacao, Gossypium arboreum, Gossypium raimondii). Although the dataset contained 21 sequences with this intron pattern, most of these were paralogs. Eleven copies were found in the genome of Theobroma cacao, and other species contained between one and three copies. Pattern III-B^* contains a single intron at position 9. This pattern has no EST support and was predicted based on genomic sequences from Fragaria nubicola and Fragaria vesca.

Intron positions and their phases are highly conserved. Position 5 is found in all three plant T2/S-RNase classes, as expected (Igić & Kohn, 2001). Position 9 introns are shared by some class I and III members. All other intron positions are class-specific. Position 2 occurs only in class I sequences and seven intron positions (3, 4, 6, 7, 8, 10, and 11) are unique to class II. Overall, merely eight distinct intron position patterns exist, seven of which are class-specific. One apparent exception—likely due to convergence—is the shared pattern of a single intron at position 5 in nine grass sequences belonging to class I and many sequences, especially S-RNases, in class III (patterns I–B and III–C).

None of the intron patterns found in land plants appear located in identical positions as those in algal, fungal, or animal genes. Sequences from the Volvox carteri and Bryopsis maxima contain four introns, three of which may be in the same ancestrally shared positions as the ones found in land plant sequences. Although their phases differ, and exact location appears to be slightly shifted, each species contains an intron possibly ancestrally shared at position 2 in plant T2/S-RNases. Introns in the second reading phase (+2)—a potential homologue to one found at position 3 in land plants—is also present, although our sequence alignment is ambiguous in this region. Another intron appears to be homologous to position 8. It is in the same phase (+0 phase; not interrupting the reading frame) as the ones found in land plant T2/S-RNases.

A limited sample of animal sequences examined contain six to eight introns, two of which are potentially in identical sites to those we find in the land plants. One of these, near position 5, as in plants, is in phase +0. The other, near position 8, is in phase +2, while plant introns at this position are in phase +0. Two out of four introns found in fungal sequences may be homologs of the ones found in plant sequences. An intron from Aspergillus oryzae is near position 8 and is in the same phase (+0) as land plant introns. The other intron from this species, near position 10, is in phase +1, compared to phase +0 of plants. It is unclear how quickly intron sliding and phase evolve, generally, and we have little statistical evidence to establish clear links between fungal, animal, and plant RNases in our dataset.

Isoelectric point (pI) values

We estimated isoelectric point (pI) values for T2/S-RNases as vague heuristic indicators of possible subcellular localization and function, with a possibly informative pattern of values across the gene family tree (Drawid & Gerstein, 2000; Kirkwood et al., 2015). The predicted pI values in our dataset range from 3.91 to 9.91. Although the predicted pI value ranges of all three RNase classes largely overlap, the values show distinct class-specific trends (Fig. 4). Class I peptides generally have acidic pI values, with a median of 5.04 (min = 3.91; max = 9.75; n = 349). Class II peptides have similarly acidic pI values, with a median of 5.90 (min = 4.52; max = 9.03; n = 125). By way of contrast, class III peptides are significantly more basic. The non-S-RNases (or unknown function peptides) have a median pI of 8.56 (min = 4.61; max = 9.91; n = 118), while S-RNases have a median pI of 9.18 (min = 8.10; max = 9.73; n = 119).

In order to examine the distribution of pI values across T2/S-RNases, we mapped the predicted sequence pI values on the corresponding gene tree (Fig. 2). Class I contains nine independent pI shifts from acidic to basic values, seven of these were associated with the loss of the active histidine residue. Class II contains eleven independent pI shifts from acidic to basic values, none of these were associated with the loss of the active histidine residue. Class III contains 22 independent pI shifts from basic to acidic values, six of these were associated with the loss of the active histidine residue. The apparent conservation and concordance of pI is fairly remarkable, given its lack of clear relationship with protein function (Drawid & Gerstein, 2000; Brett, Donowitz & Rao, 2006).

F-box domain-containing genes near T2/S-RNase loci

We searched the available genome assemblies for the newly identified class III T2/S-RNase family members, without a known function, in an attempt to find whether they are co-located with F-box motif-containing genes (within 2 Mb). We reasoned that such associations may have comprised—or still comprise—a functional S-locus. But the resulting picture is complex. Many class I and class II RNases also contain F-box motif-containing genes within 2 Mb, which is perhaps unsurprising given their abundance in plant genomes (Wang et al., 2004). Nevertheless, the structure of putative S-loci ought to resemble the canonical pattern: a class III T2/S-RNase, accompanied by a more than a few F-box-containing genes.

We identified genomic regions that contain a class III RNase and multiple F-box loci in five eurosid species, not including previously characterized S-loci (Fig. 5 and Table S2). Citrus clementina (Rutaceae) genome contains a class III RNase and six F-box loci on scaffold 5. This genomic region, however, was largely unresolved (46% composed of ambiguous characters), which may have prevented the discovery of more F-box loci. The genomes of Theobroma cacao, Gossypium raimondii, and Gossypium arboreum (Malvaceae) contain a class III RNase and 15, 14, and 14 F-box loci on scaffold 10r, chromosome 11, and chromosome 10, respectively. Phaseolus vulgaris (Fabaceae) contains a class III RNase and eight F-box loci on chromosome 4, as in C. clementina, however, 45% of the genomic region investigated was unresolved. Two or more F-box loci in each of these haplotypes contain in-frame stop codons, which is also the case in the non-functional haplotype of domesticated tomato, Solanum lycopersicum (Table S2). Unlike canonical S-haplotypes, however, all of these contain more than one class III RNase; Citrus clementina, Theobroma cacao, and Gossypium raimondii contain two, while Gossypium arboreum and Phaseolus vulgaris contain three.

We also examined evidence for an alternative cause of co-location of RNases and F-boxes, an unspecified functional constraint that causes an association of RNases and F-boxes. Interestingly, as is the case in Arabidopsis thaliana (Wang et al., 2004), genomic regions flanking class I and II T2/S-RNases from several species also contain multiple F-box loci (bottom two haplotypes in Fig. 5). Genomic segments (2 Mb) flanking the well-characterized Solanum lycopersicum RNases LX and LE (class I T2/S-RNases), which occur in tandem on chromosome 5, contain at least seven F-box loci. Similarly, Gossypium arboreum class II locus on chromosome 4 has at least six F-box loci. Specifically, all F-box sequences flanking class III RNases in Fig. 5 cluster within F-box groups that contain known S-locus F-box genes (Fig. S3), while the F-box genes associated with Solanum lycopersicum class I and Gossypium arboreum class II RNases cluster outside this clade.

T2/S-RNases in land plants

T2/S-RNases comprise a ubiquitous family of endoribonucleases, often found in low copy number in the genomes across all domains of life, with the sole exception of Archaea (Irie, 1999; MacIntosh, 2011). Their wide distribution alone suggests conservation of important function(s). A common and possibly ancestral role appears to be ribosomal RNA decay and recycling (Hillwig et al., 2009; Ambrosio et al., 2014). In a variety of species, T2/S-RNases are induced under oxidative stress and function in tRNA cleavage, which appears to be a conserved response in eukaryotes (Thompson et al., 2008; MacIntosh et al., 2001; Ambrosio et al., 2014). In plants, their main function appears to be phosphate harvesting from degraded RNA (MacIntosh, 2011). Indeed, phosphorus is a limiting nutrient for plants, so intracellularly abundant rRNA and senescing tissues comprise important recyclable resources (Bariola et al., 1994).

We find that seed plant genomes feature an expanded repertoire of T2/S-RNases, compared with other organisms, so that each diploid genome contains four or more members of this gene family. The causes of this expansion are difficult to infer, but there is an intriguing possibility that it accompanied the invasion of land and the subsequent development of vasculature and increase in stature, especially as plants moved away from steady sources of dissolved available inorganic phosphates. It is inviting to further speculate that unicellular plants are unlikely to have many members of the T2/S-RNases. Multicellular plants, especially those with a variety of metamers and displaying diverse organ identity, may rely on separate subfunctionalized paralogs for efficient phosphate recycling and recruitment. They would also be far more prone to diversification into neofunctionalized paralogs with unrelated roles, such as sexual self-recognition and defense. The exact distribution and patterns of diversification of T2/S-RNases in land plants remain unclear, but increasing genome sequencing coverage across plants is likely to help focus the study on the processes than may have been responsible for shaping the many roles taken on by this gene family.

Implications for the evolution of self-incompatibility

The distribution of both the S-RNase-based SI and class III RNases has so far been restricted within the core eudicots. Combined with the shared intron presence–absence patterns found in S-RNases and similarities in the male components of SI, this provides considerable evidence for the single origin of S-RNase-based SI. Class III sequences are not found in all core eudicots, and do not appear to be strictly essential in any organism in which they were studied to date.

It is perhaps significant that they are absent from several sequenced species in the genera Arabidopsis and Brassica, which both express a sporophytic SI system (or are otherwise self-fertile). This observation suggests that class III RNases originated before the divergence of core eudicots, and may be maintained due to their role in RNase-based SI. It is possible then, that novel functions of proteins found within class III originated from S-RNase paralogs, with current utility unrelated to sexual systems. Another possibility has gained traction.

Starting with Sassa et al. (1996), many have proposed that S-RNase-based SI in Rosids and Asterids evolved independently, based on the non-monophyly (polyphyly) of S-RNases, although their conclusions contrasted with two contemporary analyses (Xue et al., 1996; Richman, Broothaerts & Kohn, 1997). This argument is consistent with repeated recruitment of class III T2/S-RNases for a role in SI. Later work involved a steadily increasing number of sequenced genes, and clarified that strict monophyly of S-RNases is not a necessary condition for shared ancestry (e.g., Steinbachs & Holsinger, 2002). For example, gene duplication, followed by functional changes, could easily account for the widespread occurrence of paralogs, and yield non-monophyly of functional S-RNases, as could a great variety of sources of error in phylogenetic inference. Moreover, subsequent studies also found that the male component of SI response in this system, SFBs, are expressed in the pollen of species from each of the well-studied families (Zhou et al., 2003; Ushijima et al., 2003; Sijacic et al., 2004), a development that appeared to affirm the single-origin hypothesis beyond doubt. The debate has now acquired a new twist.

It has come to light that an increasing number of differences distinguish the mechanism of action of S-RNase-based SI in Prunus (Rosaceae) from that acting in the relatively closely related subtribe Malinae (Rosaceae) and in the euasterid families Solanaceae and Plantaginaceae. The inferred dissimilarities include the mode of recognition (self- vs. non-self; e.g., Ushijima et al., 2004; Fujii, Kubo & Takayama, 2016), phenomenology and causes of breakdown of SI (Golz et al., 2001; Ushijima et al., 2004; Hauck et al., 2006; Xue et al., 2009), magnitude of selection and sites experiencing it (Ashkani & Rees, 2016), S-locus structure (Akagi et al., 2016), as well as the patterns of divergence and relationships (Kohn, 2008; Akagi et al., 2016) among both S-RNases and F-box-containing genes (SLF/SFB). A range of novel evolutionary scenarios have been proposed to explain them, ranging from surprising wholesale convergence to homology with divergence (Morimoto, Akagi & Tao, 2015; Aguiar et al., 2015a; Akagi et al., 2016).

On the other hand, it seems particularly likely that precise models describing the evolution of RNase-based SI remain elusive in part due to its apparent variation and complexity. The problem is compounded by a notable lack of detailed functional studies outside of a few species in Solanaceae, Plantaginaceae, and Rosaceae, which ensures that we have very little ability to shape reasonable expectations regarding the capacity of such systems to undergo the kinds of divergence observed over the relevant timescales.

Even the most sophisticated evolutionary analyses can merely supply a scorecard of similarities and differences in the discovered components of RNase-based SI, in the absence of context provided by extensive comparative data. The question is somewhat academic, in the sense that there are well-documented examples of the incredible capacity for change among molecular components underlying conserved traits. It is possible to imagine a total turnover of genetic components, which would eliminate many outward diagnostic signs of homology. This process is clearly exemplified by the repeated independent co-option of unrelated genes as lens crystallins in different vertebrate lineages (True & Carroll, 2002). Replacement or modification of one or more genetic components of a complex system by co-option of unrelated genes does not necessarily interrupt its genealogical continuity and function. Therefore, just as we discarded the expectation of strict monophyly among the molecular components of SI, due to the vast capacity of genes to undergo duplication and subsequent changes, perhaps we should do the same with the expectation of identity for all functional details of SI response for all lineages. It seems that the main difficulty concerns the development of a common framework of approaches that can delineate ‘deep homology’ (Shubin, Tabin & Carroll, 1997) from convergence, comprised of re-recruitment of similar components.

The presently employed framework or haphazard data collection—and often flawed analyses—from a number of unrelated lineages is insufficient for detailing this instance of deep homology. The growing list of differences and interacting units that cause RNase-based SI, each with possible unknown pleiotropic effects, is increasing the complexity of our task. In Brassicaceae, where a distinct kinase-based mechanism operates, we are aware of one complicating instance within the family. Species in the genus Leavenworthia appear to contain two paralogous loci, one of which apparently encodes S-allele phenotypes, and another with an as yet undetermined function (Chantha et al., 2013). As the authors point out, such a finding illustrates the vastness of the problems before us, but it does not necessitate rejection of homologous ancestry of the trait across the family. In the present study, with perhaps the most extensive collection of data on female and male components of RNase-based SI at hand (to date), we likewise find that there is little evidence to overturn the long-standing hypothesis of a single RNase-based GSI system origination predating the common ancestor of rosid and asterid eudicots. This conclusion, like those of its detractors, is necessarily rife with potential problems.

Challenges for the inference of history of T2/S-RNases and self-incompatibility

Phylogenetic analyses of the components that interact to affect SI responses are vulnerable to a number of sources of error and bias (Felsenstein, 2004). Most studies that posit some manner of convergent re-evolution of RNase-based SI principally rely on the precise relationships of S-RNase (or SLF/SFB) gene trees. It is trivially unsurprising that a particular group of functional S-RNases with a shared history may be recovered as para- or polyphyletic, even under perfectly specified inference models. Sequence-based inferences of evolutionary relationships among very distantly related genes should be viewed with some skepticism. In the particular case of the T2/S-RNase family in angiosperms, our focus here, the available gene sequences are relatively short and highly divergent, so that significant loss of information is expected due to the accumulation of multiple changes per site, even under correctly specified models of sequence evolution. Under such conditions there may simply not be enough information to accurately recover the historical relationships. This uncertainty may not be reflected in node support values due to methodological artifacts, such as long branch attraction—a spurious, yet confident, association of distantly related sequences.

More seriously, inference bias and error can result from a vast range of unspecified evolutionary processes, such as gene duplication, large (multi-nucleotide) subsequent changes or loss. Phylogenetic tree inference is conditional on both the correct data (sequences and their alignment) and the model of sequence evolution. While we may be somewhat confident that the gene sequences are adequate, the assessment of alignment accuracy is less trivial, especially with high sequence divergence (Felsenstein, 2004; Kumar et al., 2012). Substitution models encompass only a narrow range of biologically possible processes, and they do not easily accommodate indels, recombination, variation of substitution rates over time (heterotachy) and between clades (Kumar et al., 2012). Perhaps critically for analyses involving the numerous SLF/SFB paralogs, the presently used models of tree inference do not accommodate gene conversion, which can result in spectacular model-misspecification, and subsequently erroneous inference. Such models exist (Song et al., 2011), but are not easily integrated into the common workflows. Despite these difficulties, molecular phylogenetic approaches remain indispensable, and are often the only hypothesis generating tools available. The task of reconstructing evolutionary events on timescales of ca. 50–100 My ought to be daunting, carefully framed, and generally include circumspect qualification of the resulting analyses.

Nevertheless, even in the presence of a variety of flaws, evaluation of a variety of protein features in phylogenetic context may help to narrow down the list of candidate T2/S-RNases for functional characterization. In this vein we provide an on-line service (http://t2.karol.is), which places user-provided sequences on the T2/S-RNase phylogeny used in this study. If amino acid sequences are provided, pI values will be calculated. Users of such automated workflows would do best not to ignore a variety of possible problems, as outlined above.