C. clementina Ciclev10004681m.g and Ciclev10027731m.g were selected as candidate genes considering the presence of SNPs in the genomic sequence of seedlessness mutants according to our genome databases (data not shown). Ciclev10004681m (XP_006421212.1) protein sequence showed the presence of conserved motifs characteristic of the Glycosyltransferase family 92 domain (IPR008166) (Supplementary Figure S1A) and therefore was named as Glycosyltransferase 1 (CcGLT1). Ciclev10027731m (XP_024035447.1) protein sequence showed the presence of conserved motifs characteristic of the family to which it belongs, such as the RNA recognition motif domain (IPR000504) and Zinc finger, RING-type (IPR001841) (Supplementary Figure S1B) and therefore named as RNA binding protein family 1 (CcRBP1).

Both candidate genes were also selected for their prevailing expression in male and female tissues according to Citrus Annotation Project (CAP, http://citrus.hzau.edu.cn/orange/) database (Figures 1A, B). In silico gene expression analysis showed that CcGLT1 transcription is much higher in ovules than in other tissues (Figure 1A). Similarly, high accumulation in ovules was found for CcRBP1 transcript (Figure 1B) together with a much higher accumulation in leaves (Figure 1B).

To explore the diversification history of the candidate proteins in relation to other dicotyledonous plant orthologs, a phylogenetic tree was constructed for each of them. Results revealed four main groups of GLT1 (Figure 1C): one group consists of GLT1 proteins from sweet orange and trifoliate orange, both parents of citrange, showing a strong closeness with CcGLT1 and a protein of Manihot esculenta; the second group contains proteins from two members of the Salicaceae family (Populus trichocarpa and Salix purpurea) and Carica papaya; the third group contains proteins from two members of the Rosaceae family (Malus domestica and Prunus persica), two members of Malvaceae family (Gossypium raimondii and G. hirsutum) and Theobroma cacao; the fourth group, showing a greater distance from the proteins of citrus (group I, Figure 1C), contains proteins from Arabidopsis thaliana, Vitis vinifera and Olea europaea. Phylogeny analysis for CcRBP1 revealed five main groups (Figure 1D): the first group includes proteins from Clementine, trifoliate and sweet orange, all showing a great closeness, and a protein from Anacardium occidentale; the second group contains proteins from T. cacao, M. esculenta, P. trichocarpa and two members of genus Gossypium (G. hirsutum and G. raimondii); the third group includes two proteins of Vitis vinifera and Castanea dentata; the fourth group contains the proteins of Rosaceae P. persica and M. domestica; the fifth group, the farthest from citrus (group I, Figure 1D), contains proteins of A. thaliana (RBP1.1 and RBP1.2).

To further explore the envisaged transcription of CcGLT1 and CcRBP1 in reproductive tissues, the 2-kb upstream sequence of both genes was analyzed in silico and compared to those of other plant orthologs (Supplementary Figure S2A). Results revealed the presence of cis-acting elements along 700-bp upstream of the coding sequence of both genes that are reported to be involved in pollen and seed expression (Supplementary Figure S2B).

Thus, the in silico expression data and the presence of cis-acting regulatory sequences specific for reproductive tissues encouraged us to further characterize CcGLT1 and CcRBP1 transcription patterns. Expression of candidate genes CcGLT1 and CcRBP1 was assayed by qRT-PCR analysis in different tissues including leaves at two developmental stages (young and mature), ovaries and anthers from flower buds, flowers at pre-anthesis and anthesis and young fruits from Carrizo citrange (Figure 2). Results showed that in female tissues both CcGLT1 and CcRBP1 transcripts accumulated in the ovary during the flower bud stage and that their expression was higher than in ovaries from later developmental stages (three and two-fold increased respect ovary during pre-anthesis and at anthesis, respectively for CcGLT1; two-fold increased for CcRBP1) and also higher than in petals (five and three-fold higher for CcGLT1 and CcRBP1, respectively) or leaves (six and three-fold higher for CcGLT1 and CcRBP1, respectively) (Figure 2B). In male tissues, both CcGLT1 and CcRBP1 transcripts accumulated at higher levels in anthers during the anthesis stage than in anthers from previous developmental stages (two and three-fold increase respect to anthers during flower bud and pre anthesis stages, respectively for CcGLT1 and four and three-fold increase compared to anthers during flower bud and pre anthesis stages, respectively for CcRBP1) (Figure 2B).

Characterization of selected candidate genes was attempted using a biotechnological approach that combines their silencing by RNAi and the concurrent overexpression of CsFT gene to induce early flowering and fruiting phenotype and thus speed up the characterization of transgenic fruit (Pons et al., 2014). Three plasmids containing both an CsFT overexpression cassette, the nptII selectable marker cassette and an intron-hairpin (ihp) cassette to silence either CcGLT1, CcRBP1 or both genes, were generated (Figure 3A) and used to transform epicotyl segments of Carrizo citrange via Agrobacterium tumefaciens-mediated transformation. The plasmid pROKII-CsFT (Figure 3A) previously generated (Pons et al., 2014) was used to transform control (CN) plants. Kanamycin-resistant regenerants obtained after transformation were screened by PCR to verify the integration of the different T-DNAs using primers specific to the CsFT and nptII transgenes as well as to each ihp cassette (Supplementary Table S1). PCR-positive shoots, designated as CN (control plants transformed with pROKII-CsFT), GLT (plants transformed with pROKII-CsFT-CcGLTi), RBP (plants transformed with pROKII-CsFT-CcRBPi) and GLT-RBP (plants transformed with pROKII-CsFT-CcGLTi-CcRBPi) lines, were grafted onto vigorous citrus rootstocks in a greenhouse. Eight lines for each construct were obtained. CsFT overexpression caused a clearly visible early flowering phenotype, with most transgenic lines starting to flower during in vitro regeneration and again when grafted on vigorous rootstocks in the greenhouse, about 3 months after the transformation experiments, mainly exhibiting a terminal flower bud indicating that the branch had transitioned to flowering. Six months after grafting in the greenhouse, transgenic plants reached full flowering (Figure 3B, CN and ihp lines) and about six months later they reached the full fruiting stage (Figure 3C). Early flowering and fruiting phenotype confirmed the overexpression of the CsFT cassette in the transgenic lines.

Differences were observed in plant size and architecture among transgenic lines when compared to the non-transformed control plants (WT) because of the effect of CsFT overexpression in plant architecture. No significant difference was found between CN and GLT plants in relation to plant architecture, evaluated by measuring the branch length/number of nodes per plant. Conversely, RBP and GLT-RBP plants showed reduced sizes when compared to the other groups, CN and GLT (Supplementary Figure S3), indicating that downregulation of CcRBP1 was likely affecting plant development, further compromising growth. RBP plants with a strong-reduced size were removed from the analysis and were evaluated plants for the CcGLT1 gene knocking-down.

For the characterization of CcGLT RNAi and its effect on reproductive development, three plants were selected among GLT and GLT-RBP lines, based on their similar architectural phenotype, meaning that they were showing early flowering and fruiting, but their vegetative development was not compromised excessively. The selected transgenic plants included: GLT.41, GLT.73 and GLT.115 for GLT; GLT-RBP.15, GLT-RBP.71 and GLT-RBP.291 for GLT-RBP.

Knocking-down of the GLT1 gene in GLT and GLT-RBP lines was evaluated by measuring endogenous GLT1 transcript abundance in ovaries from flower buds (Figure 4A) and anthers from flowers at anthesis (Figure 4B) compared to that of CN plants. About ten-fold decrease of the GLT1 gene was observed in GLT and GLT-RBP lines in the ovary at the flower bud stage (Figure 4A) with a significant difference compared to the control CN line. In anthers from flowers at anthesis, downregulation of the GLT1 was significant in all GLT-RBP and GLT lines analyzed, with reductions up to twelve-fold (GLT-RBP.291) (Figure 4B).

Transgenic plants were characterized by microscopic and histological observations to evaluate male and female performance. Male performance was evaluated by monitoring pollen viability, pollen grain germination and pollen tube growth by in vitro culture tests. Results from FDA assays showed a significant reduction of pollen viability in transgenic lines GLT and GLT-RBP compared to the CN control, while no significant differences were observed in RBP lines, suggesting that RBP1 downregulation in RBP and GLT-RBP lines was not enough to compromise pollen survival (Figure 5, Table 1).

Similar results were obtained for pollen germination rate in vitro, with a significant reduction of the percentage of pollen germination specifically in transgenic lines GLT.41 and GLT.73 and all GLT-RBP lines compared to the CN (Table 1), consistent with pollen viability results. However, no difference was found in the pollen tube length from transgenic lines, independently of their genetic background (Table 1).

Given the reduced pollen viability and germinability detected in some of the ihp lines, cross-tests were conducted on flowers from Clementine mandarin using pollen from CN, one GLT line (GLT.115) and one GLT-RBP line (GLT-RBP.15). Although seed number in Clementine fruit was lower when pollen from GLT and GLT-RBP lines was used, differences were not significant compared to CN pollen, as the number of seeds formed per fruit was high in all cases, varying between 19 and 23 (Supplementary Figure S4). These results indicated that pollen viability was not so compromised to preclude fecundation of compatible ovaries from regular citrus flowers and consequent seed development.

Female performance was assessed by measuring ovule degeneration from self-pollinated pistils stained with aniline blue. Results showed a strong accumulation of callose in the chalazal region from GLT and GLT-RBP transgenic lines, whose fluorescence was observed by aniline blue staining, indicating ovule degeneration with a significant difference compared with CN (Figures 6A, B).

Seed development was assessed by counting the number of aborted ovules/seeds and regular seeds per fruit. Results showed a significantly reduced number of seeds per fruit in all the GLT and GLT-RBP transgenic lines compared to the CN control (Figures 7A, B).

Moreover, some of the ovules/seeds had aborted in GLT and GLT-RBP lines while most of them were fully formed in CN fruits (Supplementary Table S2), indicating that seed viability was partially compromised, likely through the action of GLT1 downregulation, though not completely as some fertilized ovules remained viable to produce regular seeds.

Genome databases from proprietary seedless and seeded varieties [irradiated mandarin (Citrus reticulata Blanco) versus non-irradiated counterparts] were used to select candidate genes affected by SNPs in their genomic sequence. Also, in silico gene expression analysis was performed to evaluate pollen and/or ovule/seed tissue distribution of candidate genes according to a citrus transcriptome database (Citrus Annotation Project: http://citrus.hzau.edu.cn/). Sequence similarity searches to candidate genes CcGLT1 (Gene ID: 18032927) and CcRBP1 (Gene ID: 18034564) were conducted using the BLASTp program on NCBI database. To explore the phylogenic relationships protein sequences of the candidate genes and orthologs were obtained from Phytozome v13 database (https://phytozome-next.jgi.doe.gov/). Protein sequences (XP_006421212.1 and XP_024035447.1) were submitted to the InterPro database (http://www.ebi.ac.uk/interpro/) to determine their belonging to a characterized family and/or the presence of functionally important domains and sites. The amino acid sequences were aligned using muscle alignment tool and phylogenetic trees were created through MEGA 11 program (Tamura et al., 2021) with the neighbor-joining (NJ) method and bootstrap was set at 1000 replications. The upstream regulatory regions (2.0-kb) of the candidate genes were obtained from Phytozome v13 database and characterized to determine cis-acting elements associated with their expression in pollen and/or ovule/seed organs using the New PLACE database (https://www.dna.affrc.go.jp/PLACE/?action=newplace). The results were submitted to the online Gene Structure Display Serve (GSDS, http://gsds.cbi.pku.edu.cn/) for visualization.

Expression of candidate genes CcGLT1 and CcRBP1 was assayed by quantitative reverse transcription PCR (qRT-PCR) analysis in different tissues from Carrizo citrange (Citrus sinensis L. Osb. x Poncirus trifoliata L. Raf.) plants including leaves at two developmental stages, ovaries and anthers from flower buds, flowers at pre-anthesis and anthesis and 2-week-old (after anthesis) 5-mm long young fruits (Figure 1A). Plant materials were collected from trees held in the germplasm collection at IVIA, in Moncada, Valencia, Spain. Collected samples were immediately frozen in liquid nitrogen, ground to a fine powder and stored at -80°C until further use. Total RNA was extracted from samples and treated with rDNase RNase-free using Nucleo Spin^® RNA Plant kit (Macherey-Nagel, Germany). Total RNA was quantified using a NanoDrop^®ND-1000 (NanoDrop products, Thermo Fisher Scientific, United States) spectrophotometer. First-strand cDNA was synthesized from 1 µg of each DNase-treated RNA using oligo(dT)18 and SuperScriptTM II Reverse Transcriptase (Invitrogen, Thermo Fisher Scientific, United States) according to the manufacturer’s instructions. qRT-PCR analysis was run on QuantStudio™ 3 Real-Time PCR System (Thermo Fisher Scientific, United States) on 500 ng of total cDNA adding 6 µL of SYBR Green PCR Master Mix (Applied Biosystems, USA) and 0.3 µM of gene-specific primers in a total volume of 12 µL. The primer pairs used for GLT1 and RBP1 expression were designed based on citrus coding sequences from Citrus clementina Hort. ex Tan. database available in Phytozome v13 database. Primer sequences are detailed in Supplementary Table S1. Citrus LCY1 (Alquézar et al., 2009) and UPL7 (Mafra et al., 2012) genes were used as housekeeping reference genes. The reactions were subjected to temperature cycling as follows: 50°C (20 s), 95°C (10 min), then 40 repeats of 95°C (15 s) and 60°C (40 s), followed by 95°C (15 s), 60°C (1 min) and 95°C (15 s). The amplification efficiencies of each primer pair and the dynamic range of all the genes analyzed were determined by monitoring the variation of ΔCT by using a 10-fold dilution series of a mix of cDNA samples from different tissues as a standard curve. To demonstrate the expression stability of the reference genes LCY1 and UPL7 under our experimental conditions, the algorithm geNorm was used (https://genorm.cmgg.be) (Vandesompele et al., 2002). The relative expression level of the target genes normalized to the expression of the housekeeping genes (LCY1 and UPL7) was calculated following the mathematical model described by Livak and Schmittgen (2001). The values reported are the mean ± SD of at least three independent assays. Statistical analyses were performed using ANOVA (p < 0.01).

Using GoldenBraid cloning system v3.0, three RNAi vectors were generated to (i) silence the expression of the CcGLT1 gene, (ii) silence the expression of the CcRBP1 gene and (iii) to silence the expression of both CcGLT1 and CcRBP1 genes. A schematic representation of the T-DNA region of the constructs used in this work is shown in Figure 3A. The 404-bp fragment corresponding to the sequence used in the ihp CcGLTi (nucleotide positions 83–486 of the 6084-bp CcGLT1 gene sequence) was PCR-amplified from a gBlocks Gene Fragment (IDT) using primers B433 and B434 for sense orientation and primers B437 and B438 for antisense orientation. The 415-bp fragment corresponding to the sequence used in the ihp CcRBPi (nucleotide positions 3914–4328 of the 8391-bp CcRBP1 gene sequence) was PCR-amplified from a gBlocks Gene Fragment (IDT) using primers B409 and B427 for sense orientation and primers B439 and B440 for antisense orientation. The 819-bp fragment corresponding to the sequence used in the ihp CcGLTi-CcRBPi (nucleotide positions 83–486 and 3914–4328 of the 6084-bp CcGLT1 and 8391-bp CcRBP1 genes sequences, respectively) was PCR-amplified from a gBlocks Gene Fragment (IDT) using primers B433 and B427 for sense orientation and primers B439 and B438 for antisense orientation (Supplementary Table S1). PCR products were gel-purified using an E.Z.N.A. Cycle pure kit (Omega, Norcross, United States). The pUPD2 plasmid and 40 ng of each purified PCR amplicon were used to perform a BsmBI GoldenBraid (GB) reaction as described by Sarrion-Perdigones et al. (2013). pUPD2 sense and pUPD2 antisense parts for each ihp cassette were subjected to further GB reaction with BsaI and assembled into the pDGB3-α2 vector separated by an intron (GB01281) and under the control of the CaMV 35S promoter (GB0552) and the nos (nopaline synthase) terminator (GB0037). Sequences of GB-Parts are accessible at GB cloning website (https://goldenbraidpro.com/) using the GB database ID. The ihp cassettes were excised from pDGB3-α2 by digestion with HindIII (New England Biolabs) and ligated in the binary vector pROKII-CsFT (Pons et al., 2014) digested with the same restriction enzyme and dephosphorylated, to obtain the final pROKII-CsFT-CcGLTi, pROKII-CsFT-CcRBPi and pROKII-CsFT-CcGLTi-CcRBPi RNAi plasmids (Figure 3A). The empty plasmid pROKII-CsFT was used to transform control plants (Figure 3A). After confirmation by plasmid restriction analysis and by sequencing, each vector was transferred to A. tumefaciens strain EHA105 by thermal shock. T4 DNA ligase was purchased from Promega, BsaI from New England Biolabs, and BsmBI from Fermentas. Plasmid extractions were made by using The E.Z.N.A. Plasmid Mini Kit I (Omega Bio-Tek). Escherichia coli XL1-Blue strain was used for gene cloning.

Epicotyl explants from five-week-old Carrizo citrange in vitro-grown seedlings were used for transformation as described in Peña et al. (2004). Briefly, epicotyl explants about 1 cm in length were incubated for 5 min in the bacterial suspension. After removing bacterial debris and drying, explants were transferred to a solid co-cultivation medium [4.4 g/L MS salts (Murashige and Skoog, 1962), 3% (w/v) sucrose, 100 mg/L myo-inositol, 1mg/L nicotinic acid, 1 mg/L pyridoxine hydrochloride, 0.2 mg/L thiamine hydrochloride, 2 mg/L indole-3-acetic acid (IAA), 1 mg/L 2-isopentenyl-adenine (2-ip), 2 mg/L 2,4-dichlorophenoxyacetic acid (2,4-D) and 8 g/L agar, pH 5.7] and maintained for 3 days in the semi-dark at 25°C. Then, the explants were transferred to solid selection medium [4.4 g/L MS salts, 3% (w/v) sucrose, 100 mg/L myo-inositol, 1mg/L nicotinic acid, 1 mg/L pyridoxine hydrochloride, 0.2 mg/L thiamine hydrochloride, 3 mg/L BAP, 8 g/L agar, pH 5.7] supplemented with 100 μg/L kanamycin for nptII selection and 250 μg/L of vancomycin and 500 μg/L of cefotaxime to control bacterial growth. The plates were maintained in the dark at 25°C for two weeks and then transferred to 16 h photoperiod at 25°C. After four-five weeks, regenerated shoots were shoot-tip grafted in vitro onto Carrizo citrange seedlings as described in Peña et al. (2008). Seeds used for seedling production came from mother seed trees located at the IVIA germplasm bank in Moncada, Valencia, Spain.

For molecular analysis, genomic DNA from regenerated plantlets was extracted using the CTAB method as described by McGarvey and Kaper (1991). The presence of the T-DNA was checked by PCR using pairs of specific primers for the different RNAi cassettes (Supplementary Table S1). For CsFT, we amplified the region encompassing the end of the 35S promoter (35Sfinal-F) and the entire CsFT transgene (FTcs2) to avoid non-specific amplification of the endogenous FT gene. Similarly, for CcGLT1 and CcRBP1 we used primer designed on the 35S promoter (35Sfinal-F) and the intron (B406 for pROKII-CsFT-CcGLTi and B427 for pROKII-CsFT-CcRBPi, and pROKII-CsFT-CcGLR-CcRBPi). Approximately 2 months after grafting in vitro, the PCR-positive plantlets were grafted in the greenhouse onto rough lemon (C. jambhiri Lush.) rootstocks. PCR was performed using FIREPol^® DNA Polymerase following the manufacturer’s instructions. Reactions for CsFT and nptII amplification were carried out under the following conditions: 95°C for 5 min, 35 cycles of 95°C for 30 s, 57°C for 30 s and 72°C for 45 sec, followed by 72°C for 10 min. Reactions for ihp constructs amplification were carried out using conditions of 95°C for 5 min, 40 cycles of 95°C for 30 s, 60°C for 30 s and 72°C for 1.5 m, followed by 72°C for 10 min. PCR products were detected by electrophoresis on 1% agarose gels. Silencing of the CcGLT1 and CcRBP1 genes in male and female reproductive tissues of the transgenic plants was assayed by qRT-PCR analysis as detailed above. Primer sequences used are indicated in Supplementary Table S1.

Histological observations were performed to evaluate pollen grain germination, pollen tube growth, and ovule degeneration by in vitro culture tests. A minimum of ten flowers per transgenic line was collected. Anthers were removed from the flowers at pre-anthesis and were dried in Petri dishes over silica gel at room temperature. Then, pollen from dehiscent anthers was used for in vitro analysis and histological observations. Pollen germination was evaluated as described by Bennici et al. (2019). Pollen viability was assessed by fluorescein diacetate (FDA) staining as described by Heslop-Harrison and Heslop-Harrison (1970). Around 300 pollen grains for each flower were manually tagged on random images. Pollen with bright green fluorescence was classified as viable, while pollen with diminished fluorescence was labeled as dead.

Cross-tests were conducted on flowers from Clementine mandarin to test pollen performance. A batch of 30 flowers randomly selected from the three Clementine plants grown in a greenhouse, located at the IBMCP, in Valencia, Spain, were cross-pollinated one day before anthesis. Flowers from Clementine plants (ten flowers per line, one Clementine plant per line) were emasculated, hand-pollinated with a small paint brush with pollen from CN, GLT.115 and GLT-RBP.15 plants and bagged in cotton tissue. Fruits form cross-pollinated flowers were collected at maturity and evaluated for seed content.

Ovule degeneration was evaluated from self-pollinated transgenic plants. Ten flowers per transgenic line was collected. Hence, anthers were removed from the flowers at pre-anthesis and were dried in Petri dishes over silica gel at room temperature. Then, dehiscent anthers were used to pollinate the emasculated flowers at anthesis of the transgenic plants. The pollinated flowers were bagged to avoid any undesired cross-pollination. The pistils from the self-pollinated flowers were collected after 7 days and fixed in FAA solution (formalin, glacial acetic acid, 70% ethanol, 1:1:18, v/v) (Johansen, 1940) and stored at 4°C until the histological and microscopic observations. The pistils fixed in FAA were washed three times in distilled water and sliced into cross sections as described by Montalt et al. (2019). Then, slices were stained with 0.1% aniline blue in 0.1 N K₃PO₄. Around 200 ovules from ten flowers per transgenic line were evaluated. Microscopic observations were performed under a Leica DM5000B microscope (Leica Wetzlar, Germany) by bright-field and fluorescence microscopy for both pollen germination and pollen tube growth and both pollen viability and ovule degeneration evaluation, respectively. Pollen viability was assessed by FDA staining using excitation and emission wavelengths of 494 and 530 nm, respectively. Callose depositions for ovule degeneration evaluation were evaluated using the excitation wavelength of 370 nm and the emission maximum of 509 nm. Images were acquired with a Leica DFC 550 digital camera and analyzed using ImageJ software (https://imagej.nih.gov/). Seed development was assayed by counting the number of seeds in 8 to 21 fruits per plant at the full-colored stage.