Adaptive reduction of male gamete number in the selfing plant arabidopsis thaliana

ABSTRACT The number of male gametes is critical for reproductive success and varies between and within species. The evolutionary reduction of the number of pollen grains encompassing the

male gametes is widespread in selfing plants. Here, we employ genome-wide association study (GWAS) to identify underlying loci and to assess the molecular signatures of selection on pollen

number-associated loci in the predominantly selfing plant _Arabidopsis thaliana_. Regions of strong association with pollen number are enriched for signatures of selection, indicating

polygenic selection. We isolate the gene _REDUCED POLLEN NUMBER1_ (_RDP1_) at the locus with the strongest association. We validate its effect using a quantitative complementation test with

CRISPR/Cas9-generated null mutants in nonstandard wild accessions. In contrast to pleiotropic null mutants, only pollen numbers are significantly affected by natural allelic variants. These

data support theoretical predictions that reduced investment in male gametes is advantageous in predominantly selfing species. SIMILAR CONTENT BEING VIEWED BY OTHERS THE IDENTIFICATION OF

THE _ROSA S_-LOCUS AND IMPLICATIONS ON THE EVOLUTION OF THE ROSACEAE GAMETOPHYTIC SELF-INCOMPATIBILITY SYSTEMS Article Open access 12 February 2021 THE EVOLUTION AND MAINTENANCE OF TRIOECY

WITH CYTOPLASMIC MALE STERILITY Article Open access 14 October 2024 OUTCROSSING RATES IN AN EXPERIMENTALLY ADMIXED POPULATION OF SELF-COMPATIBLE AND SELF-INCOMPATIBLE _ARABIDOPSIS LYRATA_

Article Open access 16 December 2021 INTRODUCTION Male gamete numbers, reflected by pollen grain (each containing two sperm cells) numbers in seed plants and sperm numbers in animals, have

been studied extensively from agricultural, medical, and evolutionary viewpoints1,2,3,4,5,6,7,8,9. Evolutionary theory predicts that the breeding system could act as a major selective force

on male gamete numbers. In highly promiscuous outcrossing species, a large number of male gametes should be produced because of male–male gamete competition, so reduced male gamete numbers

are considered to be deleterious1. In contrast, they may be advantageous at lower outcrossing rates because of the high cost of their production, decreasing fitness. In an agricultural

context, low pollen number may have been selected during domestication10, but may serve as a barrier for hybrid breeding of wheat and other species11,12. In flowering plants, the transition

from an outcrossing to a selfing breeding system through loss of self-incompatibility is one of the most prevalent evolutionary trends6,13. Selfing populations or species generally show

lower pollen grain numbers per flower (hereafter, pollen number) as well as reduced flower size. There has been a sustained debate on whether the reduced pollen number is a result of

adaptive evolution or the accumulation of deleterious mutations owing to reduced selection6,14,15, but little was known about the genetic basis of pollen number variation to assess molecular

signatures of selection on it. To unravel the genetic basis of quantitative natural variation in pollen number, we here focus on the predominantly selfing plant _Arabidopsis thaliana_6,16.

Studies have shown that the evolution of predominant selfing in _A. thaliana_ occurred much more recently than its evolutionary divergence from outcrossing relatives6,16,17. Thus, in

addition to fixed, genetically based differences from these outcrossing relatives, we expect that variation with regard to pollen number may still be segregating among current accessions. By

harnessing the genetic and genomic resources available in _A. thaliana_, we conduct a genome-wide association study on pollen number variation. We find that natural variants of the _RDP1_

gene confer variation in pollen number without detectable pleiotropy. Signatures of selection at the top genome-wide association study (GWAS) peaks, including the _RDP1_ locus, support the

theoretical prediction that reduced investment in male gametes should provide an advantage in selfing species. RESULTS GENOME-WIDE ASSOCIATION STUDY AND SIGNATURES OF SELECTION To examine

variability in pollen number on a species-wide scale, we determined pollen number per flower in 144 natural _A. thaliana_ accessions (Fig. 1a–d; Supplementary Table 1 and Supplementary Data 

1) and found approximately fourfold variation (average ~4000) (Fig. 1e). Histological sections of stamens from representative accessions confirmed pollen number variation among accessions

(Fig. 1c, d). We also measured the number of ovules per flower (Supplementary Table 2). We did not find significant correlations between numbers of pollen grains and ovules (_P_ = 0.5164),

although negative correlations have often been reported in between-species comparisons, as expected on theoretical grounds owing to trade-offs in resource allocation to male versus female

function14,18. Furthermore, we found that pollen number per flower was not significantly correlated with any of the 107 published phenotypes of flowering, defense-related, ionomic, and

developmental traits (Supplementary Table 3; Supplementary Note 1)19, nor with climate variables, geographic location, or _S_-haplogroups across the 144 accessions (Supplementary Tables 4

and 5, Supplementary Fig. 1)20,21. These data suggest that variation in pollen number is largely independent of other traits. To evaluate genome-wide signatures of natural selection on loci

associated with gamete numbers, we first performed GWAS for pollen and ovule numbers using a genome-wide single-nucleotide polymorphism (SNP) data set for these lines, which was obtained by

imputation based on genome-wide resequencing data and 250 k SNP data (Fig. 1f, h; Supplementary Fig. 2)22. In total, 68 peaks of association were identified (10-kb windows having SNPs with

_P_ < 10–4), although only one pollen number-associated peak remained significant after Bonferroni correction. Focusing on the identified GWAS peaks, we performed an enrichment analysis

to ask whether pollen and ovule number-associated peaks are enriched in long-haplotype regions, which could be owing to partial or ongoing sweeps of segregating polymorphisms23,24,25. To

identify long-haplotype regions, we first calculated the extended haplotype homozygosity (EHH), which measures decay of haplotypes that carry a specified core allele as a function of

distance23. We then obtained the integrated haplotype score (iHS) statistic for each SNP, which compares EHH between two alleles of the SNP by controlling for the allele frequency of each

SNP23. We found that 10-kb windows including pollen number-associated loci were significantly enriched in extreme iHS tails (_P_ < 0.05, permutation test; Fig. 1i, j; Supplementary Fig. 

3). These loci showed generally high iHS scores, and two of the top five GWAS peaks were outliers of the genome-wide iHS distribution (Supplementary Table 6). The enrichment was robust to

changes in sample composition, allele frequency cutoffs, and the use of windows (Supplementary Figs. 4–6; see Supplementary Note 2 for details). Ovule number also showed enrichment, albeit

less than pollen number (Fig. 1j). In principle, the iHS enrichment could be confounded by recombination rate and the accuracy of imputation. To deal with such potential confounding factors,

we compared these results with the results of an iHS enrichment analysis for GWAS peaks (_P_ < 0.0001) for 107 other phenotypes, as these confounding factors should also influence the

enrichment for other traits. We found that the iHS enrichment for pollen number GWAS peaks (_P_ = 0.002 for the top 1% iHS tail; Supplementary Fig. 3) was among the highest, compared with

that for many known adaptive traits included in the 107 phenotypes, such as leaf number at flowering time and resistance to _Pseudomonas_ pathogens19,25 (Supplementary Table 7). In addition,

