Higher Rates of Protein Evolution in the Self-Fertilizing Plant Arabidopsis thaliana than in the Out-Crossers Arabidopsis lyrata and Arabidopsis halleri.

The common transition from out-crossing to self-fertilization in plants decreases effective population size. This is expected to result in a reduced efficacy of natural selection and in increased rates of protein evolution in selfing plants compared with their outcrossing congeners. Prior analyses, based on a very limited number of genes, detected no differences between the rates of protein evolution in the selfing Arabidopsis thaliana compared with the out-crosser Arabidopsis lyrata. Here, we reevaluate this trend using the complete genomes of A. thaliana, A. lyrata, Arabidopsis halleri, and the outgroups Capsella rubella and Thellungiella parvula. Our analyses indicate slightly but measurably higher nonsynonymous divergences (dN), synonymous divergences (dS) and dN/dS ratios in A. thaliana compared with the other Arabidopsis species, indicating that purifying selection is indeed less efficacious in A. thaliana.

In plants, the transition from out-crossing to self-fertilization is quite common, and is generally seen as a dead-end due to accumulation of deleterious mutations (Stebbins 1957). Population genetics theory predicts that selfing organisms will have a lower effective population size (Ne) than their outcrossing congeners with the same population size (Pollak 1987). Reduced Ne is expected to result in a reduced efficacy of natural selection (Charlesworth et al. 1993; Charlesworth and Wright 2001), thus allowing the fixation of slightly deleterious mutations (Ohta 1973). As a result, selfing organisms are expected to exhibit accelerated rates of protein evolution (Kimura 1983; Charlesworth and Wright 2001; Glémin 2007) and less codon usage bias (Qiu et al. 2011). These predictions are supported by some empirical evidence: natural selection is reduced in selfing species of the family Triticeae (Haudry et al. 2008; Escobar et al. 2010) and the genera Capsella (Qiu et al. 2011; Slotte et al. 2013), Eichhornia (Ness et al. 2012), Collinsia (Hazzouri et al. 2013), and Mimulus (Brandvain et al. 2014) (for review, see Hough et al. [2013] and Shimizu and Tsuchimatsu [2015]). In addition, an analysis of polymorphism data for a number of plant species revealed a weak increased in the nonsynonymous to synonymous polymorphism ratio (π/π) of selfers (Glémin et al. 2006).

The plant Arabidopsis thaliana is thought to have shifted to self-fertilization 150,000–1,000,000 years ago (Charlesworth and Vekemans 2005; Bechsgaard et al. 2006; Tang et al. 2007; Tsuchimatsu et al. 2010; Shimizu and Tsuchimatsu 2015; Durvasula et al. 2017; Tsuchimatsu et al. 2017). In agreement with the predicted reduction in the efficacy of natural selection, this species exhibits less codon bias than the out-crosser Arabidopsis lyrata (Qiu et al. 2011). However, analysis of 16 genes did not detect significantly higher rates of protein evolution in A. thaliana compared with A. lyrata (Wright et al. 2002). In addition, comparison of 13 pairs of orthologous genes in these two species revealed no differences in the ratios of nonsynonymous to synonymous polymorphisms or in the ratios of nonsynonymous to synonymous fixations (Foxe et al. 2008). A comparison of 675 A. thaliana and 73 A. lyrata nonorthologous genes found higher ratios of nonsynonymous to synonymous polymorphisms and higher ratios of nonsynonymous to synonymous fixations in A. thaliana (Foxe et al. 2008); however, these results may have been affected by biases in the data set—for example, seven of the A. lyrata genes were chosen due to their high levels of expression, and highly expressed genes tend to evolve under strong purifying selection (Pál et al. 2001; Drummond et al. 2005).

These analyses, in any case, were limited by the very small amount of genomic information available at the time. Here, we revisit the prediction that A. thaliana should exhibit faster rates of protein evolution than A. lyrata or than Arabidopsis halleri taking advantage of the now completely sequenced genomes of A. thaliana (The Arabidopsis Genome Initiative 2000), A. lyrata (Hu et al. 2011), A. halleri (Briskine et al. 2017) and the outgroup Capsella rubella (Slotte et al. 2013). A. thaliana diverged 6–13 Ma ago from the A. lyrata/A. halleri clade (Beilstein et al. 2010; Hohmann et al. 2015) and 8–14 Ma ago from C. rubella (Koch and Kiefer 2005; Hohmann et al. 2015) (fig. 1).

