Efficient purging of deleterious mutations in plants with haploid selfing.

In diploid organisms, selfing reduces the efficiency of selection in removing deleterious mutations from a population. This need not be the case for all organisms. Some plants, for example, undergo an extreme form of selfing known as intragametophytic selfing, which immediately exposes all recessive deleterious mutations in a parental genome to selective purging. Here, we ask how effectively deleterious mutations are removed from such plants. Specifically, we study the extent to which deleterious mutations accumulate in a predominantly selfing and a predominantly outcrossing pair of moss species, using genome-wide transcriptome data. We find that the selfing species purge significantly more nonsynonymous mutations, as well as a greater proportion of radical amino acid changes which alter physicochemical properties of amino acids. Moreover, their purging of deleterious mutation is especially strong in conserved regions of protein-coding genes. Our observations show that selfing need not impede but can even accelerate the removal of deleterious mutations, and do so on a genome-wide scale.

Introduction

Heritable genetic variation, the substrate of natural selection (Maynard Smith and Szathmáry 1995), is ultimately caused by mutations. Because few mutations are beneficial and most are deleterious when they first arise, populations of organisms carry a genetic load of deleterious mutations in their genomes. Whether natural selection can remove such deleterious mutations from a population depends on various factors, such as how strongly deleterious the mutations are, how large the population is, whether an organism is diploid or haploid, and whether it reproduces sexually or asexually (Glémin 2007; Otto and Gerstein 2008; Charlesworth and Willis 2009; Glémin and Galtier 2012). Among sexually reproducing organisms, the mating system—whether organisms reproduce primarily through selfing or outcrossing—is a key factor that determines how effectively deleterious mutations can be removed.

On the one hand, selfing reduces a population’s effective size Ne, and thus selection’s ability to remove weakly deleterious mutations (Charlesworth and Wright 2001; Wright et al. 2008; Wright et al. 2013). In particular, in organisms with a diploid-dominant life phase, selfing increases genome-wide homozygosity, which in turn reduces the effective number of independently sampled gametes in a population, an effect that can cut the effective population size in half (Nordborg 2000). In addition, high levels of homozygosity cause low effective recombination rates (Charlesworth and Wright 2001) that further decrease Ne and promote the accumulation of deleterious mutations (Felsenstein 1974). Finally, the very population structure of selfers, which is characterized by small subpopulations in which founder events can be frequent, can also reduce the effective population size (Charlesworth and Wright 2001).

On the other hand, selfing can facilitate the purging of recessive deleterious mutations (Glémin 2007; Glémin and Galtier 2012). Such mutations can be completely masked in diploids if they are paired with a wild-type allele in a heterozygote. Selfing increases the incidence of homozygotes where two recessive deleterious alleles are paired (Charlesworth and Willis 2009). It can thus expose deleterious mutations and facilitate their purging.

Both theoretical and experimental evidence suggest that the Ne-reducing effects of selfing are usually stronger than the enhanced exposure of deleterious mutations caused by selfing (Bustamante et al. 2002; Slotte et al. 2010, 2013; Qiu et al. 2011; Ness et al. 2012; Hazzouri et al. 2013). In other words, selfers with a diploid-dominant life cycle should overall remove deleterious mutations less well than outcrossers (Charlesworth and Wright 2001; Glémin 2007).

However, this may not be the case in plants with a haploid-dominant life cycle or with a free-living haploid phase for two reasons. First, theory suggests that effective population size reduction caused by selfing should be less pronounced in haploid-dominant than in diploid-dominant plants because their life cycle poorly fit the assumptions underlying Wright–Fisher populations (Stenøien and Såstad 2001; Bengtsson and Cronberg 2009). More specifically, in bryophytes, the automatic 2-fold reduction of Ne caused by selfing does not operate. Ne is thus primarily reduced by a decrease in effective recombination rates (Hedrick 1987a, 1987b).

In addition to that, plants with an independent haploid phase or with a haploid-dominant life cycle can undergo “intragametophytic selfing,” which can be more effective in purging deleterious mutations than regular selfing in diploid-dominant plants. Intragametophytic selfing occurs between genetically identical male and female gametes produced by the same gametophyte (bisexual gametophyte) (fig. 1; Shaw 2000). It leads to a complete loss of heterozygosity in merely one generation (Shaw 2000; Barner et al. 2011; Billiard et al. 2012). In the resulting diploids, all deleterious mutations are exposed to selection. In this way, intragametophytic selfing can facilitate the purging of deleterious mutations (Eppley et al. 2007). This extreme form of selfing is not possible in plants with unisexual gametophytes, where female and male gametes are always produced on genetically separate haploid individuals. Such plants must undergo outcrossing or intergametophytic selfing—a less extreme form of selfing among siblings produced by the same diploid (an analogous situation to selfing in a hermaphroditic flowering plant, see fig. 1), which reduces heterozygosity only by a factor one-half per generation (Charlesworth and Wright 2001).

F

Schematic representation of the mating system of selfer and outcrosser organisms with a haploid-dominant life cycle. Solid black arrows refer to meiotic events and the dashed horizontal line separates the haploid and diploid phases. Thickness of the gray arrows indicates the relative frequency of selfing and outcrossing events in selfers and outcrossers.

Taken together, the available body of theory (Hedrick 1987a, 1987b; Holsinger 1987) suggests that haploid-dominant plants with frequent intragametophytic selfing might purge their genetic load more effectively than outcrossers. This prediction is supported by some indirect experimental evidence, for example, less severe inbreeding depression in plants capable of intragametophytic selfing (Klekowski 1973, 1982; Lloyd 1974; Taylor et al. 2007; Eppley et al. 2007).

However, the proportion of the genome subject to purging via homozygous exposure is not known. Yet, this information is crucial to understand the extent at which intragametophytic selfing can affect the accumulation of deleterious mutations on a genome-wide scale. If we assume that only diploid-specific genes are subject to purging, then the genome-wide effect of intragametophytic selfing on the accumulation of deleterious mutations should be minor (Szövényi et al. 2011). This is because only very few genes are diploid-specific (2–3%; Szövényi et al. 2011, 2013) and because genes also expressed in the haploid phase may not be affected by purging due to effective haploid selection (Klekowski 1973, 1982; Lloyd 1974; McCauley et al. 1985; Taylor et al. 2007; Eppley et al. 2007; Billiard et al. 2012). Alternatively, if genes expressed in both phases can also be affected by purging via homozygous exposure, intragametophytic selfing is expected to have a genome-wide effect because the majority of the genes is shared by the two phases (97%; Szövényi et al. 2011, 2013). Therefore, no experimental evidence specifically addresses the question: Do intragametophytic selfers accumulate a lower number of deleterious mutations on a genome-wide scale?