iHS enrichment of the ovule number GWAS peaks was also significant (_P_ = 0.030 for the top 1% iHS tail; Supplementary Fig. 3). These enrichments support polygenic selection on a

considerable number of loci associated with male and female gamete numbers throughout the genome. ISOLATION OF THE _REDUCED POLLEN NUMBER1_ GENE To further understand the molecular basis of

pollen number variation and to examine the nature of the putative targets of selection, we tried to identify the genes underlying pollen number variation; however, the top five peaks of

association did not contain any genes with known functions in early stamen or pollen development. To obtain experimental evidence concerning genes underlying pollen number variation, we

conducted functional analyses of the genes under the highest pollen number GWAS peak, which explains ~20% of the total phenotypic variance between accessions and satisfies the criterion for

genome-wide significance (–log10 _P_ = 7.60). This region is of particular interest because it also satisfies the criteria for genome-wide significance of the iHS statistic (_P_ = 0.0149;

Fig. 1i, Supplementary Fig. 7), suggesting a selective sweep. To test whether the signature of selection in this region might be owing to traits other than pollen number, we examined whether

there is an association signal for any of the 107 published phenotypes, ovule number, or variants showing climatic correlations19,20. In the 10-kb window including the SNP of the highest

GWAS score for pollen number, we found no genotype–phenotype associations below _P_ < 10–5 or climate–SNP correlations below an empirical _P_ < 0.01, i.e., there is no significant

evidence for selection on traits other than pollen number in this region. We also found that accessions with the long-haplotype variants produced lower pollen numbers than those with

alternative haplotype variants (_P_ = 2.152 × 10–6, _t_ test; population structure-corrected GWAS _P_ = 2.95 × 10–6; Fig. 1k), as expected if this haplotype was under selection for reduced

pollen number. Of the three genes in this chromosomal region with the highest GWAS scores (AT1G25250, AT1G25260, and AT1G25270; Fig. 1g), the expression level of AT1G25260, a gene of unknown

function, was much higher in flower buds than that of the other two genes (Supplementary Fig. 8). Therefore, we obtained two T-DNA insertion mutants of AT1G25260 in the standard Col-0

accession from the Nottingham _Arabidopsis_ Stock Centre26. These mutants showed a 32% reduction in pollen number (Fig. 2a; Supplementary Table 8). We hereafter refer to AT1G25260 as

_REDUCED POLLEN NUMBER1_ (_RDP1_). Because both _rdp1-1_ (insertion in the 5′ UTR) and _rdp1-2_ (insertion at the end of the coding sequence) (Fig. 2b) homozygotes showed low levels of

expression, they are likely to be hypomorphic mutants (Supplementary Fig. 8). We generated two amorphic (null) frameshift mutants of _RDP1_ (_rdp1-3_ and _rdp1-4_) using the CRISPR/Cas9

system27,28 in the Col-0 background. These mutants indeed showed an even greater reduction in pollen number, but still produced about half the number of pollen grains of the corresponding

wild type (53% for _rdp1-3_; Fig. 2a), suggesting a quantitative nature of the effect of _RDP1_. Pollen size was slightly increased, in agreement with the well-known negative relationship

between pollen number and size, even within the same genotype (Supplementary Fig. 9; Supplementary Table 9)14. The mutant phenotype was complemented by transforming a 4.3-kb genomic fragment

of the Col-0 accession encompassing _RDP1_ (Fig. 2a, Supplementary Fig. 9). In contrast to _RDP1_ mutants, CRISPR/Cas9 induced null mutants in AT1G25250 or AT1G25270 did not result in any

significant change in pollen number (Supplementary Fig. 10). The phenotype of four independent mutants of _RDP1_ together with successful complementation using the wild-type allele thus

demonstrated that _RDP1_ is involved in the control of pollen number. Based on phylogenetic analysis, RDP1 is a putative homolog of the yeast mRNA turnover 4 protein (Mrt4p) (Supplementary

Figs. 11 and 12). The _MRT4_ gene is nonessential in yeast but null mutants show a phenotype of slightly slower growth. The Mrt4p shares similarity with the ribosome P0 protein and is

necessary for the assembly of the P0 protein into the ribosome. The human ribosome P0 gene is reported to have an extra-ribosomal function in cancer by modulating cell proliferation29,30.

During anther development, sporogenous cells first divide and differentiate into microsporocytes31. Following meiosis of the microsporocyte, four microspores are formed, each of which

undergoes two mitotic divisions to form a mature pollen grain containing the male gametes. The null _rdp1_-_3_ mutant produced fewer microsporocytes than the wild type (Fig. 2c–e),

indicating a reduction in cell numbers before meiosis. Consistent with this, in situ mRNA hybridization experiments detected strong expression of _RDP1_ in sporogenous cells and the

microsporocytes derived from them, but not in microspores (Fig. 2f, g; Supplementary Fig. 13). _RDP1_ was also expressed in other proliferating cells, including those in inflorescences,

floral meristematic regions, and ovules (Supplementary Figs. 8 and 13), supporting a more widespread role in proliferating cell types with putatively high demands for ribosome biogenesis32.

Furthermore, fusing the _RDP1_ promoter to the _uidA_ reporter gene encoding β-glucuronidase (GUS) to assess its activity confirmed the _RDP1_ expression pattern in stamens (Supplementary

Fig. 14) and demonstrated marked expression in root tips and young leaf primordia during the vegetative phase; these data are supported by quantitative reverse transcription PCR experiments

(Supplementary Fig. 8b). Consistent with _RDP1_ expression in proliferating tissues, the _rdp1-3_ null mutant showed pleiotropic phenotypes, including slower vegetative growth and reduced

ovule numbers per flower (Supplementary Fig. 15). Because these pleiotropic phenotypes would be deleterious in natural environments, these data indicate that natural alleles of _RDP1_ are

not null variants (see below). In summary, these data suggest that _RDP1_ is required in proliferating _A. thaliana_ cells, yet the natural variants we have identified predominantly affect

the proliferation of sporogenous cells in the anthers, consistent with the function of its yeast homolog in cell proliferation. NATURAL VARIANTS OF _RDP1_ CONFER POLLEN NUMBER VARIATION It

has been difficult to experimentally determine whether a particular gene has natural alleles with subtle phenotypic effects on quantitative traits33. When allelic effects are subtle,

transgenic analysis of natural alleles is not sufficiently powerful because the phenotypes of transgenic _A. thaliana_ plants tend to be highly variable as a result of the variation between

lines, e.g., owing to different transgene insertion sites. In contrast, a quantitative complementation test can identify responsible genes by testing the effect of natural alleles in a

heterozygous state with a null allele if the effects of other loci are small, although this may be confounded by polygenic effects in the genetic background33,34. To conduct such

quantitative complementation, we took advantage of the CRISPR/Cas9 technique to generate frameshift null alleles in nonstandard natural accessions, in which no prior mutant was available. We

used Bor-4 and Uod-1, which have high and low pollen number phenotypes, respectively (Fig. 2b; Fig. 3). There were a number of sequence differences between the two accessions in the region