. 1.

—Phylogenetic relationships among the species used in the current study. The tree topology and divergence times were obtained from Hohmann et al. (2015).

For each C. rubella gene, we identified the most likely orthologs in A. thaliana and A. lyrata. For each of the 18,107 identified trios, protein sequences were aligned, and the resulting alignments were used to guide the alignment of the corresponding coding sequences (CDSs). To reduce the impact of annotation errors, we removed all alignments for which >5% of positions included gaps. For each of the resulting 12,994 alignments, PAML (free-ratios model; Yang 2007) was used to estimate the nonsynonymous divergence (dN), synonymous divergence (dS) and the nonsynonymous to synonymous divergence ratio (ω = dN/dS) in each of the branches of the phylogeny (fig. 1). The ratio dN/dS is expected to be lower than 1 when nonsynonymous mutations are under purifying selection (with values closer to 0 indicating stronger selection), equal to 1 when protein sequences evolve neutrally, and higher than 1 for genes under positive selection (for review, see Alvarez-Ponce [2014]).

In the A. thaliana branch, the median of the values estimated by the free-ratios model were dN = 0.0108, dS = 0.0757, and dN/dS = 0.1427, and the mean values were dN = 0.0133, dS = 0.0805, and dN/dS = 0.1865 In the A. lyrata branch, the median values were dN = 0.0085, dS = 0.0612, and dN/dS = 0.1389, and the mean values were dN = 0.0107, dS = 0.0667, and dN/dS = 0.1880 (supplementary table S1, Supplementary Material online and fig. 2). A Mann–Whitney U test showed significant differences in the dN (P = 1.964 × 10−119), d (P < 10−300), and dN/dS (P = 0.0127) of both species. In 8,572 of the cases, dN was higher in A. thaliana than in A. lyrata, and in 4,396 of the cases dN was higher in A. lyrata, indicating that rates of protein sequence evolution are often higher in A. thaliana (binomial test, P = 1.20 × 10−324). In 8,938 of the cases, dS was higher in A. thaliana, and in 4,055 of the cases dS was higher in A. lyrata, indicating faster rates of evolution of synonymous sites in A. thaliana (binomial test, P = 4.94 × 10−324); these results are consistent with prior studies reporting higher mutation rates in A. thaliana (Yang et al. 2013). In 6,625 of the cases, dN/dS was higher in A. thaliana, and in 6,161 of the cases dN/dS was higher in A. lyrata, indicating that purifying selection on protein sequences is often less effective in A. thaliana (binomial test, P = 4.22 × 10−5). Differences were stronger when analyses were restricted to genes that are highly expressed in A. thaliana (supplementary table S2, Supplementary Material online).

. 2.

—Distribution of dN, dS, and dN/dS values in the A. thaliana and A. lyrata branches. Values above the 90th percentile are not represented.

For each alignment, Tajima’s relative rate test (Tajima 1993) was used to contrast whether the number of substitutions accumulated in A. thaliana and A. lyrata was significantly different. Statistically significant differences were detected in 1,363 and 1,333 genes for synonymous and nonsynonymous sites, respectively. Of the 1,363 genes with significant differences in synonymous rates of evolution, there were more unique synonymous changes in A. thaliana in 1,222 genes compared with 141 genes where A. lyrata had more unique synonymous changes. Of the 1,333 genes with an asymmetry in the number of nonsynonymous sites, A. thaliana and A. lyrata had more unique changes in 1,077 and 256 cases, respectively (table 1).

Table 1

Tajima's Relative Rate Tests

	All Substitutions	Synonymous Substitutions	Nonsynonymous Substitutions