To answer this question, we here use mosses as model systems for plants with a haploid-dominant life cycle capable of intragametophytic selfing. We compare the genome-wide accumulation of deleterious mutations in protein-coding genes between two pairs of moss species that either predominantly self or outcross. We find that selfing mosses accumulate fewer deleterious mutations, and especially so in conserved positions of protein-coding genes. Thus, intragametophytic selfing in haploid-dominant plants can indeed facilitate effective purging of deleterious mutations.

Materials and Methods

Study Species

We use transcriptomes of two pairs of moss species to compare the number of deleterious mutations accumulated between species that predominantly self or outcross. The first pair comprises Physcomitrella patens and its relative Funaria hygrometrica. Both P. patens and F. hygrometrica have haploid-dominant life cycles, and their haploid individuals with combined sexes primarily undergo intragametophytic selfing (Shaw 1991; Eppley et al. 2007; Perroud et al. 2011). Because they undergo more selfing, and a more extreme form of selfing than the other species pairs, we refer to them as the “selfers.” Both species are genetically and functionally haploid and neither of them has an allo- or autopolyploid origin (Rensing et al. 2007; Shaw et al. 2010; Liu et al. 2012).

Our pair of outcrossing organisms (hereinafter referred to as “outcrossers”) includes two peat moss species Sphagnum subsecundum and Sphagnum cribrosum. Their haploid life stages can be either females or males, that is, individual plants can bear either male or female gametangia, but never both. They can outcross or undergo intergametophytic selfing, that is, fertilization between products of a single meiosis, comparable to selfing in diploid plants (Shaw 2000). Although they may therefore experience some amount of selfing (Eppley et al. 2007; Szövényi et al. 2009), we refer to them as outcrossers to distinguish them from the more extreme selfers P. patens and F. hygrometrica. Importantly, the divergence times between the selfer and outcrosser species pairs are similar (McDaniel et al. 2010; Shaw et al. 2010; Liu et al. 2012; Szövényi et al. 2013).

Available evidence suggests that the selfer and outcrosser species pairs have comparable range size and population structure. Both the selfer and the outcrosser pairs consist of a broadly and a more narrowly distributed species. The selfer F. hygrometrica is world-wide distributed, whereas P. patens occurs primarily in Europe. Similarly, the outcrosser S. subsecundum is widely distributed across the Northern Hemisphere whereas S. cribrosum only occurs in the United States (Anderson et al. 2009). Both the selfer and the outcrosser species pairs have well-interconnected populations (via long-distance spore dispersal) which is evidenced by weak among-population genetic differentiation both at the local and global scales (Shaw 1991; Shaw et al. 2008; Szövényi et al. 2008). Within-population genetic diversity is considerably higher in the outcrossers indicating greater effective population size of local populations (Shaw 1991; Shaw et al. 2008; Szövényi et al. 2008). Finally, both species pairs are effective colonizers with small, temporary populations (Shaw 1991; Shaw et al. 2008; Szövényi et al. 2008; Anderson et al. 2009). Therefore, range size and population structures of the selfer and the outcrosser species pairs are comparable, rendering mating system differences the primary candidate factor to affect mutation purging.

We focus our analysis on among-species divergence data because molecular traces of selfing are predicted to be more pronounced for among-species divergence data than for within-population polymorphism data when mating systems have been stable for a long time (Glémin 2007). Phylogenetic evidence suggests that the breeding systems of our selfer and outcrosser species pairs have indeed been stable, so the mating system had enough time to affect the accumulation of mutations (Shaw et al. 2010; Liu et al. 2012).

Sequence Data and Read Filtering