encompassing the _RDP1_ gene from 777 bp upstream of the start codon to 643 bp downstream of the stop codon. We found one non-synonymous and six synonymous substitutions in the coding

region, and 62 substitutions and six indel mutations in the non-coding region (Supplementary Fig. 16). Yet, both accessions did not reveal obvious loss-of-function mutations (Supplementary

Fig. 16), and _rdp1_ CRISPR null mutants of each accession showed reduced pollen number compared with the corresponding wild type (_P_ < 2.2 × 10–16 for Bor-4, _P_ = 9.84 × 10–7 for

Uod-1; Fig. 3a, b; Fig. 2h). These results show that both of the naturally occurring variants of the _RDP1_ gene are not null mutants but rather encode a functional protein. Disruption of

_RDP1_ had a stronger effect on pollen number in Bor-4 than in Uod-1 (analysis of variance (ANOVA) interaction effect _P_ = 1.07 × 10–5; Fig. 2h). This finding supports the notion that the

Bor-4 allele has a stronger promotive effect on pollen number than the Uod-1 allele, although other loci in the genetic backgrounds of Bor-4 and Uod-1 may contribute to this difference

through epistasis. To test the allelic effect of _RDP1_, we utilized a quantitative complementation test that controls for genetic background (Fig. 3). Among F1 plants obtained by crossing

heterozygotes for the frameshift mutation in each genetic background, we compared two genotypes: _RDP1_Bor_/rdp1_Uod vs. _rdp1_Bor_/RDP1_Uod. These F1 genotypes are identical except for the

differences at _RDP1_, where they both carry a frameshift allele but differ with respect to the functional allele; because of the crossing design, any independently segregating off-target

effects, resulting from CRISPR/Cas9 mutagenesis would be equally distributed between the two genotype cohorts of interest. We found that pollen number in plants with _RDP1_Uod was

significantly lower than in plants with _RDP1_Bor (nested ANOVA, 468 flowers from 26 individuals of _RDP1_Bor/_rdp1_Uod and 368 flowers from 20 individuals of _rdp1_Bor_/RDP1_Uod, _P_ = 4.85

 × 10–8; Fig. 3c). The significant difference cannot be attributed to stochastic individual differences, because the significant difference between plants bearing functional _RDP1_Bor and

_RDP1_Uod alleles was also observed in an individual-based test using averaged data of each individual separately (_P_ = 0.0331). Thus, in an otherwise identical genetic background, the

respective functional _RDP1_ haplotypes cause differences in pollen number. We also measured rosette leaf size, flowering date, ovule number, dry weight, and seed weight of plants in the two

cohorts, but none of these traits showed significant differences between _RDP1_Bor/_rdp1_Uod and _rdp1_Bor/_RDP1_Uod (Supplementary Fig. 17). Thus, in contrast to the experimentally

generated null mutants that showed pleiotropic growth defects, these results indicate that natural allelic differences at _RDP1_ affect pollen number in the absence of any detectable

deleterious pleiotropy. This finding is also supported by no genotype–phenotype associations for other traits, as described above. DISCUSSION We here isolated the _RDP1_ gene underlying

natural variation in male gamete numbers. Our study provides evidence for polygenic selection on pollen number-associated loci, including _RDP1_. Even though _RDP1_ encodes a

ribosome-biogenesis factor that would be required globally for proliferative growth, the naturally selected alleles predominantly confer reduced pollen number. This is analogous to a

hypomorphic allele of the human _G6PD_ gene, which encodes an enzyme in the pentose phosphate pathway and features a long haplotype because of selection for malaria resistance24. The mean

number of pollen grains of an outcrossing population of _A__rabidopsis_ _lyrata_ is ~18,000 (ref. 35) and thus several times higher than our counts in _A. thaliana_ (~2000–8000; Fig. 1e).

Although the evolutionary split between the lineages leading to _A. thaliana_ and _A. lyrata_ is estimated to have occurred ~5 million years ago or before36,37,38, several studies suggested

that predominant selfing in _A. thaliana_ evolved much more recently (0–0.413 million years ago based on the timing of the loss of a self-incompatibility gene, ~0.5 million years ago based

on the abundance of transposable elements, and ~0.3–1 million years ago based on the pattern of genome-wide linkage disequilibrium)6,39,40. Thus, in addition to presumably fixed differences

to distantly related outcrossing congeners, it is quite conceivable that some underlying loci, not limited to but including _RDP1_, are still segregating within _A. thaliana_. These might

reflect an ongoing selection process for the further reduction in pollen number, which has been considered a hallmark of the so-called selfing syndrome6,8,14. Although we note the

possibility that partial sweeps of _RDP1_ and other segregating loci are not directly related to the transition to predominant selfing, our analysis did not find evidence of other selective

forces including local adaptation, climate association, or pleiotropic selection on other traits. Therefore, our study supports the theoretical predictions that reduced investment in male

gametes is advantageous in predominantly selfing species. Our work also illustrates that a combination of GWAS and functional analysis using a quantitative complementation test based on the

CRISPR/Cas9-based alleles provides a powerful approach to dissect allelic differences underlying quantitative natural variation. METHODS POLLEN AND OVULE COUNTING FOR GENOME-WIDE ASSOCIATION

STUDIES To perform GWAS, numbers of pollen grains and ovules per flower were counted for 144 and 151 world-wide natural accessions, respectively (Supplementary Tables 1, 2 and Supplementary

Data 1). Plants were grown at 21 °C under a 16 h light/8 h dark cycle without vernalization. We grew four plants per accession. Three flower buds per plant were harvested from the main

inflorescence, and each flower bud was collected into a 1.5 mL tube and dried at 65 °C overnight. We sampled individual flower buds of young main inflorescences but avoided the first and

second flowers of the inflorescence because these flowers tend to show developmentally immature morphologies. We collected flower buds with mature pollen but before the anthers were opened

(flower stage 13), and added 30 μL of 5% Tween 20 (Sigma-Aldrich, St. Louis, MO, USA) to each tube. The tubes were sonicated using a Bioruptor (Diagenode, Seraing, Belgium) in high power

mode with 10 cycles of sonication-ON for 30 s and sonication-OFF for 30 s so that the pollen grains were released from the anther sacs. After a short centrifugation and vortexing, 10 μL of

the solution was mounted on a Neubauer slide. We took three images per sample using a light microscope. The number of pollen grains per image was counted using the particle counting

