Four accessions of A. thaliana, Col-0 (N22625), Bor-4 (N22591), Mz-0 (N22636), and Uod-1 (N22612), as well as A. lyrata ssp. lyrata (CS22696), were used in this study. 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 with insecticide (Kendo Gold, Syngenta Agro). The frameshift mutants of RDP1 created by CRISPR/Cas9, i.e., rdp1-3 (Col-0) and rdp1-6 (Bor-4), were generated previously (Tsuchimatsu et al., 2020). The frameshift of rdp1-3 and rdp1-6 were identical. These frameshift mutants were genotyped by PRIMA using primers (Forward primer, 5′–TAGGCACAATGGAAAGTTAG–3′; Reverse primer, 5′–TT AACATAAAAGAACCATTGTAAG–3′) and 40-mer probe (5′–GGAGACTTTGTAGATACCAGAGTTCATCTTCTGCAGAAC G–3′) as described previously (Kakui et al., 2021). T-DNA fragments containing CRISPR/Cas9 and GFP marker were removed from each strain of rdp1-3 and rdp1-6 by selecting non-fluorescence seeds (Tsutsui and Higashiyama, 2017). Bolting time was determined when the main inflorescence stem reached 3 cm. Flowering time was recorded when the first flower opened. Time-lapse movie of plant growth from Col-0 and rdp1-3 was obtained using a TLC200 Pro camera (Brinno, Taiwan).

For histological analysis, inflorescences were fixed with FAA (formaldehyde:acetic acid:70% ethanol = 1:1:18), 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 color camera (DMC2900, Leica).

The RDP1 sequence of Col-0, including the 1,960 bp upstream region from the coding sequence to 661 bp downstream from the stop codon, was amplified using the following primers: 3795_At1g25260F (5′–TTTCTCCCCACATTTCTC–3′) and 3988_At1g25260R (5′–TATGTTATCAAAATTCATAAAATG–3′). The RDP1 sequences from different accessions were amplified from Mz-0, Col-0, Bor-4, and A. lyrata by PCR (PrimeSTAR GXL, Takara Bio, Japan). These PCR products were cloned into pFAST-R01 (Shimada et al., 2010). Each construct was independently transformed into rdp1-3 (Col-0 background) plants using the floral-dip method (Clough and Bent, 1998) with Agrobacterium tumefaciens (GV3101).

Arabidopsis thaliana plants with heterozygous RDP1 alleles from two accessions, i.e., the Col-0 accession (RDP1/rdp1-3, termed RDP1^Col/rdp1^Col hereafter) and the Bor-4 accession (RDP1/rdp1-6, termed RDP1^Bor/rdp1^Bor hereafter), were prepared (Tsuchimatsu et al., 2020). F₁ plants were generated by crossing RDP1^Col/rdp1^Col and RDP1^Bor/rdp1^Bor. Subsequently, F1 progenies are genotyped and analyzed.

We sampled flower buds and counted pollen numbers using a cell counter as described previously (Kakui et al., 2020; Tsuchimatsu et al., 2020). In summary, flowers of stage 12 [unopened anther with mature pollen grains (Smyth et al., 1990)] were collected and incubated overnight at 60°C. Flowers were collected from the side stem because flowers of the main stem have significantly larger pollen numbers (Tsuchimatsu et al., 2020). 1st and 2nd flowers of flower buds were excluded from sampling because they tend to show abnormal flower shapes (Kakui et al., 2020). We collected flowers from the start of flowering until approximately 2 weeks later. Subsequently, 30 μL of 5% Tween-20 was added, and the mixture was sonicated using a sonicator (Bioruptor Plus, Diagenode, Belgium), to release the pollen grains. The pollen suspension was mixed with a pollen counting solution (CASYton, OMNI Life Science, Germany), and particles were counted on a cell counter (CASY cell counter, OMNI Life Science). Pollen-number data were statistically analyzed and plots were constructed in R (R Core Team, 2013).