As a proxy of the incidence of mutation purging, we use the accumulation of likely deleterious mutations between selfer and outcrosser species pairs. To estimate the incidence of such mutations, we used high-throughput sequencing to generate two genome-wide transcriptome data sets for this study (S. subsecundum and S. cribrosum) and employed another such data set from our previous work (F. hygromerica) (Szövényi et al. 2013). Briefly, for F. hygrometrica we had collected samples of three haploid and three diploid developmental stages, extracted RNA, and performed sequencing on an Illumina GAIIx machine that generated 173,542,944 raw reads. We then discarded all reads containing low quality base calls (Q < 20 and one or more undefined nucleotides [N], FASTQ Quality Trimmer, http://hannonlab.cshl.edu/fastx_toolkit/ , last accessed May 20, 2014) and/or showing adaptor contamination (Tagdust, false discovery rate threshold of 0.01; Lassmann et al. 2009). We used 133,851,302 quality-filtered reads in the final assembly. Detailed description of data collection, read filtering and assembly are described in Szövényi et al. (2013). Sequence data are available in the ArrayExpress database (www.ebi.ac.uk/arrayexpress, last accessed May 16, 2014) under accession number E-MTAB-1664.

For S. subsecundum and S. cribrosum, we used the Spectrum Plant Total RNA kit (Sigma) to extract total RNA from 100 mg of gametophytic tissue (capitulum) obtained from 6-month-old cultivated plants originally collected in Swanville, Waldo County, Maine, USA. We note that at the time of collection no sporophytic tissues were available and thus we extracted RNA from purely gametophytic material. Because only few genes are expected to show sporophyte-specific expression, differential tissue sampling in the selfer and outcrosser species pairs is unlikely to significantly bias our analyses (Szövényi et al. 2011, 2013). We subjected Poly-A selected RNA to paired-end sequencing (insert size: 300 bp) on one full lane and half a lane of an Illumina Hiseq2000 sequencer for S. subsecundum and S. cribrosum, respectively. The sequencing runs provided 163,776,364 and 94,818,996 raw paired end reads for S. subsecundum and S. cribrosum, respectively. Prior to assembly, we filtered out duplicate reads using an in-house Python script, trimmed the 3′-end of each read in a pair based on quality using the FASTX-Toolkit (FASTQ Quality Trimmer, http://hannonlab.cshl.edu/fastx_toolkit/, last accessed May 16, 2014) with a quality score threshold of 25, clipped off Illumina sequencing adapters, when present, with the FASTX-Toolkit (FASTQ/A Clipper, http://hannonlab.cshl.edu/fastx_toolkit/, last accessed May 16, 2014), and finally kept only paired-end reads for which both sequences were longer than 40 bp after trimming and clipping. After these steps, 117,836,947 and 81,917,656 cleaned paired end reads were left for S. subsecundum and S. cribrosum, respectively. Raw sequence data obtained are available in the ArrayExpress database (www.ebi.ac.uk/arrayexpress, last accessed May 16, 2014) under accession number E-MTAB-2482.

For the fourth species, P. patens, we retrieved high quality genomic data (Rensing et al. 2008) for transcript sequences and their corresponding protein translations from a public database (http://www.phytozome.net/version8.0, last accessed October 8, 2013).

Transcriptome Assembly

We assembled the S. subsecundum, S. cribrosum, and F. hygrometrica RNA-seq data into virtual transcripts using the assembler Trinity with default options (Grabherr et al. 2011), which has high sensitivity and leads to the most contiguous transcriptome assembly, for example, the highest number of full-length transcripts for nonmodel species. In S. subsecundum and S. cribrosum, we assembled a total of 117,836,947 and 81,917,656 cleaned paired end reads into 541,929 and 508,499 contigs, respectively (371,788 and 261,128 “components” according to Grabherr et al. 2011). In F. hygrometrica, we obtained 133,851,302 cleaned single end reads and assembled them into 192,665 transcripts (102,856 “components” according to Grabherr et al. 2011).

Protein Translation, Orthology, and Pairwise Alignment of Orthologs

In order to predict coding frame and to generate protein translations, we compared all contigs of the assembly against the National Center for Biotechnology Information (NCBI) nr database plant division (ftp://ftp.ncbi.nlm.nih.gov/blast/db/nr.xx.tar.gz , last accessed May 10, 2013) with the aid of BLASTX (Altschul et al. 1997). We filtered the Blast (ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.2.25+-ia32-linux.tar.gz, last accessed May 10, 2013; Altschul et al. 1997) output and kept only transcripts with BLAST hits showing at least 30% identity over at least 150 amino acids. Using these filtering criteria, 21,101 (S. subsecundum), 16,991 (S. cribrosum), and 34,119 (F. hygrometrica) transcripts had a valid match against the database. After that we provided each transcript and its best hit protein as input to the program GeneWise2 (http://www.ebi.ac.uk/∼birney/wise2/Wise2.4.0, last accessed March 8, 2013; Birney et al. 2004), and used the best scoring protein-to-cDNA alignment to predict the encoding gene’s ORF (Open Reading Frame) and generate a protein translation. We also discarded all transcripts with internal stop codons. The various filtering steps and the homology-guided translation ensured a high-confidence transcript set (we discarded transcripts without homology to the database) and helped exclude artifactual or contaminant sequences that are expected to occur in de novo transcriptome assemblies (Vijay et al. 2013).

We next identified orthologous gene pairs from the sets of proteins thus identified in each species, and did so separately for the outcrosser (S. subsecundum vs. S. cribrosum) and for the selfer species pair (F. hygrometrica vs. P. patens). In order to establish orthology, we first kept the transcript with the longest protein translation for each putative gene (referred to as a “component” in the software Trinity). We then used the reciprocal best hits approach (Bork et al. 1998) to establish orthologous protein pairs, keeping only pairs showing 30% identity over at least 150 amino acids in a BLASTP search (ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.2.25+-ia32-linux.tar.gz, last accessed May 10, 2013; Altschul et al. 1997).

To generate DNA alignments of the coding sequence of orthologous proteins identified, we used Muscle (Edgar 2004) with default parameters. After that we used PAL2NAL (Suyama et al. 2006) to map protein alignments onto nucleotide alignments. We used the “remove gaps” option of PAL2NAL in order to generate codon alignments without gaps.

Estimating the Accumulation of Deleterious Mutations

To determine whether deleterious mutations accumulate at different rates in selfers compared with outcrossers, we used three complementary approaches.

First, we followed well-established precedence (Wright and Andolfatto 2008; Gossmann et al. 2010; Slotte et al. 2011; Yang and Gaut 2011) by assuming that the majority of nonsynonymous mutations is deleterious, and estimated their incidence relative to silent mutations. Specifically, we used our pairwise alignments to obtain estimates for the number of nonsynonymous substitutions per nonsynonymous sites (dN), as well as the number of synonymous substitutions per synonymous (dS) sites, and computed their ratio dN/dS using the KaKs_Calculator (Zhang et al. 2006) with the YN (Yang and Nielsen) model of sequence evolution (Yang and Nielsen 2000). We filtered the pairwise alignments by excluding entries with dS > 2, dS = 0, and dN/dS > 2. We used these thresholds because alignments with greater dS values are likely to become saturated with synonymous substitutions (Axelsson et al. 2008; Mank et al. 2010). Furthermore, visual inspection indicated that alignments with dS > 2 or dN/dS > 2 are frequently unreliable and thus cannot be trusted. Subsequently, we asked whether the dN/dS ratios differed between selfer and outcrosser species pairs. Furthermore, we also tested whether any such difference persists after excluding genes potentially under positive selection. For the latter analysis, we kept only gene pairs with dN/dS < 1, following the method of Subramanian (2013). With this approach we tried to distinguish relaxed from positive selection, because both can accelerate the accumulation of nonsynonymous mutations and thus lead to elevated dN/dS ratios.

Different amino acids differ in their physicochemical properties, such as size, charge, or aromaticity. Conservative amino acid changes, which alter an amino acid into one with similar physicochemical properties, are usually less detrimental to protein function than radical changes, which alter these properties (Zhang 2000; Eyre-Walker et al. 2002; Hughes and Friedman 2009). The ratio of the number of radical to conservative amino acid changes that a protein has experienced can thus be used as a proxy for the incidence of deleterious mutations (Zhang 2000; Eyre-Walker et al. 2002; Smith 2003). We used a complex classification of amino acids established by Miyata et al. (1979), which group amino acids into acidic, basic, neutral small, and neutral large categories. We treated all amino acid replacements that change an amino acid into one of a different category as radical, whereas we classified category-preserving changes as conservative. This classification has performed well in multiple studies and provides a reliable way to describe how strongly an amino acid change affects amino acid properties (Hanada et al. 2007; Wernegreen 2011). We estimated the proportion of radical (Dr) and conservative (Dc) amino acid changes using the hon-new program (Zhang 2000), taking into account the transition/transversion bias for each pair of genes, as estimated by PAML (pairwise option, model: F3×4, kappa). Before calculating these statistics, we again filtered alignments based on dN/dS ratios, as explained above (0 < dS < 2, dN/dS < 2). Furthermore, we discarded all alignments with no amino acid change (Dr = Dc = 0) or with no radical change (Dr = 0). We also repeated the analysis by excluding genes that may be subject to positive selection, that is, by retaining only alignments with dN/dS < 1 (Subramanian 2013). We determined significant differences in dN/dS and Dr/Dc ratios among outcrossers and selfers with a Mann–Whitney U test (Sokal and Rohlf 2012).

Finally, we also used a third approach to estimate the accumulation of deleterious mutations. It relies on the observation that the lower effective population size (Ne) potentially caused by selfing will allow slightly deleterious mutations accumulate more rapidly, and thus lead to an elevated dN/dS ratio between-species. Theory predicts that the effect of Ne is more pronounced on sites that are under more intense purifying selection, where the average selection coefficient of deleterious mutations is greater (Subramanian 2013). This prediction is a consequence of the theoretical relationship between ω (dN/dS), Ne, and the selection coefficient (s), ω = S/(1 − e−) where S = 4Ne × s (Kimura 1983; Nielsen and Yang 2003; Subramanian 2013).

Because the rate at which mutations accumulate depends on the intensity of selection, the number of accumulated deleterious mutations is expected to be lower in conserved than in nonconserved parts of proteins (sconserved > snonconserved). Furthermore, the relationship among Ne, s, and ω predicts that a unit increase in purifying selection (s) will result in a steeper decrease in ω in the species pair with the larger effective population size (Ne). That is, given that Ne (selfer) considerably differs from Ne (outcrosser) the difference in the number of accumulated deleterious mutations between outcrossers and selfers should be more pronounced in conserved parts of proteins, which are under more intense purifying selection. Furthermore, differences in the dN/dS and Dr/Dc ratios between conserved and nonconserved regions within genes should also differ between outcrossers and selfers. The effect of purging via homozygous exposure should be also more pronounced in conserved than in nonconserved domains of proteins, because deleterious mutations with large effect can be more effectively purged (Wang et al. 1999). For example, if purging is more efficient in selfers, we expect to see not only lower dN/dS and Dr/Dc ratios in selfers but also the differences to outcrossers should be more pronounced in conserved protein regions. Furthermore, dN/dS and Dr/Dc ratios between conserved and nonconserved regions within genes should differ more markedly in selfers if they experience more effective purging.

To perform this analysis, we searched protein translations of assembled transcripts against the conserved domain database (Marchler-Bauer and Bryant 2004; Marchler-Bauer et al. 2013) to identify highly conserved regions of the genes in our data set. We used the default search parameters and kept only the best match for each protein. After that, we merged multiple conserved regions per protein into a contiguous stretch of amino acids, likewise for nonconserved regions. To ensure that our analysis is based on sufficiently many amino acid changes, we only used those pairwise alignments where each protein contained a conserved and nonconserved region of at least 100 amino acids each.

Furthermore, as in the previous analyses, we discarded uninformative or unreliable protein pairs (dS = 0, dS > 2, dN/dS > 2). Finally, we used the Ka/Ks estimator (Zhang et al. 2006) and the hon-new program (Zhang 2000) to obtain the ratios dN/dS and Dr/Dc, and did so separately for the conserved and nonconserved parts of each gene. We determined significant differences in dN/dS and Dr/Dc ratios among outcrossers and selfers for conserved and nonconserved regions separately with a Mann–Whitney U test (Sokal and Rohlf 2012). One limitation of this test is that we treat conserved and nonconserved regions as independent entities. This is because no meaningful pairing exists when analyzing the total data set. However, dependency of conserved and nonconserved regions was fully accounted for (Wilcoxon matched pairs tests, see below) when analyzing the one-to-one ortholog data set. We also used Wilcoxon matched pairs tests (Sokal and Rohlf 2012) to determine whether conserved and nonconserved regions (the matched pairs) within each gene showed significantly different dN/dS and Dr/Dc ratios for each species pair separately. We performed all statistical analyses in R (R Development Core Team 2011).

Comparing functionally different gene sets could obscure the effect of mating system on the evolutionary rate of genes. To account for this bias, we identified proteins that showed one-to-one orthology across all four of our study species because orthologous proteins tend to be more similar in function than paralogs (Altenhoff et al. 2012). We identified orthologs using the algorithm provided in ProteinOrtho (Lechner et al. 2011) and repeated all the above analyses on this data set. To establish orthology across the four species gene set, we executed ProteinOrtho using all predicted proteins of S. subsecundum, S. cribrosum, F. hygrometrica, and P. patens (http://www.phytozome.net/version 8.0, last accessed October 8, 2013) with the following parameters: e value threshold of 10−4, 30% identity, adaptive best alignments similarity of f = 0.95, algebraic connectivity > 0.1. Subsequently, we kept only those co-orthologs containing one single protein per species in order to reduce the risk of including paralogous relationships. Finally, we used the four-way orthology relationship to identify corresponding pairwise orthologous genes (S. subsecundum vs. S. cribrosum; F. hygrometrica vs. P. patens).

Results

Selfers Purge Deleterious Mutations Better than Outcrossers

Our analysis of complete quality-filtered transcriptomes (see Materials and Methods) identified 10,038 orthologous gene pairs in the selfer species (F. hygrometrica and P. patens) and 6,089 orthologous gene pairs in the outcrosser species pair (S. subsecundum and S. cribrosum) (table 1).

Table 1

Accumulation of Deleterious Mutations among the Selfer and the Outcrosser Species Pairs across the Full Coding Region of Proteins

Type of Analysis	Mating System	Species Pair	Number of Pairwise Alignments	Length of Alignment in Base Pairs, Mean (median)	Lower–Upper Quartiles	dN or Dr, Mean (median)	Lower–Upper Quartiles	dS or Dc, Mean (median)	Lower–Upper Quartiles	dN/dS or Dr/Dc, Mean (median)	Lower–Upper Quartiles	Significance (selfer vs. outcrosser)
Synonymous–nonsynonymous mutations (dN/dS)
	Selfer (S)	Funaria hygrometrica versus Physcomitrella patens
			10,038	1,065 (879)	624–1,293	0.0539 (0.0359)	0.0202–0.0601	0.2926 (0.2649)	0.2218–0.3270	0.1957 (0.1339)	0.0742–0.2174	O (P < 0.00001)
		dN/dS < 1	10,038	1,065 (879)	624–1,293	0.0539 (0.0359)	0.0202–0.0601	0.2926 (0.2649)	0.2218–0.3270	0.1957 (0.1339)	0.0742–0.2174	O (P < 0.00001)
	Outcrosser (O)	Sphagnum subsecundum versus S. macrophyllum
			6,089	1,195 (963)	675–1,452	0.0727 (0.0226)	0.0096–0.0571	0.2455 (0.0516)	0.0291–0.1214	0.4222 (0.4163)	0.2376–0.6030
		dN/dS < 1	5,780	1,213 (975)	681–1,470	0.0703 (0.0220)	0.0095–0.0577	0.2540 (0.0537)	0.0307–0.1234	0.3757 (0.3866)	0.2334–0.5657
Radical-conservative amino acid changes (Dr/Dc)
	Selfer (S)	F. hygrometrica versus P. patens
			9,913	1,075 (885)	624–1,308	0.0421 (0.0360)	0.0200–0.0580	0.0740 (0.0680)	0.0430–0.0940	0.5590 (0.5526)	0.4000–0.7391	O (P < 0.00001)
		dN/dS < 1	9,900	1,082 (891)	630–1,317	0.0419 (0.0360)	0.0200–0.0580	0.0737 (0.0680)	0.0430–0.0940	0.5590 (0.5526)	0.4000–0.7391	O (P < 0.00001)
	Outcrosser (O)	S. subsecundum versus S. macrophyllum
			5,643	1,084 (876)	633–1,305	0.0539 (0.0220)	0.0100–0.0610	0.0810 (0.0300)	0.0150–0.1030	0.7576 (0.6364)	0.4930–0.8889
		dN/dS < 1	5,338	1,233 (987)	690–1,500	0.0530 (0.0220)	0.0100–0.0580	0.0806 (0.0290)	0.0140–0.1000	0.7475 (0.6309)	0.4880–0.8468

Note.—dS, proportion of synonymous mutations per synonymous sites; dN, proportion of nonsynonymous mutations per nonsynonymous sites; Dc, proportion of conserved amino acid replacements; Dr, proportion of radical amino acid replacements; Significance, statistical significance according to a Mann–Whitney U test. We conducted two separate tests: One for the unfiltered and one for the filtered (dN/dS < 1) data sets.

Overall, we found that the mean dN/dS ratios among orthologous genes in both species pairs fell into the interval 0.1–0.4, a range that is typical for such genome-wide analyses (Slotte et al. 2011; Yang and Gaut 2011). Importantly, the dN/dS ratio in the outcrosser species pair was twice as high than in the selfer species pair, and this difference was statistically significant, according to a Mann–Whitney U test (fig. 2 and table 1; zdf=1 = 53.48, P < 0.00001). After removing gene pairs that may be subject to positive selection (dN/dS > 1), the dN/dS ratio was still much higher in the outcrosser species (table 1; zdf=1 = 50.01, P < 0.00001). Repeating the above analyses on the data set containing only genes that are one-to-one orthologous across our four study species led to the very same conclusions (supplementary table S1, Supplementary Material online). In sum, deleterious mutations are purged more efficiently from the genomes of the selfer species.

F

Ratio of nonsynonymous mutations per nonsynonymous sites to synonymous mutations per synonymous sites (dN/dS, mean, and SE) in the selfer and in the outcrosser species pairs. dN/dS ratios are significantly greater in the outcrosser than in the selfer species pair (P < 0.00001).

Selfers Purge Radical Amino Acid Changes More Efficiently than Outcrossers

Many radical changes are deleterious and eliminated by natural selection (Zhang 2000). In a complementary analysis which focused on the ratio of radical to conservative amino acid changes Dr/Dc, we found that the Dr/Dc ratio was about 36% greater in the outcrosser than in the selfer species pair (fig. 3 and table 1; zdf=1 = 25.62, P < 0.00001). We repeated this analysis by excluding those genes that may be subject to positive selection (dN/dS > 1) which led to the same outcome (table 1, selfer vs. outcrosser zdf=1 = 23.59, P < 0.00001). Finally, repeating these analyses on one-to-one orthologs led to the very same conclusions (supplementary table S1, Supplementary Material online). Analysis of radical and conserved amino acid changes between selfer and outcrosser species pairs thus shows that selfers are more efficient in purging deleterious mutations than outcrossers, which is consistent with our observations based on dN/dS ratios.

F

Ratio of radical to conserved amino acid changes (Dr/Dc, mean, and SE) in the selfer and in the outcrosser species pairs. Dr/Dc ratios are highly significantly greater in the outcrosser than in the selfer species pair (P < 0.00001).

Comparison of Conserved and Nonconserved Stretches of Genes

In our next analysis, we contrasted dN/dS or Dr/Dc ratios in conserved and nonconserved regions of genes (see Materials and Methods). We first investigated whether the difference in dN/dS values between conserved and nonconserved regions within each gene differs more considerably in selfers than in outcrossers, which would imply more effective purging in selfers. We found that dN/dS ratios were smaller in conserved than in nonconserved regions of genes in both the selfer and the outcrosser species pairs. Importantly, in selfing species, this difference was more pronounced (table 2 and fig. 4). Specifically, mean dN/dS ratios were on average 116% greater in nonconserved regions of selfers (zS[conserved vs. nonconserved] = 40.03, P < 0.00001), and only 73% greater in outcrossers (zO[conserved vs. nonconserved] = 21.97, P < 0.00001), which is also evident from the z scores of the Wilcoxon matched pairs test (zS[conserved vs. nonconserved] = 40.03, zO[conserved vs. nonconserved] = 22.33; table 2 and fig. 4). We repeated this analysis after exclusion of genes that may have been subject to positive selection (dN/dS > 1), which leads to the same conclusion (table 2; zS[conserved vs. nonconserved] = 40.03, P < 0.00001; zO[conserved vs. nonconserved] = 22.17, P < 0.00001), as did the same analysis for one-to-one orthologous genes (supplementary table S2, Supplementary Material online).

Table 2

Accumulation of Deleterious Mutations (dN/dS) between the Selfer and the Outcrosser Species Pairs in Conserved Domain and Nonconserved Domain Stretches of Proteins

Sequence Stretch	Mating System	Species Pair	Number of Pairwise Alignments	Length of Alignment in Base Pairs, Mean (median)	Lower– Upper Quartiles	dN, Mean (median)	Lower–Upper Quartiles	dS, Mean (median)	Lower–Upper Quartiles	dN/dS, Mean (median)	Lower–Upper Quartiles	Significance (selfer vs. outcrosser; conserved and nonconserved domains separately tested)
Conserved domains
	Selfer (S)	Funaria hygrometrica versus Physcomitrella patens
			2,995	768 (639)	459–930	0.0321 (0.0219)	0.0124–0.0362	0.2835 (0.2504)	0.2123–0.3120	0.1260 (0.0814)	0.0454–0.1496	O (P < 0.00001)
		dN/dS < 1	2,995	768 (639)	459–930	0.0321 (0.0219)	0.0124–0.0362	0.2835 (0.2504)	0.2123–0.3120	0.1260 (0.0814)	0.0454–0.1496	O (P < 0.00001)
	Outcrosser (O)	Sphagnum subsecundum versus S. macrophyllum
			1,699	786 (666)	465–969	0.0420 (0.0079)	0.0041–0.0199	0.1679 (0.0361)	0.0246–0.0702	0.2651 (0.1904)	0.1036–0.3575
		dN/dS < 1	1,665	804 (678)	468–987	0.0309 (0.0078)	0.0040–0.0184	0.0922 (0.0361)	0.0246–0.0702	0.2448 (0.1889)	0.1020–0.3527
Nonconserved domains
	Selfer (S)	F. hygrometrica versus P. patens
			2,995	1,376 (765)	456–1500	0.0661 (0.0513)	0.02883–0.0843	0.2761 (0.1559)	0.2109–0.2991	0.2629 (0.1975)	0.1092–0.3296	O (P < 0.001)
		dN/dS < 1	2,995	1,376 (765)	456–1500	0.0661 (0.0513)	0.02883–0.0843	0.2761 (0.2560)	0.2109–0.2991	0.2629 (0.1975)	0.1092–0.3296	O (P < 0.001)
	Outcrosser (O)	S. subsecundum versus S. macrophyllum
			1,699	970 (720)	456–1152	0.0760 (0.0204)	0.0104–0.0591	0.1860 (0.0467)	0.0283–0.1273	0.4523 (0.4170)	0.2515–0.6220
		dN/dS < 1	1,665	1,016 (736)	465–1215	0.0638 (0.0202)	0.0099–0.0590	0.1153 (0.0461)	0.0282–0.1274	0.4034 (0.4083)	0.2460–0.6183

Note.—dS, proportion of synonymous mutations per synonymous sites; dN, proportion of nonsynonymous mutations per nonsynonymous sites; Significance, statistical significance according to a Kruskal–Wallis ANOVA. We conducted two separate tests: One for the unfiltered and one for the filtered (dN/dS < 1) data sets. We show results of the two tests for conserved domain and nonconserved domain stretches of proteins separately.

F

Ratio of nonsynonymous mutations per nonsynonymous sites to synonymous mutations per synonymous sites (dN/dS, mean, and SE) for conserved and nonconserved domains of genes in the selfer and in the outcrosser species pairs. dN/dS ratios in conserved and nonconserved domains of genes are highly significantly different both within the selfer and the outcrosser species pairs studied (P < 0.00001). However, the difference between conserved and nonconserved domains is less pronounced in the outcrosser comparison. Finally, dN/dS ratios are significantly (P < 0.00001) greater in the outcrosser than in the selfer species pair both in conserved and in nonconserved domains with a greater effect in conserved domains.

Because the ratio of radical to conservative amino acid changes (Dr/Dc) also reflects mutational purging, we reasoned that it should show a pattern similar to dN/dS in conserved and nonconserved gene regions (table 3 and fig. 5). Indeed, the proportion of radical amino acid changes was significantly smaller (15%) in conserved regions for selfers (zS[conserved vs. nonconserved] = 18.38, P < 0.00001). This was also true for outcrossers but to a significantly lesser extent (10%, zO[conserved vs. nonconserved] = 5.51, P < 0.00001). This conclusion did not change when we excluded genes that may have been subject to positive selection (dN/dS > 1; table 3; zS[conserved vs. nonconserved] = 18.38, P < 0.00001; zO[conserved vs. nonconserved] = 5.36, P < 0.00001) and when we restricted our analysis for one-to-one orthologous genes (supplementary table S3, Supplementary Material online). Altogether, this means that conserved regions accumulate fewer radical amino acid changes in the selfers than in the outcrossers. This supports the hypothesis that purging of deleterious mutations through homozygous exposure is more efficient in selfers than in outcrossers.

Table 3

Accumulation of Deleterious Mutations (Dr/Dc) in the Selfer and the Outcrosser Species Pair in Conserved Domain and Nonconserved Domain Stretches of Proteins

Sequence Stretch	Mating System	Species Pair	Number of Pairwise Alignments	Length of Alignment in Base Pairs, Mean (median)	Lower– Upper Quartiles	Dr, Mean (median)	Lower– Upper Quartiles	Dc, Mean (median)	Lower–Upper Quartiles	Dr/Dc, Mean (median)	Lower– Upper Quartiles	Significance (selfer vs. outcrosser; conserved and nonconserved domains separately tested)
Conserved domains
	Selfer (S)	Funaria hygrometrica versus Physcomitrella patens
			2,746	768 (639)	459–930	0.0254 (0.0160)	0.0070–0.0280	0.0534 (0.0390)	0.0240–0.0570	0.5304 (0.4130)	0.2389–0.6027	O (P < 0.001)
		dN/dS < 1	2,746	768 (639)	459–930	0.0254 (0.0160)	0.0070–0.0280	0.0534 (0.0390)	0.0240–0.0570	0.5304 (0.4130)	0.2389–0.6027	O (P < 0.001)
	Outcrosser (O)	Sphagnum subsecundum versus S. macrophyllum
			1,203	858 (750)	504-1,050	0.0389 (0.0090)	0.0050–0.0310	0.0637 (0.0170)	0.0090–0.0540	0.6677 (0.5455)	0.3514–0.8462
		dN/dS < 1	1,168	866 (759)	504–1,065	0.0388 (0.0090)	0.0050–0.0290	0.0638 (0.0160)	0.0090–0.0530	0.6605 (0.5455)	0.3571–0.8358
Nonconserved domains
	Selfer (S)	F. hygrometrica versus P. patens
			2,746	1,374 (765)	456–1,500	0.0551 (0.0400)	0.0190–0.0680	0.0893 (0.0730)	0.0490–0.1100	0.6160 (0.5294)	0.3871–0.7086	O (P < 0.00001)
		dN/dS < 1	2,746	1,374 (765)	456–1,500	0.0551 (0.0400)	0.0190–0.0680	0.0894 (0.0730)	0.0490–0.1100	0.6160 (0.5294)	0.3871–0.7086	O (P < 0.00001)
	Outcrosser (O)	S. subsecundum versus S. macrophyllum
			1,203	1,047 (783)	486–1,263	0.0680 (0.0220)	0.0100–0.0690	0.1004 (0.0300)	0.0160–0.1150	0.7383 (0.6667)	0.4750–0.8726
		dN/dS < 1	1,168	1,061 (792)	492–1,281	0.0681 (0.0210)	0.0100–0.0680	0.1013 (0.0320)	0.0160–0.1150	0.7285 (0.6667)	0.4737–0.8587

Note.—Dc, proportion of conservative amino acid replacements; Dr, proportion of radical amino acid replacements; Significance, statistical significance according to a Kruskal–Wallis ANOVA. We conducted two separate tests: One for the unfiltered and one for the filtered (dN/dS < 1) data sets. We show results of the two tests for conserved domain and nonconserved domain stretches of proteins separately.

F

Ratio of radical to conserved amino acid changes (Dr/Dc, mean, and SE) for conserved and nonconserved domains of genes in the selfer and in the outcrosser species pairs. Dr/Dc ratios in conserved and nonconserved domains of genes are significantly different both within the selfer (P < 0.00001) and in the outcrosser species pairs (P < 0.00001). In both conserved and nonconserved domains, Dr/Dc ratios are significantly greater in the outcrosser (P < 0.00001) than in the selfer species pair. Nevertheless, this difference is more pronounced in the conserved than in the nonconserved domain stretches.

The preceding analysis focused on differences in dN/dS and Dr/Dc ratios among conserved and nonconserved regions “within” genes for the selfer and outcrosser species pairs. The following analysis focuses on analogous predictions “between” species pairs. We found that the dN/dS ratio was smaller in the selfer than in the outcrosser species pair and this difference was more pronounced in conserved gene regions (table 2 and fig. 4). Specifically, in the conserved domains of genes, the dN/dS ratio was approximately 116% (PS vs. O < 0.00001) higher in the outcrosser than in the selfer species pair. In contrast, in nonconserved regions of genes the difference in dN/dS ratios between selfer species and outcrosser species, while significant (PS vs. O < 0.00001), was less pronounced (73% greater in outcrosser species than in selfing species). Repeating this analysis by excluding genes that may have been subject to positive selection (dN/dS > 1) did not affect this observation (table 2; dN/dS in conserved stretches: PS vs. O < 0.00001; dN/dS in nonconserved stretches: PS vs. O < 0.001), and neither did the analysis of one-to-one orthologous gene sets (supplementary table S2, Supplementary Material online).

Finally, we performed an analogous analysis for the ratio of radical to conservative amino acid changes (Dr/Dc). Again we found that the difference in Dr/Dc values between selfer and outcrosser species pairs was more pronounced in conserved than in nonconserved regions of genes (table 3 and fig. 5). Specifically, Dr/Dc ratios in conserved regions of genes were approximately 26% (PS vs. O < 0.00001) greater in the outcrosser species pair and this difference was statistically highly significant. In contrast, in nonconserved regions of genes Dr/Dc values were only approximately 18% (PS vs. O < 0.0001) greater between outcrossers and selfers. This observation was not affected by excluding genes that may have been subject to positive selection (table 3; Dr/Dc in conserved stretches: PS vs. O < 0.00001; Dr/Dc in nonconserved stretches: PS vs. O < 0.00001), nor by considering only one-to-one orthologs (supplementary table S3, Supplementary Material online). All these observations imply that deleterious mutations are more effectively purged in the selfers than in the outcrossers.

Discussion

Selfing has two opposing effects on the efficacy of purifying selection. On the one hand, selfing reduces effective population size, because it increases genome-wide homozygosity (Nordborg 2000) and lowers effective recombination rates (Wright et al. 2008, 2013). These effects reduce the efficacy of purifying selection. On the other hand, selfing exposes recessive deleterious mutations in a homozygous state and can thus also enhance purging (Glémin 2007). In selfers with a diploid-dominant life cycle the first effect dominates, thus leading to a net reduction in purifying selection (Slotte et al. 2013; Wright et al. 2013). Here we show that in plants with a haploid-dominant life cycle undergoing intragametophytic selfing, an extreme form of selfing, homozygous exposure of deleterious mutations can lead to enhanced purging on a genome-wide scale.

This assertion is supported by three lines of evidence. First, we found that selfers accumulate fewer nonsynonymous mutations (dN/dS) and radical amino acid replacements (Dr/Dc) than outcrossers. Second, within genes, significantly fewer deleterious mutations (dN/dS and Dr/Dc) accumulate in conserved than in nonconserved regions, and this effect is greater in selfers than in outcrossers. Finally, by contrasting the accumulation of deleterious mutations (estimated again through dN/dS and Dr/Dc) among the two species pairs, we found that the difference between selfers and outcrossers is more pronounced in conserved than in nonconserved gene regions.

Previous studies provided indirect evidence that intragametophytic selfing can facilitate purging, but did not directly assess the effect of intragametophytic selfing on the accumulation of deleterious mutations at the level of genes and genomes. For instance, they found a lower incidence of lethality in homozygous embryos of ferns and a lack of inbreeding depression in homozygous sporophytes of mosses capable of intragametophytic selfing (Klekowski 1973, 1982; Lloyd 1974; Hedrick 1987a, 1987b; Taylor et al. 2007). By studying constrained sequence evolution on a genome-wide scale, our study goes beyond previous phenotypic observations and provides direct evidence that purifying selection under intragametophytic selfing can balance and even outweigh the effect of effective population size reduction caused by selfing (Hedrick 1987a, 1987b; Holsinger 1987). Furthermore, our observations also suggest that the correlation between intragametophytic selfing and a lack of both homozygous lethality and inbreeding depression may be explained by a decreased accumulation of deleterious mutations in the genome. Any direct connection between deleterious mutations and observed phenotypic effects, however, has to be experimentally demonstrated in future work.

It is generally assumed that the potential masking of deleterious mutations is restricted to diploid-specific genes of a plant’s lifecycle. The reason is that deleterious mutations in haploid-specific genes (as well as in genes needed in both stages) are always exposed to selection. Previous work showed that only approximately 2% of genes expressed in F. hygrometrica and P. patens show diploid-specific expression (Szövényi et al. 2011; O’Donoghue et al. 2013). Such a minute fraction of the genome subject to effective purging by intragametophytic selfing might go undetected. However, the significant effect of intragametophytic selfing on purging we observed shows otherwise. Together with our previous finding that diploid- and haploid-specific genes experience similar efficacy of purging (Szövényi et al. 2013), it contradicts the notion that only diploid-specific genes are affected by intragametophytic selfing.

To explain this apparent contradiction, we note that a gene’s exposure to selection in the haploid phase does not necessarily lead to more effective selection (Gerstein et al. 2011; Gerstein and Otto 2011; Billiard et al. 2012; Arunkumar et al. 2013; Szövényi et al. 2013; Zörgö et al. 2013). Therefore, contrary to previous assumptions, genes expressed in both phases can experience the effect of intragametophytic selfing on purging. One possible mechanism is that homozygous exposure of deleterious mutations via intragametophytic selfing may extend the time span over which deleterious mutations are directly exposed to selection (once as haploids and once as diploid homozygotes). Deleterious mutations, which may escape effective purifying selection in the haploid stage, may no longer escape natural selection when exposed in both stages. Moreover, the same mutation may have different deleterious effects in haploid and diploid stages, and their compound effect may be more severe, which can further contribute to its efficient purging. Because existing theory does not fully account for complications like this, more work is needed to fully understand the evolutionary dynamics of deleterious mutations in organisms with complex life cycles.

Because our findings are in stark contrast to the observation that selfing leads to reduced selection efficacy (Glémin and Galtier 2012), it is important to address the effect of confounding factors that may bias our analyses. The first of two major factors is that selfer and outcrosser species pairs may express different gene sets that could bias our analysis (Yang and Gaut 2011). This gene set bias may be partially caused by differential tissue sampling in the selfer and outcrosser species pairs: RNA was obtained only from gametophytic tissues for the outcrosser whereas both gametophytic and sporophytic tissues were used for the selfer species pair. For this reason we repeated our analyses with genes that were one-to-one orthologous across all four study species, which yielded the same conclusions (see supplementary material, Supplementary Material online). Therefore, gene set bias cannot explain the differences we see between the selfer and the outcrosser species pairs.

Second, it is well-known that estimates of synonymous divergence dS (a measure of divergence time) and dN/dS are negatively correlated. That is, dN/dS values are higher between recently diverged species, because intraspecific polymorphisms contribute prominently to dN/dS (Rocha et al. 2006; Kryazhimskiy and Plotkin 2008; Wolf et al. 2009; Mugal et al. 2014). Being aware of this issue, we selected species pairs (selfer and outcrosser) that exhibit similar divergence times. More specifically, both the outcrosser and the selfer species pairs are estimated to have diverged about 14–20 Ma (Shaw et al. 2010; Szövényi et al. 2013) which corresponds to an approximate divergence time of 23–33 Ne generations (assuming an Ne of 600,000 [McDaniel, van Baren et al. 2013] and a mutation rate of 1.8 × 10−8 [Rensing et al. 2007]). Previous work suggests that the contribution of polymorphisms to dN/dS estimates within this time range is small and the overall influence of divergence time on total dN/dS is thus weak (Mugal et al. 2014). Because this argument depends on the parameter estimates used, we also conducted a second, complementary analysis. Specifically, we divided gene alignments into five categories according to their dS values in both species pairs and compared dN/dS values for each category separately. We found that dN/dS ratios were significantly greater for the outcrosser species pair in all dS categories investigated, both for the full and for the one-to-one ortholog data set (supplementary table S4, Supplementary Material online). All this suggests that the difference in dN/dS ratios between selfers and outcrossers is not caused by the difference in synonymous divergence.

Conclusions, Limitations, and Open Questions

Although our study provides evidence that purging under intragametophytic selfing can be more effective than under intergametophytic selfing, multiple questions remain to be resolved. First, our current analysis is based on the comparison of one selfer and one outcrosser species pair. Therefore, with the current data set in hand it is difficult to assess the generality of our findings. Including multiple selfer and outcrosser lineages in future studies will be necessary to test the generality of our findings. Second, the extent to which diploid-specific genes or genes expressed in both ploidal phases are affected by purging has to be determined. This is critical to understand the evolutionary mechanism of purging. Third, our data suggest that intragametophytic selfing is advantageous, because it can efficiently prevent the accumulation of deleterious mutations in organism with extreme selfing. Theory predicts, however, that extreme linkage disequilibrium established by repeated events of intragametophytic selfing may hinder adaptive evolution (Felsenstein 1974; Birky and Walsh 1988; Charlesworth 2012). Nevertheless, McDaniel, Atwood et al. (2013) showed that hermaphroditic lineages of haploid-dominant plants have an elevated diversification rate compared with dioecious lineages. Therefore, the better purging in haploid selfers may also have long-term macroevolutionary consequences that need to be explored further. Thus, the extent to which intragametophytic selfing constrains the adaptive potential of species remains to be determined.

Supplementary Material

Supplementary tables S1–S4 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).