implemented in ImageJ (http://imagej.nih.gov/ij/) and in Fiji (http://fiji.sc/Fiji). We then estimated the total pollen number per flower based on the image size and the total volume. Ovule

numbers were counted by dissecting young siliques (5.3 siliques per accession on average) under a dissecting microscope. Because of limited chamber space, we split the plants into two

batches. The two batches were treated under the same conditions in the same chambers, but at different times. We controlled for this potential batch effect for pollen number by setting equal

medians and standard deviations for the two batches and used as the GWAS input of pollen number phenotype (Source Data file, Supplementary Table 1). Sometimes, there were no or very few

pollen grains per image. This was mainly in situations where anthers did not open. To eliminate these artefacts, we discarded flowers with pollen counts of <10 per image. We confirmed

that such extremely low pollen numbers did not occur in specific accessions, indicating that this is not heritable. PLANT MATERIALS AND GROWTH CONDITIONS FOR FUNCTIONAL ANALYSES For

functional analyses, we mainly used _Arabidopsis thaliana_ wild-type and mutant plants of the Col-0, Bor-4, and Uod-1 accessions. The T-DNA lines SALK_064854/N666274 (_rdp1-1_) from the Salk

collection41 and GK-879G09/N484369 (_rdp1-2_) from GABI-Kat collection42 were obtained from the European _Arabidopsis_ Stock Centre26. The T-DNA insertion in each line was confirmed using

PCR with the primers listed in Supplementary Table 10 as described at http://signal.salk.edu/tdnaprimers.2.html. DNA was extracted from young leaves using the cetrimonium bromide (CTAB)

method. We generated single-nucleotide insertion/deletion lines in the Col-0, Uod-1, and Bor-4 accessions using the FAST-CRISPR-Cas9 construct (see the section CRISPR mutant) and designated

them _rdp1-3_ and _rdp1-4_ (in Col-0), _rdp1-5_ (in Uod-1), and _rdp1-6_ (in Bor-4). _Arabidopsis_ seeds were sown on soil mixed with the insecticide ActaraG (Syngenta Agro, Switzerland) and

stratified for 3–4 days at 4 °C in the dark. The plants were grown under 16 h of light at 22 °C, and 8 h of dark at 20 °C, with weekly treatments of insecticide (Kendo Gold, Syngenta Agro)

unless noted otherwise (for GWAS, see above). STATISTICAL ANALYSIS Unless stated otherwise, statistical and population genetic analyses were performed using the statistical software R43. In

boxplots, bars indicate the median, boxes indicate the interquartile range, and whiskers extend to the most extreme data point that is no >1.5 times the interquartile range from the box,

with outliers shown by dots. HISTOLOGICAL ANALYSIS OF ANTHERS For histological analysis, inflorescences were fixed with formaldehyde:acetic acid:70% ethanol = 1:1:18 (FAA), dehydrated and

embedded in Technovit 7100 according to the manufacturer’s instructions (Heraeus Kulzer GmbH, Wehrheim, Germany). Five-micrometer sections were cut with a microtome (RM2145, Leica, Germany)

and stained with toluidine blue before observation under a Leica microscope (DM5000, Leica) equipped with a black-and-white camera (DFC345, Leica). CORRELATIONS WITH PUBLISHED GWAS RESULTS

AND CLIMATIC DATA To examine whether pollen and ovule numbers were correlated with any of the other 107 published phenotypes19, or with climate and geographic variables20, Pearson’s

correlation coefficients were calculated (Supplementary Tables 3 and 5). We also surveyed whether there were any SNPs significantly associated with climate variables in the 10-kb window

including the SNP of the highest GWAS score for pollen number (Chr1:8,850,000–8,860,000). The significance was based on the genome-wide empirical _P_ values20,25. Focusing on the same 10-kb

window, we also surveyed whether there were significant SNPs (_P_ < 10–5; minor allele frequency (MAF)>0.1) for the 107 published phenotypes in the GWAS. The correlation of pollen

number and the _S_-haplogroups21 was also tested. SNP IMPUTATION To perform GWAS, we generated a set of dense SNP markers that overlapped with the phenotyped accessions using imputation22.

First, we constructed a reference set of 186 haplotypes from resequencing data44,45. MaCH version 1.0.16.c46 was then used to impute non-genotyped SNPs for 1311 accessions in the 250 k SNP

array data set25. The command line used for each overlapping bin was: mach1 --dosage --greedy -r -d [sample_bin].dat -p [sample_bin].ped -s [ref_bin].snps -h [ref_bin].haplos --prefix

[output_prefix] The output was then merged and converted into the homozygous SNP data set. GENOME-WIDE ASSOCIATION STUDY GWAS was performed to identify loci associated with pollen number

variation in 144 natural accessions. We also performed GWAS for ovule number and for 107 published phenotypes19 using the same SNP set as for pollen number. The median values of pollen

number were used to represent each accession. We used the log-transformed pollen number values for GWAS, as they did not deviate significantly from a normal distribution (Shapiro–Wilk

normality test: _P_ = 0.438). To deal with confounding effects as a result of population structure, we employed the mixed model implemented in the software Mixmogam, in which a genome-wide

kinship matrix was incorporated as a random effect (population structure-corrected GWAS)47. By using Mixmogam, we also generated a quantile–quantile (Q–Q) plot, which shows the relationship

between the observed and expected negative logarithm of _P_ values (Fig. 1h). For the Manhattan plot (Fig. 1f, g), we removed SNPs with minor allele frequencies <0.15, leaving 1,115,178

SNPs overall. QUANTITATIVE REVERSE TRANSCRIPTION PCR For gene expression analysis, total RNA was isolated using the ZR Plant RNA MiniPrep (Zymo Research, Irvine, CA) according to the

manufacturer’s instructions. Total RNA was treated with a DNA-free DNA removal kit (Thermo Fisher Scientific, Waltham, MA), and cDNA was synthesized using a High-Capacity RNA-to-cDNA kit

(Thermo Fisher Scientific, Waltham, MA). For the quantitative PCR step, cDNA was used with SYBR Select Master Mix (Thermo Fisher Scientific, Waltham, MA) and gene-specific primers

(Supplementary Table 10). Data were collected using the StepOnePlus Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA) in accordance with the instruction manual. Expression levels

were normalized using _EF1-α_ (AT5G60390), which was used as an internal control. CRISPR/CAS9-BASED MUTANT ALLELES To generate _RDP1_ frameshift mutants, the CRISPR/Cas9 system was used28. A

20-nucleotide target sequence (5′-GCAGAAGATGAACTCCGTTC-3′) was selected using the CRISPR-P tool48. The target sequence was subcloned into a pTTK194 vector (pUC19 U6.26pro::sgRNA). Then, the

_Hin_dIII sgRNA fragment was subcloned into a pTTK182 (pFAST-R02 35Spro::Cas9) vector or pTTK227 (pFAST-R01 RPS5Apro::Cas9). For plant transformation, the binary vector was first introduced

into _Agrobacterium tumefaciens_ (GV3101) and then into _A. thaliana_ by the floral-dip method. TRANSGENIC EXPERIMENT The _RDP1_ genomic sequence was amplified by PCR using Phusion High

Fidelity PCR polymerase (New England Biolabs, Beverly, MA). We used the primers 3550_At1g25260F1 (5′-TTTCTCCCCACATTTCTC-3′) and 3551_At1g25260R1 (5′-GTTTAAAATGAGAGAACCCG-3′) to amplify the

full length of _RDP1_ and the surrounding genomic sequence that may contain promoter and terminator regions (ca. 4.3 kbp, positions 8,857,725 to 8,853,420 on chromosome 1). The _RDP1_

promoter sequence was amplified by PCR using the following primers: 3550_At1g25260F1 (5′-TTTCTCCCCACATTTCTC-3′) and 4155_At1g25260R (5′-AGCTGAGGTTTCAAATGTTGTATG-3′). These sequences were