Total RNA was isolated from four replicates of flower buds of the wild-type and rdp1-3 plants using the RNeasy plant mini kit (Qiagen, Germany). The samples encompassed a wide range of developmental stages (flower stages 1–12) (Sanders et al., 1999), excluding opened flowers. Next, 75-bp single-end read sequencing was performed on a NextSeq 500 sequencer (Illumina, San Diego, CA, United States). The sample information and the number of reads are listed in Supplementary Table 1. The transcriptome data obtained were analyzed using the SUSHI framework (Hatakeyama et al., 2016). We defined DEGs based on a false discovery rate (FDR) < 0.1. All DEGs data are listed in Supplementary Table 2. The gene ontology (GO) enrichment analysis was performed using the ShinyGO v0.66 software (Ge et al., 2020).

The CRISPR/Cas9-generated frameshift rdp1-3 mutant in the Col-0 background produced about half the number of pollen grains compared with the wild-type counterpart (Tsuchimatsu et al., 2020). We further observed the section of the developing anther and confirmed that the number of microspores was already decreased in the mutant at the anther stage 8 (Sanders et al., 1999; Figures 1A,B). Tsuchimatsu et al. (2020) reported that the frameshift rdp1-3 mutants exhibited other pleiotropic phenotypes, such as reduced ovule number and delayed growth, in contrast with the effect of natural variants. Here, we found that the bolting and flowering times were significantly delayed (Wilcoxon rank-sum test; P = 4.33e–05 for bolting time, P = 1.08e–05 for flowering time, Figures 1C,D and Supplementary Video 1). These results support the notion that RDP1 functions as an Mrt4 homolog in A. thaliana, in accordance with the slower growth observed in the Mrt4 null mutant in yeast (Rodriguez-Mateos et al., 2009a).

To narrow down the causal sequence/region underlying the allelic effect of RDP1, we performed functional experiments using natural accessions. We first performed a traditional transgenic test to detect the allelic effect of RDP1 (Figure 2A), but found that the resolution was not sufficient to detect it, as follows. The prepared RDP1 sequences ranged from the upstream region (1,960 bp in Col-0) to the downstream region (661 bp in Col-0, Supplementary Figure 1) of three A. thaliana accessions (Col-0, Mz-0, and Bor-4), as well as A. lyrata ssp. lyrata. Mz-0 was among the accessions with the lowest pollen number, and the predominantly outcrossing species A. lyrata produced a much greater number of pollen grains than the A. thaliana accessions (Tsuchimatsu et al., 2020). Each sequence was introduced into the CRISPR/Cas9-generated rdp1 null (non-functional) mutant (rdp1-3/rdp1-3, Col-0 background), and more than 10 independent transformants were obtained for each of the four constructs. However, the variation in pollen number was high, even among the transformants of the same construct, and the range of the number of pollen grains largely overlapped between the four constructs (Supplementary Figures 2A,B). When these data were compared with those obtained for wild-type Col-0, three of them showed no significant difference, whereas, unexpectedly, the construct with Col-0 showed a slight increase in pollen number, although a previous complementation test on another rdp1 mutant background (rdp1-1) showed no significant difference compared with Col-0 (Tsuchimatsu et al., 2020). We speculate that the positional effect of the insertion site was too large to detect a subtle allelic difference. In Arabidopsis transgenic experiments, it is difficult to control the insertion sites of transgenes. Depending on the insertion position, the expression level, timing, and/or tissue type of the transgenes may be affected and other genes can be disrupted. Although a drastic increase in the number of independent transgenic lines may eventually result in significant differences, the construction of these lines and the measurement of phenotypes would be highly tedious.

