The three F₁ populations, Okitsu No. 46 × ‘Kara’ (O46-K), ‘Sweet spring’ × Okistu No. 56 (SS-O56) and ‘Harehime’ × Okistu No. 63 (H-O63), and the varieties/selected strains used in this study were maintained in the Division of Citrus Research, Institute of Fruit Tree and Tea Science, NARO (Shizuoka, Japan) ( Tables S1 , S2 ). The individuals of the F₁ populations were grafted onto trifoliate orange rootstocks as a single replicate in April 2012 and 2013 ( Table S1 ), while the NPG of individuals in the O46-K population had been previously evaluated in 2015 ( Table S1 ) (Goto et al., 2016; Goto et al., 2018). The NPG in the varieties/selected strains in 2017 and 2018, O46-K in 2016, SS-O56 in 2015 and 2016, and H-O63 in 2015 and 2016 was evaluated. The protocols used for NPG evaluations are available at protocols.io (dx.doi.org/10.17504/protocols.io.q78dzrw).

Genomic DNA was extracted from fresh leaves of the F₁ populations and the 85 varieties/selected strains through a modified protocol using cetyltrimethylammonium bromide and a high-salt precipitation solution (1.2 M NaCl, 0.8 M sodium citrate; protocol available at dx.doi.org/10.17504/protocols.io.dm6gpj8jpgzp/v1). Six DNA markers (Marker No. 1, 9, 12, 13, 15, 19) in the MS-P1 locus were used in a previous report ( Table S3 ) (Ollitrault et al., 2010; Shimizu et al., 2016), while 13 new SSR markers (Marker No. 2–8, 10, 11, 14, 16–18) in the MS-P1 locus ( Table S3 ) were identified for this study. The Clementine genome sequence between TSRF161 and SSR08B66, which corresponds to scaffold 8, 5,153,769−19,761,190 bp in the physical map of C. clementina genome v1.0 (JGI) (Goto et al., 2018), was obtained from Phytozome (https://phytozome-next.jgi.doe.gov/info/Cclementina_v1_0) (Wu et al., 2014). The corresponding scaffolds of Satsuma mandarin from MiDB (https://mikan.dna.affrc.go.jp/) were identified using BLAST (https://mikan.dna.affrc.go.jp/blast/) (Kawahara et al., 2020). The identified scaffolds were screened for dimeric and trimeric SSR sequences with ≥ 7 repeats using the Simple Sequence Repeat Identification Tool (https://archive.gramene.org/db/markers/ssrtool) (Temnykh et al., 2001). The primers for the 13 markers identified based on the detected SSRs were designed using Primer3 v.4.1.0 (https://bioinfo.ut.ee/primer3/) ( Table S3 ) (Kõressaar et al., 2018). Dimeric SSRs were preferred for developing markers that can detect regions with high polymorphism, and the multicolored post-labeling method was used for genotyping analysis (Shimizu and Yano, 2011), as described previously (Shimizu et al., 2016).

Genotyping was performed with 19 markers (Marker No. 1–19 in Table S3 ) for 34 individuals in O46-K, 31 individuals in SS-O56, and 50 individuals in H-O63 ( Table S1 ). Subsequently, recombinant individuals within the MS-P1 locus were identified based on the genotype segregation pattern. A graphical genotype of the MS-P1 locus was constructed for the physical map of C. clementina genome v1.0 between 5,153,708 and 19,761,190 bp in scaffold 8. The male sterile genotypes within the MS-P1 locus in recombinant individuals were identified through comparisons with the genotype of male sterile varieties/selected strains ( Figure S1 ), whereas the male sterile phenotypes of the recombinant individuals were evaluated through comparison of the anthers from the male sterile and fertile varieties/selected strains ( Figure S1 ). The MS-P1 region was fine mapped considering the associations between the male sterile phenotypes and genotypes in the MS-P1 locus of recombinant individuals.

The gene locus, annotations, and protein sequences at the MS-P1 region of C. clementina genome v1.0 (JGI) were obtained from Phytozome (https://phytozome-next.jgi.doe.gov/info/Cclementina_v1_0) (Wu et al., 2014). The DNA sequences of maize (Zea mays) Rf2 (U43082), rice (Oryza sativa) Rf2 (AB583697), rice Rf17 (Os04g0475900), and sugar beet (Beta vulgaris) Rf1 (AB646135), all of which belong to the non-PPR family of Rf, were used as the query to perform BLAST analysis against the sequences in the Mikan Genome DB under default settings. The protein sequences of the PPR family at the MS-P1 region were analyzed using TargetP v.2.0 (https://services.healthtech.dtu.dk/service.php?TargetP-2.0) (Almagro Armenteros et al., 2019), Predotar v.1.04 (https://urgi.versailles.inra.fr/predotar/) (Small et al., 2004), and MitoProt II (https://ihg.helmholtz-muenchen.de/ihg/mitoprot.html) (Claros and Vincens, 1996) to predict their subcellular localization and targeting sequence. To classify the PPR subfamily, PPR protein sequences were analyzed using PPRfinder v.507b0fb (https://ppr.plantenergy.uwa.edu.au/) (Gutmann et al., 2020). Protein sequence alignments were performed using MAFFT v.7 (https://mafft.cbrc.jp/alignment/server/) (Katoh et al., 2019), while PPR motifs were identified using ScanProsite release 20.0 (https://prosite.expasy.org/scanprosite/), and phylogenetic trees were constructed using the neighbor-joining method of MEGA v.10.0.5 (https://www.megasoftware.net) with protein sequences of rice (Oryza sativa) Rf1a (DQ311053), rice Rf1b (DQ311054), rice Rf4 (KJ680249), Sorghum (Sorghum bicolor) Rf1 (Klein et al., 2005), Chinese cabbage (Brassica napus) Rfp1 (KX671967), Petunia (Petunia hybrida) Rf-PPR592 (AY10271), Ciclev10030242m, and Ciclev10030361m. Bootstrap values were calculated through a 1,000-permutation test.

Bulked stamen from a sterile selected strain (KyOw14) and fertile variety (‘Shiranuhi’) were collected seven days before flowering (DBF) and 1 DBF in a field from the Division of Citrus Research, Institute of Fruit Tree and Tea Science, NARO (Shizuoka, Japan), as two biological replicates. Total RNA was isolated using an RNeasy Plant Mini Kit (Qiagen, Hilden, Germany), while the RNA-seq analysis was carried out by Hokkaido System Science Co., Ltd. (Sapporo, Japan). RNA-seq libraries were generated from the total RNA using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England BioLabs, Inc., Ipswich, MA). The 150 bp paired-end sequencing of RNA-seq libraries was performed using the NovaSeq 6000 system (Illumina, San Diego, CA), and the sequenced reads were trimmed using fastp v.0.22.0 and mapped to C. clementina genome v1.0 (JGI) (https://phytozome-next.jgi.doe.gov/info/Cclementina_v1_0) using HISAT2 v.2.1.0. Transcripts per million (TPM) values were obtained to measure gene expression using StringTie v.2.1.1. TPM values were imported to Subio Platform v.1.24.5853 (Subio Inc., Aichi, Japan) for normalization. Statistical analysis of the comparisons between ‘Shiranuhi’ and KyOw14 was carried out using the “compare 2 groups” tool of the Subio platform. All RNA-seq data was deposited in the DDBJ Sequence Read Archive under accession number DRA015326.

To confirm the reproductivity of the transcriptome analyses concerning the PPR genes localized in the mitochondria, three replicates of the total RNA from the stamen of KyOw14 and ‘Shiranuhi’ at 7 DBF and 1 DBF were used for quantitative RT-PCR and semi-quantitative RT-PCR. Reverse-transcription was carried out using QuantiTect Reverse Transcription Kit (Qiagen, Hilden, Germany). Quantitative PCR was performed QuantStudio3 Real Time PCR System (Thermo Fisher Scientific, Waltham, USA) using Power SYBR Green PCR Master Mix (Thermo Fisher Scientific) under 10 min at 95°C, followed by 40 cycles of 15 s at 95°C and 60 s at 60°C. The primers for the gene expression of Ciclev10030242m (Fd: GATAAGAGAAATGTAATGCCAGACG, Rv: CCTTCACACGTTGTAATAAGAATGG), Ciclev10030361m (Fd: AGAGAAATGTAATGCCAGACGC, Rv: CAGCTCATTGACTTACTTGTGTCTA), and Ciclev10029947m (Fd: ACCCAATTGTGTCATATTTACTACGC, Rv: AATCAACATTCAACTCACCCACTTAC) were designed at a specific region or 3’UTR referred by Phytozome (https://phytozome-next.jgi.doe.gov/info/Cclementina_v1_0) (Wu et al., 2014). The expression of these genes were normalized by that of eEF1a (Fd: TCAAGGATCTCAAGCGTGGTT, Rv: CTTCCCTGGCCGGATCAT) (Shimada et al., 2018). Specific amplification in each PCR reaction was confirmed by melting curve analysis and agarose gel electrophoresis.

The sources of cytoplasm in the 85 varieties/selected strains were determined referring to the pedigree chart (Imai et al., 2017).

Both 00918-2 and 00918-3 were designed on the scaffolds of C_unshiu_00918, as described above (Marker No. 20 and 21 in Table S3 ). Haplotypes were defined based on segregation of the genotyped alleles at 00918-1, 00918-2, and 00918-3 in the three F₁ populations (O46-K, SS-O56, and H-O63) ( Table S1 ), before being classified into seven groups (HT1–HT7). HT3, derived from sweet orange, and HT6, derived from Dancy tangerine, Iyo, Ponkan, or Willowleaf mandarin, were differentiated based on the pedigree chart (Imai et al., 2017), although HT3 and HT6 are indistinguishable by 00918-1, 00918-2 and 00918-3. Additional haplotypes were defined and classified into four groups (HT8–HT11) based on the pedigree chart (Imai et al., 2017) of genotyped alleles, with 00918-1, 00918-2, and 00918-3 in 85 varieties/selected strains ( Table S2 ). Finally, the diplotype (combination of haplotypes) in the 85 varieties/selected strains and F₁ population individuals were determined using the 11 haplotypes and pedigree chart (Imai et al., 2017).

Statistical analyses were performed with EZR v.1.37, which is a graphical user interface for Rcmdr v.2.4-0 (https://www.jichi.ac.jp/saitama-sct/SaitamaHP.files/statmedEN.html) (Kanda, 2013).

Previously, we identified an MS-P1 locus between 33.5cM (TSRF161) and 69.1cM (SSR08B66) on linkage group 8 (LG8) in the corresponding linkage map of Okitsu No. 46 × Okitsu No. 56 (Goto et al., 2018). For identifying individual recombinations within the MS-P1 locus, genotyping was carried out for 19 DNA markers between TSRF161 and SSR08B66 (Marker No. 1−19; Table S3 ) in individuals of three F₁ populations derived from SS-O56, O46-K, and H-O63 ( Table S1 ). The following three recombinant individuals were identified: Recombinant line 1, male fertile; Recombinant line 2, male sterile; and Recombinant line 3, male sterile ( Figure 1A ). Several DNA marker genotypes of these recombinant lines indicated the presence of the same genotypes as that in Okitsu No. 46, Satsuma mandarin, and ‘Kiyomi’ ( Figure 1A ; Table S2 ). Notably, there were no pollen grains on the anthers of these varieties/selected strains, indicative of male sterility ( Figure S1 ); therefore, the corresponding genotypes were considered male sterile. Association analysis between the male sterile phenotype and genotypes in the recombinant lines enabled isolation of the MS-P1 region within 00220-2 and 00432-1, corresponding to a 920kb genomic region from 6,184,904 to 7,104,658 bp of scaffold 8 of C. clementina genome v1.0 (JGI) ( Figure 1A ).

Based on the annotation information for C. clementina genome v1.0 (JGI), 66 functional genes, including 10 PPR family genes, were predicted on the 920 kb genomic region, corresponding to the newly refined MS-P1 region ( Figure 1B ). The PPR family genes were highly redundant in this genomic region, and a tandem repeat-like cluster was observed ( Figure 1B ; Table 1 ). Genes homologous to maize Rf2 (U43082), rice Rf2 (AB583697), rice Rf17 (Os04g0475900), or sugar beet Rf1 (AB646135), which are non-PPR families of Rf, were not found among the 66 genes here. Bioinformatic analysis using PPRfinder revealed that all PPR family genes at the MS-P1 region belonged to the P-class subfamily ( Table 1 ). The proteins encoded by three genes in the PPR family, Ciclev10030242m, Ciclev10029947m, and Ciclev10030361m, were predicted to localize in the mitochondria using TargetP 2.0 (Almagro Armenteros et al., 2019), Predotar v.1.04 (Small et al., 2004), or MitoProt II (Claros and Vincens, 1996) ( Table 1 ).

Assuming that male sterility in the Kishu-cytoplasm fits into the CMS genetic model, the Rf corresponding to the Kishu-cytoplasm (Rf-MS-P1) should be expressed more in a male fertile variety/selected strain than in a male sterile variety/selected strain. Transcriptomic analyses of KyOw14 (sterile selected strain) and ‘Shiranuhi’ (fertile variety) stamens were performed at 7 and 1 DBF through RNA-seq analysis. RNA extracted from the stamens (two replicates) were used to build libraries for RNA-seq ( Figure 2A ). An average of 26 million reads per sample were obtained ( Table S4 ), among which 97% of the reads were passed through quality thresholds, and 95% were mapped to the reference genome ( Table S4 ). These data showed that the RNA sequencing quality was highly valid and suitable for comparative analysis. Comparative transcriptomic analysis between the PPR family genes at the MS-P1 region of KyOw14 and ‘Shiranuhi’ was performed ( Figure 2B ; Table S5 ). The expression of Ciclev10030242m, Ciclev10029914m, and Ciclev10030361m was significantly higher in ‘Shiranuhi’ than in KyOw14 at both 7 and 1 DBF (p< 0.05) ( Figure 2B ; Table S5 ). Among them, Ciclev10030242m and Ciclev10030361m were predicted to be localized in the mitochondria ( Table 1 ). Expression of Ciclev10029947m, which was also predicted to be localized in the mitochondria ( Table 1 ), was hardly detected in either KyOw14 or ‘Shiranuhi’ ( Figure 2B ; Table S5 ). The expression of Ciclev10030242m and Ciclev10030361m was reproduced by quantitative RT-PCR ( Figure S3A ). The expression of Ciclev10029947m was also reproduced by quantitative RT-PCR (( Figure S3A ) and semi-quantitative RT-PCR ( Figure S3B ). Predicted subcellular localization and transcriptomic analysis indicated Ciclev10030242m and Ciclev10030361m as transcriptionally plausible candidates for Rf-MS-P1.

The protein sequence analysis showed that Ciclev10030242m and Ciclev10030361m were highly similar proteins, with 86.66% overlap ( Figure 3A ), maintaining 14 PPR motifs with a mitochondrial targeting sequence ( Figure 3A ). Neighbor-joining phylogenetic tree analysis was carried out using protein sequences of Ciclev10030242m, Ciclev10030361m, rice Rf1a (DQ311053), rice Rf1b (DQ311054), rice Rf4 (KJ680249), sorghum Rf1 (Klein et al., 2005), Chinese cabbage Rfp1 (KX671967), and petunia Rf-PPR592 (AY10271), all of which are PPR family Rfs. The result showed that Ciclev10030242m and Ciclev10030361m were clustered together ( Figure 3B ). Notably, the cluster was closer to that of Chinese cabbage Rfp1 and petunia Rf-PPR592 than to that of rice Rfs ( Figure 3B ). Further, sorghum Rf1 was far from the other clusters because it alone was included in the PLS-class subfamily ( Figure 3B ) (Dahan and Mireau, 2013). Therefore, Ciclev10030242m and Ciclev10030361m were potential candidates for Rf-MS-P1, as they belonged to the P-class PPR family, had mitochondrial targeting sequences, and were expressed more in ‘Shiranuhi’ (fertile variety) than in KyOw14 (sterile strain). In addition, the existence of Ciclev10030242m and Ciclev10030361m at the MS-P1 region provides evidence that molecular mechanisms underlying citrus male sterility fit into the general plant CMS model.

The molecular mechanism of citrus male sterility fits within the general plant CMS model; however, it remains unclear how the wide phenotypic variations of male sterility and fertility observed in Japanese breeding germplasms are controlled. It was hypothesized here that male sterility is influenced by a combination of allelic haplotypes in which Rf or rf are located; therefore, we defined allelic haplotypes at the MS-P1 region among the founders of Japanese breeding germplasm. The nearest markers to candidates for Rf-MS-P1 were narrowed down to 00918-1 through fine mapping of the MS-P1 region ( Figure 1A ). We focused on the region containing 00918-1 and developed two markers, 00918-2 and 00918-3 in the proximity of 00918-1 using the genome sequence of C_unshiu_00918 in MiDB ( Figure 1 ; Marker No. 20 and 21 in Table S3 ). Each allelic haplotype was determined based on the genotype of these three markers. The 85 varieties/selected strains ( Tables S2 ; S6A ) and the three F₁ populations (SS-O56, O46-K, H-O63) were genotyped ( Tables S1 ; S6B ) with the 00918-1, 00918-2, and 00918-3 markers. The parent-offspring diagnosis did not show discrepancies in the obtained genotypes for any of the varieties/selected strains and the three F₁ populations except for ‘Willking’ (see legend in Table S6A ). A total of 11 haplotypes were defined (HT1–HT11) among the 12 founder varieties, with three genotyped alleles at 00918-1, 00918-2, and 00918-3 ( Table 2 ). The reduced NPG was linked to the haplotype block derived from kunenbo in the MS-P1 locus (Goto et al., 2018), which was identified as HT1 in this study ( Table 2 ). Okitsu No. 46 harbored HT1 and HT3 derived from kunenbo and sweet orange, respectively ( Table S6A ). Okitsu No. 46 showed reduced NPG ( Figure S1 ; Table S6A ), and individuals with reduced NPG appeared in the progeny of Okitsu No. 46 (Goto et al., 2016; Goto et al., 2018); therefore, HT3 was associated with the reduced NPG and was defined as a haplotype that harbored the “248” allele at 00918-1, the “234” allele at 00918-2, and the “221” allele at 00918-3 derived from sweet orange ( Table 2 ). HT6 was defined as a haplotype comprising same alleles as HT3 derived from Dancy tangerine (C. tangerina hort. ex Tanaka), Iyo (C. iyo hort. ex Tanaka), Ponkan (C. reticulata Blanco), or Willowleaf mandarin (C. deliciosa Ten.) ( Table 2 ). The three alleles comprising each haplotype were completely linked to each other in all varieties/selected strains and individuals in the three F₁ populations, indicating that they were inherited as a haplotype blocks. The diplotype in the 85 varieties/selected strains and the three F₁ populations was determined based on the 11 haplotypes ( Tables S6A, B ). HT3 and HT6 were distinguished based on the pedigree chart (Imai et al., 2017); however, several varieties/selected strains were indistinguishable from each other ( Table S6A )

Evaluating the appearance frequency of haplotypes at the MS-P1 region showed that HT1 was the most frequent haplotype, followed by HT6 and HT2 among the 85 varieties/selected strains ( Table 2 ). Japanese breeding programs primarily use Satsuma mandarin and ‘Kiyomi’ as seed parents ( Figure S2 ) (Imai et al., 2017), both of which have HT1, which is derived from kunenbo. ‘Kiyomi’ also has HT2, which is derived from sweet orange ( Figure S2 ); thus, HT1 or HT2 was present in all varieties used in the present study during 1st and 2nd generation breeding ( Figure S2 ). This constitutes one of the main reasons why HT1 and HT2 were frequently observed haplotypes. The diplotypes in the varieties/selected strains with Kishu-cytoplasm predominantly showed combinations of HT1 and HT2 ( Table S6A ); therefore, focus was placed on HT1 and HT2 at the MS-P1 region, and the association between the NPG and combinations of HT1/HT2 as well as another haplotype in the F₁ populations of SS-O56, O46-K, and H-O63 and the varieties/selected strains harboring the Kishu-cytoplasm (breeding germplasms) were investigated ( Tables S6A, B ). Subsequently, NPGs in the breeding germplasms were evaluated for two years ( Tables S6A, B ). First, to investigate the association between the NPG and the combination of HT1 and another haplotype, breeding germplasms which have HT1/HT1, HT1/HT2, HT1/HT3, HT1/HT4, HT1/HT5, HT1/HT6, or HT1/HT7 were selected based on the diplotype at the MS-P1 region ( Figure 4A ). Association analysis revealed that the NPGs in HT1/HT1, HT1/HT2, HT1/HT3, HT1/HT4, and HT1/HT5 were significantly lower than those in HT1/HT6 and HT1/HT7 (p< 0.05) ( Figure 4A ). Near zero pollen grains were detected in HT1/HT1 ( Figure 4A ), where the average NPG (7 pollen grains) was lower than that in HT1/HT2 (138), HT1/HT3 (1162), HT1/HT4 (636), and HT1/HT5 (688) ( Figure 4A ). Additionally, the NPG in HT1/HT2 was lower than that in HT1/HT3, HT1/HT4, and HT1/HT5 ( Figure 4A ). The NPG in HT1/HT6 and HT1/HT7 was relatively high, although they harbored HT1. Two varieties, ‘Nishinokaori’ and ‘Kara’, deviated from these results ( Figure 4A , white and black arrows, respectively; Table S6A ). Recombination was not observed between 00220-2 and 00432-1, which were flanking markers for 00918-1, -2, and -3 at the MS-P1 region ( Figure 1 ) of ‘Nishinokaori’ and ‘Kara’ (data not shown). Further, to investigate the association between the NPG and the combination of HT2 and another haplotype, the breeding germplasms that have HT2/HT2, HT2/HT3, HT2/HT5, HT2/HT6, or HT2/HT7 were selected ( Figure 4B ). Here, HT1/HT1 was used as the control. Association analysis showed that the average NPG in HT2/HT2 (922) was significantly lower than that in HT2/HT6 (4349) and HT2/HT7 (5466) (p < 0.05) ( Figure 4B ). Second, the average NPG in HT2/HT2 was higher than that in HT1/HT1, though this difference was not significant. ‘Asumi’ harbored HT2/HT3, with an average NPG of 4079; whereas Kuchinotsu No. 44 harbored HT2/HT5, and its average NPG was 7216 ( Figure 4B ; Table S6A ). Third, to investigate the association between NPG and combinations of haplotypes without HT1 and HT2, breeding germplasms that have HT3/HT4, HT5/HT5, and HT5/HT10 were selected, while HT1/HT1 and HT2/HT2 were used as the controls ( Figure 4C ). The average NPG in HT3/HT4 (8383) and HT5/HT5 (7148) was significantly higher than that in HT1/HT1 (p < 0.05) ( Figure 4C ); whereas the average NPG in HT3/HT4 and HT5/HT5 was higher than that in HT2/HT2 ( Figure 4C ). Kishu mandarin harbored HT5/HT10, and its average NPG was 4070 ( Figure 4C ; Figure S6A ). These results reveal that male sterility in Kishu-cytoplasm was influenced by a combination of allelic haplotypes at the MS-P1 region. In addition, HT1 and HT2 reduced NPG levels, while HT6 and HT7 increased NPG.

Proteins encoded by the major Rfs contain PPR motifs, localize to the mitochondria, and belong to the P-class subfamily (Lurin et al., 2004; Gaborieau et al., 2016). Ciclev10030242m, Ciclev10029947m, and Ciclev10030361m were annotated to the P-class PPR family genes, and their protein sequences predicted localization to the mitochondria ( Table 1 ). Ciclev10030242m and Ciclev10030361m were expressed more in the male fertile variety, ‘Shiranuhi’ (which harbored HT1/HT6), than in the male sterile selected strain, KyOw14 (which harbored HT1/HT1) ( Figure 2B ; Figure S3A ; Table S6A ). KyOw14 was obtained from the ‘Kiyomi’ × Satsuma mandarin cross and ‘Shiranuhi’ from ‘Kiyomi’ × Ponkan ( Table S2 ). Therefore, KyOw14 and ‘Shiranuhi’ have the same female parent and their male parents, Satsuma mandarin and Ponkan are classified as admixture mandarins based on genomic and phylogenetic analysis (Wu et al., 2018). These show that the genetic background of KyOw14 and ‘Shiranuhi’ are very close. In addition, the diplotype at MS-P1 in KyOw14 was HT1/HT1, which is homozygous non-functional restorer-of-fertility, and that in ‘Shiranuhi’ was HT1/HT7, which is non-functional restorer-of-fertility/functional restorer-of-fertility ( Table 2 ; Table S6 ). These data show that KyOw14 and ‘Shiranuhi’ are the optimal combination for transcriptome comparison in the genes at the MS-P1 region, suggesting that Ciclev10030242m and Ciclev10030361m are certainly functional in the fertile individual. Taken together, these observations show that Ciclev10030242m and Ciclev10030361m are the most likely potential candidates for Rf-MS-P1. Moreover, three nuclear genes that encode PPR proteins contribute to the cytonuclear interactions in citrus at Chr 3, Chr 4, and Chr 7 (Wang et al., 2022). In contrast, the candidates for Rf-MS-P1 identified in this study were located on Chr 8; therefore, to the best of our knowledge, this study is the first to identify potential candidates for Rf in citrus.

Ciclev10030242m and Ciclev10030361m were identified as candidates for Rf-MS-P1. The presence of these genes at the MS-P1 region provides evidence that male sterility in the Kishu-cytoplasm fits the CMS genetic model. In this model, the presence or absence of expressed male sterility is explained based on the normal cytoplasm, sterile cytoplasm, rf, and Rf (Chen and Liu, 2014). Male fertility is exhibited in individuals with normal cytoplasm, regardless of the nuclear gene. In addition, male sterility is expressed in individuals with male sterile cytoplasm and homozygous rf; therefore, even if individuals have male sterile cytoplasm, male sterility is suppressed, and male fertility is expressed in cases harboring Rf. A previous study has shown that heterozygosity for rf and semi-functional Rf partially restore male sterility, resulting in the expression of partial male sterility; however, homozygous individuals for the semi-functional Rf completely restored male sterility (Arakawa et al., 2019). The present study showed that the diplotype at the MS-P1 region influenced the NPG of breeding germplasms, as indicated by the association analysis ( Figures 4A–C ). The NPG was negligible in those harboring homozygous HT1 at the MS-P1 region, indicating that these individuals exhibited complete male sterility ( Figure 4A ). Considering the CMS genetic model, this result shows that HT1 was the rf haplotype ( Table 2 ). In contrast, > 2000 pollen grains per anther were detected in almost all individuals harboring HT1/HT6 and HT1/HT7 ( Figure 4A ), indicative of their exhibited male fertility. Considering HT1 as the rf haplotype, this result shows that HT6 and HT7 are the Rf haplotypes ( Table 2 ). In addition, the average NPG was 138 for individuals harboring HT1/HT2 and 922 for those harboring HT2/HT2 ( Figures 4A, B ), showing that HT2 was the less-functional Rf haplotype ( Table 2 ). NPG values > 2000 were detected in individuals harboring HT2/HT6 and HT2/HT7, supporting that HT6 and HT7 are Rf haplotypes ( Figure 4B ). Furthermore, the NPG range in individuals harboring HT1/HT3, HT1/HT4, and HT1/HT5 ( Figure 4A ) was 636−1162, indicating that they exhibited partial male sterility. Approximately 2000−15,000 NPG were detected in individuals harboring HT2/HT3, HT2/HT5, HT3/HT4, and HT5/HT5 ( Figures 4B, C ). Considering the results of Arakawa et al. (Arakawa et al., 2019), findings from the present study suggest that HT3, HT4, and HT5 are semi-functional Rf haplotypes ( Table 2 ). Comparatively, the function of restorer-of-fertility in HT8, HT9, HT10, and HT11 could not be determined owing the lack of sufficient data ( Table 2 ). The physiological and diplotype analysis together provided enough evidence to support that the molecular mechanism underlying male sterility in the Kishu-cytoplasm fits the CMS genetic model. In addition, these results also show that the NPG in the breeding germplasms were influenced by the combination of allelic haplotypes at the MS-P1 region.

The results here also suggest that Rf-MS-P1 exists in a genome region including HT6 or HT7, which are the Rf haplotypes ( Table 2 ). A region on C. clementina genome v1.0 (JGI) corresponding to the MS-P1 region was derived from Willowleaf mandarin, which harbors HT6 ( Table 2 ) (Ollitrault et al., 2012); therefore, the genome region must contain Rf-MS-P1. Accordingly, the use of the sequence of C. clementina genome v1.0 (JGI) to identify Rf-MS-P1 candidates was a reasonable strategy.

‘Nishinokaori’ and ‘Kara’ exhibited a discrepancy between the NPG and diplotype at MS-P1 ( Figure 4A ; Table S6A ); however, recombinations between markers around the MS-P1 region were not observed (data not shown). ‘Nishinokaori’ and ‘Kara’ exhibited an increased NPG, although ‘Nishinokaori’ harbored HT1/HT3 and ‘Kara’ harbored HT1/HT4, both of which are partial male sterile diplotypes ( Figure 4A ; Table S6A ). Further, two independent Rf loci corresponding to one CMS-associated gene exist in sorghum, radish, wheat, and rice (Yasumoto et al., 2009; Jordan et al., 2010; Castillo et al., 2014; Tang et al., 2017). Minor QTLs associated with the NPG (MS-P2 and MS-P3) were detected in our previous study (Goto et al., 2018); thus, these loci may distort the NPG in ‘Nishinokaori’ and ‘Kara’.

The NPG varied from 1 to 1,800 in the Satsuma mandarin strains from 2017 and 2018 ( Table S6A ). Therefore, they exhibited partial male sterility according to the CMS genetic model; however, the pollen grains were not released from the inside of anthers. Accordingly, it was proposed here that Satsuma mandarin has male sterility contributing to seedless fruits (Goto et al., 2018). Considering these previous finding, along with those of the present study, an NPG of < 1,800 pollen grains per anther was proposed as a criterion of male sterility for seedless citrus breeding. Here, an NPG of < 1,800 was detected in breeding germplasms harboring HT1/HT1, HT1/HT2, HT1/HT3, HT1/HT4, HT1/HT5 ( Figure 4A ), and HT2/HT2 ( Figure 4B ); thus, these breeding germplasms maintain male sterility for seedless citrus breeding. Male sterile individuals can be selected from F₁ populations through screening for individuals harboring HT1/HT1, HT1/HT2, HT1/HT3, HT1/HT4, HT1/HT5, and HT2/HT2, which are diplotypes genotyped using 00918-1, 00918-2, and 00918-3 at the MS-P1 region.

HT1 was the most frequent haplotype, followed by HT6 and HT2 among the 85 varieties and selected strains ( Table 2 ). This may have been because Japanese breeders primarily used Satsuma mandarin and ‘Kiyomi’ as seed parents ( Figure S2 ). Additionally, the selection of seedless lines could have resulted in a higher frequency of HT1 and HT2, as both contribute to seedlessness ( Table 2 ). HT6 is an Rf haplotype; however, it was the second most frequent haplotype. ‘Encore’, Ponkan, and their progenies, which have the HT6 haplotype, have been frequently used as parents for citrus breeding programs in Japan (Imai et al., 2017), likely indicating why HT6 was the second most frequent haplotype observed here.

In hybridization-based citrus breeding programs, only the screening of large numbers of seedlings can guarantee the identification of new varieties with good quality traits; however, citrus has a long juvenile phase and requires extended periods for evaluating fruit quality. Therefore, the selection of novel variety candidates is both time- and cost-prohibitive (Raveh et al., 2020); although this limitation can be overcome by applying marker-assisted selection at the seedling stage (Roose, 2007). Accordingly, the findings of the present study can aid in marker-assisted selection for seedless citrus breeding through screening male sterile individuals in F₁ populations.

This study has a limitation that should be addressed in future research. Rf-MS-P1 candidates were predicted to localize in the mitochondria by only bioinformatic programs. In addition to further fine mapping of MS-P1 region, cellular localization analysis of Rf-MS-P1candidates could help identify the actual Rf-MS-P1. In addition, the results here suggest that Rf loci are also present, in addition to Rf-MS-P1. Accordingly, identifying all Rf loci, and identifying the Rf genes at the Rf loci will enable a deeper understanding of the CMS mechanism in citrus