cloned into a pCR8 vector using the pCR8/GW/TOPO TA Cloning kit (Thermo Fisher Scientific, Waltham, MA). pCR8-RDP1 complete sequences were cloned into a pFAST-G01 vector49, and the pCR8-RDP1

promoter vector was cloned into a pGWB3 vector50 using Gateway LR clonase II (Thermo Fisher Scientific, Waltham, MA). We designated these constructs pFAST-RDP1 and pGWB3-RDP1pro,

respectively. Plants were grown at 22 °C under a 16 h light/8 h dark cycle until transformation via _A. tumefaciens_ strain GV3101 using the floral-dip method51. pFAST-RDP1 was transformed

into _rdp1-1_ homozygous mutants. Transformed seeds were selected using a fluorescence stereomicroscope with a GFP filter (Olympus SZX12, Japan). pGWB3-RDP1pro was transformed to Col-0.

Seedlings carrying pGWB3-RDP1pro were selected on 1/2 Murashige and Skoog medium supplemented with 50 μg mL−1 kanamycin and 50 μg mL−1 hygromycin. MEASURING POLLEN NUMBER AND SIZE WITH A

CELL COUNTER To expedite pollen number counting, we established a rapid method using a cell counter (CASY TT, OMNI Life Science GmbH, Germany)52. We found that the pollen numbers of flowers

on side inflorescences and side branches of the main inflorescence were similar, but those of flowers on the main inflorescence were higher (Supplementary Fig. 18). To obtain a large number

of replicates, we sampled flowers from the former two positions. We sampled during the first 3 weeks of flowering and excluded the first and second flowers on each branch; this yielded up to

40 flowers per individual. Collecting, suspending, and sonicating of flowers for GWAS were conducted as described above. All pollen solutions were suspended in 10 mL of CASYton (OMNI Life

Science GmbH), and pollen numbers were counted with a CASY TT cell counter53. Three 400 μL aliquots of each pollen solution were counted. We counted particles within a size range of 12.5–25 

μm (estimated diameter) as pollen. Samples with a clear peak at 7.5–12.5 μm were discarded as broken or unhealthy samples. COUNTING MICROSPOROCYTE NUMBER Flower buds were collected and fixed

in 3:1 (vol:vol) ethanol:acetic acid for 16 h. The fixed flower buds were rehydrated in 99.5%, 70%, and 50% ethanol for 30 min each. Flower buds were transferred to 1 n NaOH for 3 h and

then stained with aniline blue solution (0.1% aniline blue, 100 mm K3PO4) for 3 h. Anthers at the microsporocyte stage were dissected under a microscope with UV illumination (DM5000, Leica,

Germany). Z-stack images were obtained with a confocal microscope (SP5, Leica, Germany). IN SITU HYBRIDIZATION Flower buds were fixed with 3.7% formaldehyde, 5% acetic acid, 50% ethanol

(FAA) and dehydrated through an ethanol series. Fixed samples were embedded in paraplast using an embedding machine (ASP200, Leica, Germany). _RDP1_ cDNA was PCR-amplified using the primer

pair (At1g25260g2F; 5′-tgcctaatcaaagcgagtagacc-3′ and At1g25260gR; 5′-cagagcaagttcagcttgaaagtagc-3′) and cloned into pCR4-TOPO (Thermo Fisher Scientific) vector. Cloned cDNA was used as a

template for in vitro transcription using a MAXIscript T7 labeling kit (Thermo Fisher Scientific) for hybridization54. GUS ASSAYS Plant samples were incubated in 90% acetone for 20 min at

room temperature, washed with 50 mm phosphate buffer containing 0.1% Triton X-100, 2 mm potassium ferrocyanide, and 2 mm potassium ferricyanide, and incubated in the same buffer supplemented

with 1 mg mL−1 X-Gluc for 5 h at 37 °C55. PHYLOGENETIC ANALYSIS Multiple sequence alignment was performed using ClustalW implemented in the CLC Workbench (version 7.7). A phylogenetic tree

was generated by the neighbor-joining distance algorithm, using the aligned region (amino acid positions 44–153) and a bootstrap value of 1000. Yvh1 of _Saccharomyces cerevisiae_ was used as

an outgroup. Accession numbers of used sequences are listed in Supplementary Table 11. SELECTION SCAN For the selection scan, we used the imputed SNP data set that was also used for GWAS.

We used 298 accessions (Supplementary Data 1), covering all the accessions used for our GWAS of pollen and ovule numbers and the GWAS of 107 phenotypes reported by Atwell et al.19,25,56. We

used the iHS statistic for the selection scan; this statistic compares the EHH between two alleles by controlling for the allele frequency of each SNP23. The iHS statistic uses the contrast

of EHH values on each SNP; iHS values strongly deviate from zero when one allele has a long haplotype (high EHH) and the other has a short one (low EHH); the R library rehh57 was used to

calculate the iHS statistic. The _Arabidopsis lyrata_ reference genome58 was used to infer the ancestral state for each SNP. We first calculated the iHS for each SNP. Then, we split the

genome into 10-kb windows and used the maximum score from the iHS scan for each window as the test statistic, as LD is known to decay within ~10-kb on average in geographically world-wide

samples of _A. thaliana_59. Empirical _P_ values were calculated for all windows and for all SNPs, based on their ranks in genomic distributions. To deal with the geographically biased

sampling, which could be a possible confounding factor, we also performed the selection scan with 144 accessions that were used for GWAS (see Supplementary Note 2 for details; Supplementary

Fig. 5). We then asked whether the GWAS-associated windows were enriched in the extreme tails of the iHS statistic. Enrichment analysis was performed across the 10-kb windows. To examine

whether enrichment of the iHS was commonly observed in GWAS peaks, we also performed this analysis for the GWAS results of the other 107 publicly available phenotypes19 and compared them

with the GWAS of pollen and ovule numbers. Four thresholds of empirical _P_ value tails were considered: 10%, 5%, 2.5%, and 1%. GWAS SNPs were considered if _P_ values were smaller than 10–4

and the MAF was >0.1. We also performed the same enrichment analysis with MAF > 0.15 (Supplementary Fig. 4). To assess the effect of using a 10-kb window size, we also performed the

iHS enrichment analysis on a per-SNP basis, in addition to the 10-kb window size (Supplementary Fig. 6). The statistical significance of the fold-enrichment was addressed based on

permutation tests that preserve the linkage disequilibrium structure in the data. A set of windows was resampled for each permutation, preserving the relative positions of the windows, but

shifting them by a randomly chosen uniformly drawn number of windows for each permutation. A similar method of permutation has been used in several population genomic studies19,25,56.

Permutation was performed 1000 times. MEASURING PLANT PHENOTYPES For measuring rosette leaf size, we took images of plants that included a ruler at 3 weeks after germination. A minimum

circumscribed circle was drawn manually on the picture using Fiji; then, the area was measured and transformed depending on the scale. The flowering date was counted as days from sowing to

flowering. The dry weight of plants was measured using aerial parts, and seed weight was determined by collecting seeds from dried plants. _P_ values for quantitative complementation test

(Supplementary Fig. 17) are shown in Supplementary Table 12. REPORTING SUMMARY Further information on experimental design is available in the Nature Research Reporting Summary linked to this