We then used the quantitative complementation test to detect the allelic effect of RDP1. Previously, we showed the allelic differences of RDP1 between Uod-1 (a variant conferring reduced pollen number) and Bor-4 (a variant conferring increased pollen number) via the quantitative complementation test (Figure 2B; Tsuchimatsu et al., 2020). Here, we tested quantitative complementation using the Bor-4 and Col-0 accessions. We aligned the RDP1 genomic region of Col-0 to that of Uod-1 and Bor-4 (Figure 3A and Supplementary Figure 1). The analysis of the alignment of genomic sequences showed that Col-0 experienced a historical recombination event near the start of the transcribed region, i.e., between its upstream and transcribed (exon/intron) regions (Figure 3B; Tsuchimatsu et al., 2020). The upstream sequence of Col-0 was close to that of Bor-4. However, the exonic and intronic sequences of Col-0 were very close to those of Uod-1, and there was only a single amino acid substitution differentiating Bor-4 from the three other accessions (S222L, Figures 3B,C). This amino acid position was deduced to be functionally important because the corresponding serine residues at the C-terminal region of human Mrt4 regulate cellular localization under the stress condition (Michalec-Wawiorka et al., 2015). In addition, the 3′UTR and downstream regions of Col-0 had several unique substitutions compared with Uod-1 and Bor-4. Phenotypically, Col-0 and Uod-1 had smaller numbers of pollen grains (Col-0, 3,355 pollen grains/flower; Uod-1, 3,277 pollen grains/flower) in contrast with Bor-4 (4,528 pollen grains/flower) (Tsuchimatsu et al., 2020). The unique sequence character of Col-0 prompted us to study whether the coding region has a variant conferring reduced pollen number.

To compare the RDP1 allele, plants were prepared as follows (Figure 2B): (1) RDP1 heterozygous plants (RDP1/rdp1) were prepared from the Bor-4 and Col-0 accessions; (2) two heterozygous plants were crossed; and (3) RDP1^Bor/rdp1^Col and rdp1^Bor/RDP1^Col plants were selected. These plants had a single different functional RDP1 allele in the original chromosomal position, but all other genome sequences were identical: a single disrupted RDP1 allele with the same frameshift position and the same genome constitution (one chromosome from Bor-4 and the other from Col-0) outside of RDP1. We found that the plants with the functional Bor-4 allele of RDP1 (RDP1^Bor/rdp1^Col) had a significantly higher pollen number than those with a functional Col-0 allele (rdp1^Bor/RDP1^Col) (P = 1.10 × 10^–5; Figure 4). We estimated allelic effect of RDP1 based on median in the quantitative complementation test. RDP1^Bor is 1,599 pollen grain increasing [RDP1^Bor/rdp1^Col (4,528 pollen grain) vs. rdp1^Bor/rdp1^Col (2,929 pollen grain)] and RDP1^Col is 1,373 pollen grain increasing [rdp1^Bor/RDP1^Col (4,302 pollen grain) vs. rdp1^Bor/rdp1^Col (2,929 pollen grain)] (Supplementary Figure 3). This result elucidated the allelic difference between Bor-4 and Col-0.

Next, we examined which biological processes are disturbed by the disruption of the RDP1 gene in flower bud tissues, and whether the expression of ribosome-related genes is altered as previously described for other Arabidopsis ribosome mutants. We performed a transcriptome analysis using wild-type of Col-0 and rdp1-3 plants. A large number of DEGs between the wild-type and rdp1-3 plants was detected using the criterion of FDR < 0.1 (3,020 genes, Supplementary Table 2). Among them, 1,284 genes were upregulated and 1,736 genes were downregulated in rdp1-3 plants. A GO term analysis showed the enrichment of many ribosome-related genes in the upregulated gene set (Figure 5 and Table 1 and Supplementary Table 3). These results strongly support the notion that RDP1 of A. thaliana functions in the biogenesis of the ribosome as a yeast Mrt4 homolog (Lo et al., 2010; Michalec-Wawiorka et al., 2015; Greber, 2016; Saez-Vasquez and Delseny, 2019).

The GO term analysis also revealed the enrichment of genes relevant for pollen development, as well as pollen-related categories in the downregulated gene set (e.g., “Pollen development,” “Pollen wall assembly,” and “Pollen tube growth”; Figure 6 and Supplementary Table 4). Among the genes relevant for pollen development, three bHLH transcription factor genes, bHLH010, bHLH089, and AMS, were downregulated (Table 2), all of which interact with the DYT1 during anther development (Feng et al., 2012; Ferguson et al., 2017). Although their fold change itself was modest (29.7–34.4% reduction; Table 2), we also observed a systematic downregulation of genes that were directly regulated by AMS (19 out of 23, P = 2.901e–15, Fisher’s exact test; Table 2; Xu et al., 2014). These data suggest that the bHLH pathway and its downstream genes were downregulated in the rdp1-3 mutant. By contrast, we did not find any enrichment of GO terms related to ovule development, although ovule number was also reduced in the rdp1-3 (Tsuchimatsu et al., 2020). This expression analysis in the flower bud suggests a distinctive role for RDP1 in anther development through transcription factor regulation.