Unique substitutions in A. thaliana	413,215	249,860	163,355
Unique substitutions in A. lyrata	343,062	205,266	137,578
Genes where A. thaliana had more substitutions	9203	8701	7575
Genes where A. lyrata had more substitutions	3207	3424	4102
Genes where P < 0.05	2008	1363	1333
Genes where P < 0.05 and A. thaliana had more substitutions	1824	1222	1077
Genes where P < 0.05 and A. lyrata had more substitutions	184	141	256
χ2 value for concatenome	6507.5	4369.4	2208.0
P value for concatenome	≪0.001***	≪0.001***	≪0.001***

P < 0.001.

For each of the 12,994 alignments, we compared the likelihood of the free-ratios model (in which each of the three branches exhibits an independent dN/dS ratio) versus that of a two-ratios model (one dN/dS ratio for A. thaliana and A. lyrata, and another for C. rubella). The free-ratios model fitted the data significantly better in 907 of the alignments (likelihood ratio test, P < 0.05), indicating that the dN/dS ratio is significantly different in A. thaliana and A. lyrata. In 477 of the 907 cases where the free-ratio model fit better than the two-ratios model, dN/dS was higher for A. thaliana, and in 430 of the cases dN/dS was higher for A. lyrata; these numbers were not significantly different from the 50%:50% (453.5:453.5) expected by chance (binomial test, P = 0.166).

Given that A. thaliana and A. lyrata are very closely related, some gene alignments may not contain sufficient information (in terms of number of substitutions) to accurately infer the strength of purifying selection acting on each branch. In order to increase the power of our analyses, we combined all 12,994 alignments into a single concatenome containing 17.8 million base pairs and repeated our analyses on it. The A. thaliana lineage exhibited higher dN, dS and dN/dS values (0.0127, 0.0759, and 0.1671, respectively) (supplementary table S3, Supplementary Material online) than the A. lyrata branch (0.0102, 0.0622, and 0.1644). These values are comparable to the mean values resulting from analysis of individual alignments. The free-ratios model fitted the data significantly better than the two-ratios model (2Δℓ = –10.213, P = 0.0014), indicating that dN/dS is significantly higher in A. thaliana, even though the differences are small. Tajima’s relative rate test (Tajima 1993) revealed an excess of synonymous and nonsynonymous changes in A. thaliana compared with A. lyrata (χ= 4369.4 and 2207.0, P ≪ 0.001 and P ≪ 0.001, respectively). The A. thaliana concatenome contained 249,860 synonymous 163,355 nonsynonymous substitutions that were not present in A. lyrata, and the A. lyrata concatenome contained 205,266 unique synonymous substitutions and 137,578 unique nonsynonymous substitutions.

It is expected that the evolution of selfing in A. thaliana may have resulted in pseudogenization of, or at least relaxation of purifying selection in, genes involved in out-crossing. If these represent a sufficiently large number of genes, this effect alone, rather than a reduction of Ne, might conceivably explain the higher average rates of protein evolution observed in A. thaliana. To discard this possibility, we repeated our analyses separately for genes of different functional categories. For all 23 KOG categories represented in the data set, the number of genes with higher dN and dS values in A. thaliana was significantly higher than the number of genes with higher dN and dS values in A. lyrata. For 19 of the categories, the number of genes for which dN/dS was higher in A. thaliana was higher than the number of genes for which dN/dS was higher in A. lyrata. For only three categories there were more genes with a higher dN/dS in A. lyrata (binomial test, P = 0.0009; table 2). These results indicate that the higher dN, dS and dN/dS values observed in A. thaliana represent a generalized trend, not specific to certain functional categories.

Table 2

Analyses of Evolutionary Rates in Different KOG Categories