paper. DATA AVAILABILITY Data supporting the findings of this work are available within the paper and its Supplementary Information files. A Reporting Summary for this Article is available

as a Supplementary Information file. The data sets generated and analyzed during the current study are available from the corresponding author upon request. _RDP1_ gene sequence data

generated in this article were registered in GenBank (National Center for Biotechnology Information) databases under the following accession numbers: LC164158 (Mz-0), LC164159 (Bor-4),

LC504218 (Uod-1), LC164160 (_A. lyrata_), and LC164161 (_Arabidopsis halleri_). _P0_ gene sequence data of _A. halleri_ generated in this Article were registered in GenBank databases under

the following accession numbers: LC164162 and LC164163. Raw and processed sequencing data for GWAS and population genetic analyses are publicly available at

https://doi.org/10.5061/dryad.jh9w0vt7z (ref. 60). The source data underlying Figs. 1e, 1j, 1k, 2a, 2e, 2h, and 3, as well as Supplementary Figs. 1, 3–6, 7a, 7d, 8, 9, 10c, 10d, 15, 17, and

18 are provided as a Source Data file. CODE AVAILABILITY Code used for selection enrichment analysis is publicly available at https://github.com/tsuchimatsu/arabidopsis_pollen_number.

REFERENCES * Harvey, P. H. & May, R. M. Out for the sperm count. _Nature_ 337, 508–509 (1989). Article  ADS  CAS  PubMed  Google Scholar  * Birkhead, T. R. & Møller, A. P. _Sperm

Competition and Sexual Selection_ (Academic Press, 1998). * Merzenich, H., Zeeb, H. & Blettner, M. Decreasing sperm quality: a global problem? _BMC Public Health_ 10, 24 (2010). Article

  PubMed  PubMed Central  Google Scholar  * Czeizel, A. E. & Rothman, K. J. Does relaxed reproductive selection explain the decline in male reproductive health? A new hypothesis.

_Epidemiology_ 13, 113–114 (2002). Article  PubMed  Google Scholar  * Kosova, G., Scott, N. M., Niederberger, C., Prins, G. S. & Ober, C. Genome-wide association study identifies

candidate genes for male fertility traits in humans. _Am. J. Hum. Genet._ 90, 950–961 (2012). Article  CAS  PubMed  PubMed Central  Google Scholar  * Shimizu, K. K. & Tsuchimatsu, T.

Evolution of selfing: recurrent patterns in molecular adaptation. _Annu. Rev. Ecol. Evol. Syst._ 46, 593–622 (2015). Article  Google Scholar  * Darwin, C. _The Effects of Cross and Self

Fertilisation in the Vegetable Kingdom_ (John Murray, 1876). * Sicard, A. & Lenhard, M. The selfing syndrome: a model for studying the genetic and evolutionary basis of morphological

adaptation in plants. _Ann. Bot._ 107, 1433–1443 (2011). Article  PubMed  PubMed Central  Google Scholar  * Barrett, S. C. H. The evolution of plant sexual diversity. _Nat. Rev. Genet._ 3,

274–284 (2002). Article  CAS  PubMed  Google Scholar  * Oka, H.-I. & Morishima, H. Variations in the breeding systems of a wild rice, _Oryza perennis_. _Evolution_ 21, 249–258 (1967).

Article  PubMed  Google Scholar  * Boeven, P. H. G. et al. Genetic architecture of male floral traits required for hybrid wheat breeding. _Theor. Appl. Genet._ 129, 2343–2357 (2016). Article

  PubMed  Google Scholar  * Langer, S. M., Longin, C. F. H. & Würschum, T. Phenotypic evaluation of floral and flowering traits with relevance for hybrid breeding in wheat (_Triticum

aestivum L_.). _Plant Breed._ 133, 433–441 (2014). Article  Google Scholar  * Goldberg, E. E. et al. Species selection maintains self-incompatibility. _Science_ 330, 493–495 (2010). Article

  ADS  CAS  PubMed  Google Scholar  * Cruden, R. W. Pollen grains: why so many? _Plant Syst. Evol._ 222, 143–165 (2000). Article  Google Scholar  * Willis, J. H. The contribution of

male-sterility mutations to inbreeding depression in _Mimulus guttatus_. _Heredity_ 83, 337–346 (1999). Article  PubMed  Google Scholar  * Tsuchimatsu, T. et al. Evolution of

self-compatibility in _Arabidopsis_ by a mutation in the male specificity gene. _Nature_ 464, 1342–1346 (2010). Article  ADS  CAS  PubMed  Google Scholar  * Tsuchimatsu, T., Kaiser, P., Yew,

C.-L., Bachelier, J. B. & Shimizu, K. K. Recent loss of self-incompatibility by degradation of the male component in allotetraploid _Arabidopsis kamchatica_. _PLoS Genet._ 8, e1002838

(2012). Article  CAS  PubMed  PubMed Central  Google Scholar  * Charnov, E. L. _The Theory of Sex Allocation_ (Princeton University Press, 1982). * Atwell, S. et al. Genome-wide association

study of 107 phenotypes in _Arabidopsis thaliana_ inbred lines. _Nature_ 465, 627–631 (2010). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Hancock, A. M. et al. Adaptation to

climate across the _Arabidopsis thaliana_ genome. _Science_ 334, 83–86 (2011). Article  ADS  CAS  PubMed  Google Scholar  * Shimizu, K. K., Shimizu-Inatsugi, R., Tsuchimatsu, T. &

Purugganan, M. D. Independent origins of self-compatibility in _Arabidopsis thaliana_. _Mol. Ecol._ 17, 704–714 (2008). Article  CAS  PubMed  Google Scholar  * The 1001 Genomes Consortium.

1,135 genomes reveal the global pattern of polymorphism in _Arabidopsis thaliana_. _Cell_ 166, 481–491 (2016). Article  CAS  Google Scholar  * Voight, B. F., Kudaravalli, S., Wen, X. &

Pritchard, J. K. A map of recent positive selection in the human genome. _PLoS Biol._ 4, e72 (2006). Article  PubMed  PubMed Central  Google Scholar  * Tishkoff, S. A. et al. Haplotype

diversity and linkage disequilibrium at human G6PD: recent origin of alleles that confer malarial resistance. _Science_ 293, 455–462 (2001). Article  CAS  PubMed  Google Scholar  * Horton,

M. W. et al. Genome-wide patterns of genetic variation in worldwide _Arabidopsis thaliana_ accessions from the RegMap panel. _Nat. Genet._ 44, 212–216 (2012). Article  CAS  PubMed  PubMed

Central  Google Scholar  * Scholl, R. L., May, S. T. & Ware, D. H. Seed and molecular resources for _Arabidopsis_. _Plant Physiol._ 124, 1477–1480 (2000). Article  CAS  PubMed  PubMed

Central  Google Scholar  * Jinek, M. et al. A programmable dual-RNA-guided DNA endonuclease in adaptive bacterial immunity. _Science_ 337, 816–821 (2012). Article  ADS  CAS  PubMed  PubMed