Although the quantitative complementation test has been validated in several animal and bacterial species, such as yeast, Drosophila, mosquito, fish, and mouse (Steinmetz et al., 2002; Stern, 2014; Turner, 2014; Podobnik et al., 2020), the preparation of mutants in a specific site of a gene represents a bottleneck regarding the application of quantitative complementation tests in plants (Stern, 2014). Recent progress in genome-editing technologies has allowed the efficient generation of frameshift non-functional mutants at the targeted gene in several different accessions. Thus, a quantitative complementation test will have a broader application of obtaining causal evidence of subtle allelic effects in plant species.

In this study, we performed a quantitative complementation test of RDP1 between Col-0 and Bor-4. A key of the method is to generate equivalent non-functional alleles in two different backgrounds by CRISPR/Cas9. We generated frameshift mutants at the same position, i.e., rdp1-3 and rdp1-6 in Col-0 and Bor-4, respectively. We consider them as null mutants because frameshift mutants are in general disruptive, and the yeast studies suggested the importance of C-terminal domains (Michalec et al., 2010), but we note that the quantitative complementation test can work even if they are not completely null. Combined with the previously reported results of the quantitative complementation test between Uod-1 and Bor-4 (Tsuchimatsu et al., 2020), our data support the notion that both Col-0 and Uod-1 have alleles that confer reduced pollen number against Bor-4. Although it is formally possible that more than two functional variants exist, the selective sweep signature of the long-haplotype variants suggests that a single mutation conferring reduced pollen number spread to many accessions, including Col-0 and Uod-1. This implies that the causal mutation of the reduced pollen number is shared between Col-0 and Uod-1, in contrast with Bor-4. The alignment of the genomic sequences of the three accessions revealed that such substitutions were mostly found in exon/intron regions, and not in upstream or downstream regions, reflecting the historical recombination inferred in Col-0. Although it is possible that mutations outside of the aligned region may have contributed to this phenomenon, substitutions in the exon/intron regions, including S222L, are promising candidate causal mutations of the reduced pollen number (Figures 3B,C), which warrants further experimental verification. In short, quantitative complementation tests using multiple accessions serve as a tool to narrow down the candidate regions of causal mutations of quantitative traits when transgenic experiments are not sufficiently powerful.

Remarkably, we detected two ontologies related to the ribosomal large subunit (“Ribosomal large subunit assembly” and “Ribosomal large subunit biosynthesis”) in the transcriptome analysis of rdp1-3. We found that 23.0% of the genes involved in the ribosomal large subunit were significantly differentially expressed (33 out of 143 genes), with most of them being upregulated (32 out of 33 genes) (Table 1). Conversely, only 8.2% of the genes involved in the ribosomal small subunit were significantly different (8 out of 98 genes). Furthermore, S222L was the only non-synonymous mutation in the mapped region revealed by the two complementation tests; thus, it was a candidate causal mutation for the pollen-number phenotype (as discussed above). Moreover, this amino acid position was shown to be important for the cellular localization and ribosomal function of human RDP1 (Michalec-Wawiorka et al., 2015). Taken together, our results strongly support the contention that RDP1 of A. thaliana functions in the biogenesis of the ribosomal large subunit as a yeast Mrt4 homolog (Lo et al., 2010; Michalec-Wawiorka et al., 2015; Greber, 2016; Saez-Vasquez and Delseny, 2019).