Categorya	Genes with Higher dN in A. thaliana	Genes with Higher dN in A. lyrata	Genes with Higher dS in A. thaliana	Genes with Higher dS in A. lyrata	Genes with Higher dN/dS in in A. thaliana	Genes with Higher dN/dS in in A. lyrata	dNP Valueb	dSP Valueb	dN/dSP Valueb
A	245	99	248	97	175	163	2.03 × 10−15***	2.2 × 10−16***	0.550
B	77	32	78	32	62	45	1.94 × 10−5***	1.36 × 10−5***	0.122
C	256	138	257	137	209	174	2.89 × 10−9***	1.53 × 10−9***	0.082
D	121	50	124	48	85	84	5.64 × 10−8***	6.13 × 10−9***	1.000
E	189	113	216	86	152	150	1.44 × 10−5***	4.63 × 10−14***	0.954
F	60	32	60	32	53	36	0.005**	0.0046**	0.089
G	500	308	566	245	389	422	1.47 × 10−11***	<2.2 × 10−16***	0.261
H	116	74	130	60	94	94	0.003**	4.16 × 10−7***	1.000
I	225	125	243	107	184	166	9.99 × 10−8***	2.69 × 10−13***	0.364
J	214	114	201	127	178	135	3.63 × 10−8***	5.20 × 10−5***	0.017*
K	697	362	720	344	541	511	<2.2 × 10−16***	<2.2 × 10−16***	0.371
L	166	68	161	73	119	113	1.23 × 10−10***	8.70 × 10−9***	0.743
M	104	67	121	50	78	92	0.006**	5.64 × 10−8***	0.319
O	710	317	743	287	507	494	<2.2 × 10 − 16***	<2.2 × 10−16***	0.704
P	340	167	367	140	262	243	1.22 × 10−14***	<2.2 × 10−16***	0.423
Q	227	115	241	101	189	153	1.40 × 10−9***	2.45 × 10−14***	0.058
S	2344	1213	2405	1156	1784	1751	<2.2 × 10−16***	<2.2 × 10−16***	0.590
T	673	335	736	272	517	479	<2.2 × 10−16***	<2.2 × 10−16***	0.241
U	353	182	389	150	269	241	1.24 × 10−13***	<2.2 × 10−16***	0.232
V	52	28	54	26	44	35	0.001**	0.002**	0.368
W	48	26	51	23	34	40	0.014*	0.001**	0.561
Y	22	3	19	6	17	8	1.57 × 10−4***	0.015*	0.108
Z	159	58	155	63	119	92	4.64 × 10−12***	3.94 × 10−10***	0.073

Category functions: A, RNA processing and modifications; B, chromatin structure and dynamics; C, energy production and conversion; D, cell cycle control, cell division, chromosome partitioning; E, amino acid transport and metabolism; F, nucleotide transport and metabolism; G, carbohydrate transport and metabolism; H, coenzyme transport and metabolism; I, lipid transport and metabolism; J, translation, ribosomal structure and biogenesis; K, transcription; L, replication, recombination, and repair; M, cell wall, cell membrane and envelope biogenesis; O, posttranslational modification; P, inorganic ion transport and metabolism; Q, secondary metabolite biosynthesis; S, function unknown; T, signal transduction; U, intracellular trafficking, secretion, and vesicular transport; V, defense mechanisms; W, extracellular structures; Y, nuclear structure; Z, cytoskeleton (Tatusov et al. 2003).

P-values determined using a binomial test comparing the total number of genes where dN/dS was higher in A. thaliana and in A. lyrata.

P < 0.05; **P < 0.01; ***P < 0.001.

Throughout the current work we have reported the comparison of the A. thaliana reference genome from the TAIR 10 release, a composite genome from 11 Columbia ecotype (Col-0) individuals, with that of A. lyrata, using C. rubella as outgroup. Nonetheless, equivalent results were obtained using another 18 A. thaliana accessions instead of the reference one (supplementary tables S4 and S5, Supplementary Material online), using A. halleri (Briskine et al. 2017) instead of A. lyrata (supplementary tables S6–S8, Supplementary Material online) or using the outcrossing and more distantly related Thellungiella parvula (Dassanayake et al. 2011) as outgroup instead of the selfing and closely related C. rubella (supplementary tables S9 and S10, Supplementary Material online).