Central  Google Scholar  * Tsutsui, H. & Higashiyama, T. pKAMA-ITACHI vectors for highly efficient CRISPR/Cas9-mediated gene knockout in _Arabidopsis thaliana_. _Plant Cell Physiol._ 58,

46–56 (2017). Article  CAS  PubMed  Google Scholar  * Rodríguez-Mateos, M. et al. Role and dynamics of the ribosomal protein P0 and its related trans-acting factor Mrt4 during ribosome

assembly in _Saccharomyces cerevisiae_. _Nucleic Acids Res._ 37, 7519–7532 (2009). Article  PubMed  PubMed Central  CAS  Google Scholar  * Bhavsar, R. B., Makley, L. N. & Tsonis, P. A.

The other lives of ribosomal proteins. _Hum. Genomics_ 4, 327–344 (2010). Article  CAS  PubMed  PubMed Central  Google Scholar  * Sanders, P. M. et al. Anther developmental defects in

_Arabidopsis thaliana_ male-sterile mutants. _Sex. Plant Reprod._ 11, 297–322 (1999). Article  CAS  Google Scholar  * Breuninger, H. & Lenhard, M. Control of tissue and organ growth in

plants. _Curr. Top. Dev. Biol._ 91, 185–220 (2010). Article  CAS  PubMed  Google Scholar  * Turner, T. L. Fine-mapping natural alleles: quantitative complementation to the rescue. _Mol.

Ecol._ 23, 2377–2382 (2014). Article  CAS  PubMed  PubMed Central  Google Scholar  * Stern, D. L. A role of _Ultrabithorax_ in morphological differences between _Drosophila_ species.

_Nature_ 396, 463–466 (1998). Article  ADS  CAS  PubMed  PubMed Central  Google Scholar  * Willi, Y. Mutational meltdown in selfing _Arabidopsis lyrata_. _Evolution_ 67, 806–815 (2013).

Article  PubMed  Google Scholar  * Koch, M. A., Haubold, B. & Mitchell-Olds, T. Comparative evolutionary analysis of chalcone synthase and alcohol dehydrogenase loci in _Arabidopsis_,

_Arabis_, and related genera (Brassicaceae). _Mol. Biol. Evol._ 17, 1483–1498 (2000). Article  CAS  PubMed  Google Scholar  * Ossowski, S. et al. The rate and molecular spectrum of

spontaneous mutations in _Arabidopsis thaliana_. _Science_ 327, 92–94 (2010). Article  ADS  CAS  PubMed  Google Scholar  * Shimizu, K. K., Kudoh, H. & Kobayashi, M. J. Plant sexual

reproduction during climate change: gene function _in natura_ studied by ecological and evolutionary systems biology. _Ann. Bot._ 108, 777–787 (2011). Article  CAS  PubMed  PubMed Central 

Google Scholar  * Tang, C. et al. The evolution of selfing in _Arabidopsis thaliana_. _Science_ 317, 1070–1072 (2007). Article  ADS  CAS  PubMed  Google Scholar  * de la Chaux, N.,

Tsuchimatsu, T., Shimizu, K. K. & Wagner, A. The predominantly selfing plant _Arabidopsis thaliana_ experienced a recent reduction in transposable element abundance compared to its

outcrossing relative _Arabidopsis lyrata_. _Mob. DNA_ 3, 2 (2012). Article  PubMed  PubMed Central  CAS  Google Scholar  * Alonso, J. M. et al. Genome-wide insertional mutagenesis of

_Arabidopsis thaliana_. _Science_ 301, 653–657 (2003). Article  ADS  PubMed  Google Scholar  * Kleinboelting, N., Huep, G., Kloetgen, A., Viehoever, P. & Weisshaar, B. GABI-Kat

SimpleSearch: new features of the _Arabidopsis thaliana_ T-DNA mutant database. _Nucleic Acids Res_ 40, D1211–D1215 (2012). Article  CAS  PubMed  Google Scholar  * R Core Team. R Foundation

for Statistical Computing (The R Foundation, 2013). * Cao, J. et al. Whole-genome sequencing of multiple _Arabidopsis thaliana_ populations. _Nat. Genet._ 43, 956–965 (2011). Article  CAS 

PubMed  Google Scholar  * Long, Q. et al. Massive genomic variation and strong selection in _Arabidopsis thaliana_ lines from Sweden. _Nat. Genet._ 45, 884–890 (2013). Article  CAS  PubMed 

PubMed Central  Google Scholar  * Li, Y., Willer, C. J., Ding, J., Scheet, P. & Abecasis, G. R. MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes.

_Genet. Epidemiol._ 34, 816–834 (2010). Article  PubMed  PubMed Central  Google Scholar  * Segura, V. et al. An efficient multi-locus mixed-model approach for genome-wide association studies

in structured populations. _Nat. Genet._ 44, 825–830 (2012). Article  CAS  PubMed  PubMed Central  Google Scholar  * Lei, Y. et al. CRISPR-P: a web tool for synthetic single-guide RNA

design of CRISPR-system in plants. _Mol. Plant_ 7, 1494–1496 (2014). Article  CAS  PubMed  Google Scholar  * Shimada, T. L., Shimada, T. & Hara-Nishimura, I. A rapid and non-destructive

screenable marker, FAST, for identifying transformed seeds of _Arabidopsis thaliana_. _Plant J._ 61, 519–528 (2010). Article  CAS  PubMed  Google Scholar  * Nakagawa, T. et al. Improved

gateway binary vectors: high-performance vectors for creation of fusion constructs in transgenic analysis of plants. _Biosci. Biotechnol. Biochem._ 71, 2095–2100 (2007). Article  CAS  PubMed

  Google Scholar  * Clough, S. J. & Bent, A. F. Floral dip: a simplified method for _Agrobacterium_-mediated transformation of _Arabidopsis thaliana_. _Plant J._ 16, 735–743 (1998).

Article  CAS  PubMed  Google Scholar  * Kakui, H., Yamazaki, M., Hamaya, N.-B. & Shimizu, K. K. in _Pollen and Pollen Tube Biology: Methods and Protocols_ (ed. Geitmann, A.) (MIMB Humana

Press, 2020). * Tedder, A. et al. Evolution of the selfing syndrome in _Arabis alpina_ (Brassicaceae). _PLoS ONE_ 10, e0126618 (2015). Article  PubMed  PubMed Central  CAS  Google Scholar 

* Peterson, K. M. et al. _Arabidopsis_ homeodomain-leucine zipper IV proteins promote stomatal development and ectopically induce stomata beyond the epidermis. _Development_ 140, 1924–1935

(2013). Article  CAS  PubMed  PubMed Central  Google Scholar  * Vielle-Calzada, J. P., Baskar, R. & Grossniklaus, U. Delayed activation of the paternal genome during seed development.

_Nature_ 404, 91–94 (2000). Article  ADS  CAS  PubMed  Google Scholar  * Nordborg, M. et al. The pattern of polymorphism in _Arabidopsis thaliana_. _PLoS Biol._ 3, 1289–1299 (2005). Article

  CAS  Google Scholar  * Gautier, M. & Vitalis, R. rehh: an R package to detect footprints of selection in genome-wide SNP data from haplotype structure. _Bioinformatics_ 28, 1176–1177