The mechanisms underlying where and how RDP1 affects the expression of bHLH transcriptional factor genes, such as AMS, bHLH010, and bHLH089, remain unknown. Previous studies have reported that the spatial expression pattern of these transcriptional factor genes and that of RDP1 overlap partly in the same cell types (Zhu et al., 2011; Zhu et al., 2015; Tsuchimatsu et al., 2020). Real-time PCR, promoter-GUS assay, and in situ hybridization studies have suggested that RDP1 is expressed broadly and strongly in proliferating cells, including sporogenous cells (Tsuchimatsu et al., 2020). A transcriptome study of dissected tissues also demonstrated that RDP1 is expressed in both tapetal cells and the pollen lineage (microspore, bicellular, and tricellular stages) (Honys and Twell, 2004; Li et al., 2017). The expressions of bHLH010 and bHLH089 were detected in both microspores and tapetum cells via in situ hybridization (Zhu et al., 2015). Moreover, in situ hybridization suggested that AMS is expressed in both tapetum cells and microspores (Zhu et al., 2011), although the defect of ams mutants is primarily found in surrounding tapetum cells that provide pollen wall to the germ lines (Xu et al., 2010). These expression data suggest that these bHLH genes are co-expressed with RDP1 in tapetum and microspore cells. It is reported that the feedback regulation of bHLH transcription factor genes (DYT1, bHLH010, bHLH089, and bHLH091) occurs during anther development (Cui et al., 2016). Because this feedback regulation requires the synthesis of the protein encoded by each gene, RDP1 may contribute to the expression of these genes through translation by ribosomes. This hypothesis was supported by the observation that two of these genes (bHLH010 and bHLH089) were downregulated in the rdp1-3 mutant. We also considered the possibility that the transcriptional reduction of these bHLH genes in tapetal cells has an indirect effect on tapetum development via the microspore reduction afforded by RDP1 mutation. However, this is unlikely because the expression patterns of other genes that are essential for tapetum development, such as DYT1, TPD1, and MS188, were not affected in rdp1-3. We could not exclude the possibility that the aberrant translational machinery of rdp1-3 may affect the mRNA abundance, or of downstream transcription independent from the feedback (Tiruneh et al., 2013; Beine-Golovchuk et al., 2018; Cheong et al., 2021) or multi-functionalization of RDP1 independent from the ribosome. Further detailed analyses are required to evaluate this scenario. In addition, previous transcriptome data obtained for the bhlh triple mutant (bhlh010/bhlh089/bhlh091) (Zhu et al., 2015) did not show an obvious transcriptional difference for ribosomal genes, including RDP1, suggesting that RDP1 or ribosome function is not located downstream of the bHLH gene regulatory network.

The ribosome is a protein biosynthesis machinery that is present in all living cells. Although the ribosome has been considered a solid component common to all cells, organ-specific phenotypes are revealed from mutant analysis of ribosome and ribosome-biogenesis genes in yeast, humans, and plants. This ribosome specialization is now suggested as not being a rare phenomenon (Martinez-Seidel et al., 2020; Norris et al., 2021). Reproductive defects of ribosomal mutants were observed, such as reduced stamen number, smaller pollen size, slower pollen tube elongation, and reduced functional ovules (Stirnberg et al., 2012; Luo et al., 2020). RDP1 is considered a homolog of yeast Mrt4, which is a key component for the 60S ribosomal assembly. We found that the expression level of RDP1 varied among different cell types and was higher in proliferating cells (Tsuchimatsu et al., 2020). Natural variants under selection supported the critical role of RDP1 in adaptive evolution (Tsuchimatsu et al., 2020). Furthermore, a non-functional RDP1 mutant exhibited reduced pollen number and slower growth speed. These data support the idea that the high expression level of RDP1 is required to respond to the high demand of ribosome production in proliferation cells. Interestingly, Tsuchimatsu et al. (2020) reported that quantitative complementation test with Uod-1 and Bor-1 showed allelic differences for pollen number but not growth speed (rosette size and bolting/flowering date). This result may suggest the ribosome specialization for pollen number. At first glance, this may be similar to cytoplasmic sterility, where the reduction of mitochondrial function results in pollen-specific phenotype, where high energy consumption may be necessary (Warmke and Lee, 1978; Chase, 2007). However, our transcriptome data does not support this idea because there was no GO term “mitochondria” in down-regulated genes. Further analysis of RDP1 will provide clues regarding how ribosome specialization affects the pollen number regulation.