In summary, all our genome-wide analyses converge at showing that, as expected from the reduced Ne due to selfing, proteins evolved faster in A. thaliana than in A. lyrata or A. halleri. Such protein sequence evolution acceleration is likely due to the combination of faster mutation rates in A. thaliana (supported by high dS values and by prior results; Yang et al. 2013) and by a weaker efficacy of natural selection on nonsynonymous mutations (supported by high dN/dS ratios). Prior analyses based on a handful of orthologous genes failed to detect differences in dN and dN/dS between A. thaliana and A. lyrata, most likely because of limited statistical power (Wright et al. 2002; Foxe et al. 2008). Indeed, the differences that we detected are subtle, consistent with the fact that A. thaliana has been selfing for a relatively short amount of time (150,000–1,000,000 years; Charlesworth and Vekemans 2005; Bechsgaard et al. 2006; Tang et al. 2007, Durvasula et al 2017) compared with the time of divergence between A. thaliana and the A. lyrata/A. halleri clade (7–13 Myr; Beilstein et al. 2010; Hohmann et al. 2015). Our analyses have compared the patterns of evolution of the A. thaliana lineage (the branch connecting A. thaliana and the most recent common ancestor of A. thaliana and A. lyrata) and the A. lyrata and A. halleri lineages (the branches connecting A. lyrata or A. halleri and the most recent common ancestor of A. thaliana and A. lyrata), and plants in the A. thaliana lineage have been selfing for only 1–17% of the length of the branch.

In addition to the recent transition to selfing of A. thaliana, other scenarios may account for the small magnitude of the differences observed between the rates of protein evolution of A. thaliana and A. lyrata. First, most proteins are under strong purifying selection in both species, in agreement with prior observations (Wright et al. 2002; Foxe et al. 2008. Yang and Gaut 2011), thus hindering the detection of strong differences. Second, selfing increases homozygosity, thus exposing recessive alleles to selection, which can reduce rates of protein evolution (see Glémin 2007). Last, population genetics analyses indicate that the Ne of A. lyrata may have also been reduced within the last 100,000 years (Mattila et al. 2017); this might have increased the rates of protein evolution in this species, thus attenuating the differences between A. thaliana and A. lyrata.

Finally, it should be noted that the fast rates of protein evolution observed in A. thaliana might be due to peculiarities of the biology of this species other than selfing. In particular, A. thaliana switched to an annual life history, whereas A. lyrata is perennial. Annual plants tend to evolve faster than perennial plants (Smith and Donoghue 2008; Gaut et al. 2011; Lanfear et al. 2013), which might account for the higher rates of synonymous evolution observed in A. thaliana. However, annual plants exhibit lower nonsynonymous to synonymous polymorphism ratios (Chen et al. 2017), and thus the annual life history of A. thaliana may not explain the observed dN/dS ratios observed in this species.

Materials and Methods

For each C. rubella gene, the longest encoded protein was chosen for analysis and orthologs in A. thaliana and A. lyrata were identified using a best reciprocal hit approach (BLASTP, E-value < 10−10). Only genes for which orthologs could be identified in both Arabidopsis species were retained. Trios of orthologous protein sequences were aligned using PRANK v.140603 (Löytynoja and Goldman 2005), and the resulting alignments were used to guide the alignments of the CDSs using an in-house script. Alignments which contained <5% gaps were retained for analyses. For each alignment, the codeml program of PAML v. 4.9 (Yang 2007) was used to estimate dN, dS and dN/dS in each of the three branches (free-ratios model) and in the A. thaliana/A. lyrata branch and the C. rubella branch separately (two-ratios model). Values of dN/dS above ten were removed from mean calculations, in order to prevent the bias introduced by these outliers, which represent artifacts due to the presence of very few mutations in the relevant lineages. The fit of both nested models was compared using a likelihood ratio test, assuming that twice the difference between the log-likelihoods of both models (2Δℓ) follow a χ2 distribution with one degree of freedom (Huelsenbeck and Crandall 1997). Tajima’s relative rate tests (Tajima 1993) were conducted using in-house scripts. A. thaliana genes were classified into different eukaryotic orthologous groups (KOG) categories using the eggNOG database v4.5.1 (Huerta-Cepas et al. 2015). Data for the 18 accessions of A. thaliana were obtained from the 1000 genomes project (Gan et al. 2011). A. thaliana gene expression data were obtained from Schmid et al. (2005) and processed as in Alvarez-Ponce and Fares (2012). All our alignments and scripts are available upon request.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.