(2012). Article  CAS  PubMed  Google Scholar  * Hu, T. T. et al. The _Arabidopsis lyrata_ genome sequence and the basis of rapid genome size change. _Nat. Genet._ 43, 476–481 (2011). Article

  PubMed  PubMed Central  CAS  Google Scholar  * Kim, S. et al. Recombination and linkage disequilibrium in _Arabidopsis thaliana_. _Nat. Genet._ 39, 1151–1155 (2007). Article  CAS  PubMed 

Google Scholar  * Tsuchimatsu, T. et al. Adaptive reduction of male gamete number in the selfing plant _Arabidopsis thaliana_. _Dryad_ https://doi.org/10.5061/dryad.jh9w0vt7z (2020).

Download references ACKNOWLEDGEMENTS We thank Matt Horton, Robert Martin, Atsushi Tajima, and Hirokazu Tsukaya for discussions, Karen Thomsen for sample preparation and pollen counting, and

the Functional Genomics Center Zurich for technical support. This study was supported by the Swiss National Science Foundation 31003A_182318, 31003A_159767, JST CREST Grant Number

JPMJCR16O3, Japan, to K.K.S., a Syngenta Fellowship through the Zurich-Basel Plant Science Center to T.T., T.S., and K.K.S., the University Research Priority Program of Evolution in Action

of the University of Zurich to K.K.S. and U.G., MEXT KAKENHI Grant Numbers 16H06469, 16K21727, 26113709 to K.K.S., 16H06467, 17H05833, 18H04813, 19H04851 to T.T., and 16H06465 to M.M.K.,

JSPS KAKENHI Grant Number 19K05976 to H.K., the Forschungskredit of the University of Zurich and an EMBO long-term fellowship to T.T., a JSPS postdoctoral fellowship for research abroad to

H.K. and T.T., and, in part, an Advanced Grant of the European Research Council to U.G. AUTHOR INFORMATION Author notes * These authors contributed equally: Takashi Tsuchimatsu, Hiroyuki

Kakui. AUTHORS AND AFFILIATIONS * Department of Evolutionary Biology and Environmental Studies & Zurich-Basel Plant Science Center, University of Zurich, 8057, Zurich, Switzerland

Takashi Tsuchimatsu, Hiroyuki Kakui, Misako Yamazaki & Kentaro K. Shimizu * Department of Plant and Microbial Biology & Zurich-Basel Plant Science Center, University of Zurich, 8008,

Zurich, Switzerland Takashi Tsuchimatsu, Hiroki Tsutsui, Afif Hedhly, Ueli Grossniklaus & Kentaro K. Shimizu * Gregor Mendel Institute, Austrian Academy of Sciences, Vienna BioCenter,

A-1030, Vienna, Austria Takashi Tsuchimatsu, Dazhe Meng & Magnus Nordborg * Department of Biology, Chiba University, Chiba, 263-8522, Japan Takashi Tsuchimatsu * Department of Biological

Sciences, University of Tokyo, Tokyo, 113-0033, Japan Takashi Tsuchimatsu * Kihara Institute for Biological Research, Yokohama City University, Yokohama, 244-0813, Japan Hiroyuki Kakui 

& Kentaro K. Shimizu * Graduate School of Science and Technology, Niigata University, Niigata, 950-2181, Japan Hiroyuki Kakui * Institute for Biochemistry and Biology, University of

Potsdam, 14476, Potsdam, Germany Cindy Marona & Michael Lenhard * Graduate School of Science, Nagoya University, Nagoya, 464-8602, Japan Hiroki Tsutsui & Masahiro M. Kanaoka * JST

ERATO Higashiyama Live-Holonics Project, Nagoya University, Nagoya, 464-8602, Japan Hiroki Tsutsui * Molecular and Computational Biology, University of Southern California, Los Angeles, CA,

90089-0371, USA Dazhe Meng * Department of Genomics and Evolutionary Biology, National Institute of Genetics, Mishima, Shizuoka, 411-8540, Japan Yutaka Sato * Institute of Integrative

Biology, ETH Zurich, 8092, Zurich, Switzerland Thomas Städler Authors * Takashi Tsuchimatsu View author publications You can also search for this author inPubMed Google Scholar * Hiroyuki

Kakui View author publications You can also search for this author inPubMed Google Scholar * Misako Yamazaki View author publications You can also search for this author inPubMed Google

Scholar * Cindy Marona View author publications You can also search for this author inPubMed Google Scholar * Hiroki Tsutsui View author publications You can also search for this author

inPubMed Google Scholar * Afif Hedhly View author publications You can also search for this author inPubMed Google Scholar * Dazhe Meng View author publications You can also search for this

author inPubMed Google Scholar * Yutaka Sato View author publications You can also search for this author inPubMed Google Scholar * Thomas Städler View author publications You can also

search for this author inPubMed Google Scholar * Ueli Grossniklaus View author publications You can also search for this author inPubMed Google Scholar * Masahiro M. Kanaoka View author

publications You can also search for this author inPubMed Google Scholar * Michael Lenhard View author publications You can also search for this author inPubMed Google Scholar * Magnus

Nordborg View author publications You can also search for this author inPubMed Google Scholar * Kentaro K. Shimizu View author publications You can also search for this author inPubMed 

Google Scholar CONTRIBUTIONS T.T., H.K., M.Y., T.S., M.L., M.N., and K.K.S. conceived and designed the study; T.T., H.K., M.Y., D.M., and K.K.S. analyzed the data with help from A.H., U.G.,

T.S., M.L., and M.N.; T.T., H.K., M.Y., C.M., H.T., A.H., Y.S., and M.M.K. performed experiments and generated the data; T.T., H.K., M.Y., and K.K.S. wrote the paper with inputs from U.G.,

T.S., M.L., and M.N.; All authors discussed the results and commented on the manuscript. CORRESPONDING AUTHOR Correspondence to Kentaro K. Shimizu. ETHICS DECLARATIONS COMPETING INTERESTS

The authors declare no competing interests. ADDITIONAL INFORMATION PEER REVIEW INFORMATION _Nature Communications_ thanks James Prendergast, and the other, anonymous, reviewer(s) for their

contribution to the peer review of this work. PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

SUPPLEMENTARY INFORMATION SUPPLEMENTARY INFORMATION REPORTING SUMMARY DESCRIPTION OF ADDITIONAL SUPPLEMENTARY FILES SUPPLEMENTARY DATA 1 SOURCE DATA SOURCE DATA RIGHTS AND PERMISSIONS OPEN

ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format,

as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third

party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the

article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright

holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Tsuchimatsu, T., Kakui, H.,

Yamazaki, M. _et al._ Adaptive reduction of male gamete number in the selfing plant _Arabidopsis thaliana_. _Nat Commun_ 11, 2885 (2020). https://doi.org/10.1038/s41467-020-16679-7 Download

citation * Received: 09 February 2018 * Accepted: 19 May 2020 * Published: 08 June 2020 * DOI: https://doi.org/10.1038/s41467-020-16679-7 SHARE THIS ARTICLE Anyone you share the following

link with will be able to read this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the Springer Nature

SharedIt content-sharing initiative