Transcriptome and Gene Editing Analyses Reveal MOF1a Defect Alters the Expression of Genes Associated with Tapetum Development and Chromosome Behavior at Meiosis Stage Resulting in Low Pollen Fertility of Tetraploid Rice

Autotetraploid rice is a useful rice germplasm for polyploid rice breeding. However, low fertility limits its commercial production. A neo-tetraploid rice with high fertility was developed from the progenies of crossing between autotetraploid lines by our research group. Our previous study showed that a myeloblastosis (MYB) transcription factor, MOF1, might be associated with the pollen development in tetraploid rice. However, little information is available about its role in pollen development in tetraploid rice. Here, we identified a new haplotype of MOF1 from neo-tetraploid rice and marked it as MOF1a. Transcriptome and qRT-PCR analysis demonstrated that MOF1a highly expressed in anthers, and displayed differential expression in neo-tetraploid rice compared to tetraploid rice line with low pollen fertility. The mutant (mof1a) of MOF1a, which was generated by CRISPR/Cas9 knockout in neo-tetraploid rice, showed low pollen fertility, and also exhibited abnormal tapetum and middle layer development, and defective chromosome behaviors during meiosis. A total of 13 tapetal related genes were found to be up-regulated in meiotic anthers of MOF1a compared with wild type plants by RNA-seq analysis, including CYP703A3, PTC1, and OsABCG26, which had been demonstrated to affect tapetal development. Moreover, 335 meiosis-related genes displayed differential expression patterns at same stage, including nine important meiosis-related genes, such as metallothionein OsMT1a. These results demonstrated that MOF1a plays an important role in pollen development and provides a foundation for understanding the molecular mechanism underlying MOF1a in reproduction of tetraploid rice.


Introduction
Rice (Oryza sativa L.) is one of the most important food crops around the world. Rice cultivars are constantly improved to ensure food security, which largely depends on abundant germplasm To understand the difference in gene expression patterns between tetraploid lines with different fertility levels, RNA-seq was employed to analyze global genes expression in the anthers during meiosis using neo-tetraploid rice, H1 with 95.15% pollen fertility, and two recombinant inbred lines, one with high pollen fertility (HF, 85.58%) and another with low pollen fertility (LF, 54.82%) ( Table 1). H1, HF, and LF indicate a neo-tetraploid line, Huaduo1, a high pollen fertility recombinant inbred line, and a low pollen fertility recombinant inbred line, respectively. Least significant difference (LSD) was used in the multiple comparison for each trait. Different letters between two samples indicated significant differences (p value < 0.01).
More than 475 million clean reads were detected from each sample. The average percentage of mapped reads to reference genome was 94.54% (Table S1). There was a high degree of correlation between three biological replicates of RNA-seq data in HF, LF, and H1 with a Pearson correlation coefficient more than 0.888 ( Figure S1).

High Expression of MOF1 Was Detected in Anther at the Meiosis Stage
Annotated by Rice Genome Annotation Project Database [37], MOF1 is comprised of 3887 base pair (bp) that included six exons and five introns. The coding sequence of MOF1 is 921 bp, and it encodes proteins with 306 amino acids, including an MYB DNA-banding domain. Resequencing and sanger sequencing showed a total of 20 single nucleotide polymorphisms (SNPs) or loci with insertion or deletion (InDels) were detected in MOF1 of H1 compared to reference genome (Nipponbare, Oryza sativa L. ssp. japonica), including two important SNPs (1078th thymine and 3463th guanine) in CDS, four loci in the 5′ UTR region, and 14 loci in introns (Table S6). Two SNPs in CDS was further verified by cloning MOF1 coding sequence from RNA of H1 anthers. The splicing of MOF1 mRNA was consistent with the reference genome. This result indicates that MOF1 in H1 encodes a protein with two different amino acids (80th Valine and 222th Serine) compared to Nipponbare (Glycine and Asparagine, respectively). This haplotype of MOF1 in H1 was designated as MOF1a in this study. The quantity of common differentially expressed genes (coDEGs) between HF and H1. (c) Venn analysis of coDEGs, DEGs in 02428-4x relative to 02428-2x [10] and DEGs preferentially expressed during meiosis [33]. Eight genes (red) overlapped among the three group DEGs.

High Expression of MOF1 Was Detected in Anther at the Meiosis Stage
Annotated by Rice Genome Annotation Project Database [37], MOF1 is comprised of 3887 base pair (bp) that included six exons and five introns. The coding sequence of MOF1 is 921 bp, and it encodes proteins with 306 amino acids, including an MYB DNA-banding domain. Resequencing and sanger sequencing showed a total of 20 single nucleotide polymorphisms (SNPs) or loci with insertion or deletion (InDels) were detected in MOF1 of H1 compared to reference genome (Nipponbare, Oryza sativa L. ssp. japonica), including two important SNPs (1078th thymine and 3463th guanine) in CDS, four loci in the 5 UTR region, and 14 loci in introns (Table S6). Two SNPs in CDS was further verified by cloning MOF1 coding sequence from RNA of H1 anthers. The splicing of MOF1 mRNA was consistent with the reference genome. This result indicates that MOF1 in H1 encodes a protein with two different amino acids (80th Valine and 222th Serine) compared to Nipponbare (Glycine and Asparagine, respectively). This haplotype of MOF1 in H1 was designated as MOF1a in this study.
In order to make sure the spatial and temporal expression of MOF1a, three internet tools (eFP Browser [38], RiceXPro [39], and GEO Profiles of NCBI) and a total of 84 RNA-seq datasets were analyzed and MOF1 displayed high expression in anthers at meiosis stage ( Figure S2, Figure S3). Moreover, qRT-PCR demonstrated high transcriptional level of MOF1a during anther development in H1, which was consistent with the RNA-seq data and three internet tools. Meanwhile, qRT-PCR results also displayed that the expression level of MOF1a in anthers continuously decreased from S7 to S11 stage (Figure 2c).
Again, the MOF1a promoter-β-glucuronidase (GUS) vector (pMOF1a::GUS) was constructed and transformed into Nipponbare calli, in which the 1619 bp sequence ahead of initiation codon of MOF1a cloned from H1 plant was used as MOF1a promoter. In agreement with qRT-PCR results, strong GUS signals were detected in anthers and continuously decreased during anther development (spikelet length from 3.5 mm to 6.5 mm). During these stages, GUS signals were also detected in stamen, lemma, and palea, but not in lodicules, pistil, sterile lemma, and rudimentary glume. The GUS signal became almost undetectable in spikelets with 7.0 mm length (Figure 2d). Taken together, MOF1a was highly expressed in stamen and its expression level decreased during stamen development.
Browser [38], RiceXPro [39], and GEO Profiles of NCBI) and a total of 84 RNA-seq datasets were analyzed and MOF1 displayed high expression in anthers at meiosis stage ( Figure S2, Figure S3). Moreover, qRT-PCR demonstrated high transcriptional level of MOF1a during anther development in H1, which was consistent with the RNA-seq data and three internet tools. Meanwhile, qRT-PCR results also displayed that the expression level of MOF1a in anthers continuously decreased from S7 to S11 stage ( Figure 2c). Again, the MOF1a promoter-β-glucuronidase (GUS) vector (pMOF1a::GUS) was constructed and transformed into Nipponbare calli, in which the 1619 bp sequence ahead of initiation codon of MOF1a cloned from H1 plant was used as MOF1a promoter. In agreement with qRT-PCR results, strong GUS signals were detected in anthers and continuously decreased during anther development (spikelet length from 3.5 mm to 6.5 mm). During these stages, GUS signals were also detected in stamen, lemma, and palea, but not in lodicules, pistil, sterile lemma, and rudimentary glume. The GUS signal became almost undetectable in spikelets with 7.0 mm length (Figure 2d).

The Mutants with Low Pollen Fertility Were Obtained by CRISPR/Cas9 in Neo-Tetraploid Rice
To understand the reproductive role of MOF1a, this gene was edited in the neo-tetraploid rice, H1, using the CRISPR/Cas9 editing system. A total of 34 T 0 transgenic plants were obtained, and three homozygous mutants (mof1a) with frameshift mutation in the targeted loci, including mof1a-1, mof1a-2, and mof1a-3, were successfully selected from T 1 . Both target loci were inserted 1 bp 'A' (adenine) in mof1a-1; 1 bp 'A' inserted into target1 of mof1a-2; 1 bp 'T' (thymine) inserted into target1 of mof1a-3; no polymorphism was detected in target2 of mof1a-2 and mof1a-3 ( Figure 2a). While predicting the protein sequence of MOF1a in mof1a according to their new genotypes, the MOF1a proteins in mof1a-1, mof1a-2, and mof1a-3 only keep 43 amino acids (aa), 56 aa, and 56 aa as same as the MOF1a protein sequence of wild type plants. The premature translation termination codons were detected at the 51st codon in mof1a-1 and the 103rd codon in mof1a-2 and mof1a-3 ( Figure 2b). mof1a-3 is T-DNA free. Five randomly selected T 1 mof1a plants from five T 0 lines and fifteen T 2 plants from mof1a-1, mof1a-2, or mof1a-3 were used for phenotyping (Table S8).
The plant type of mof1a in T 1 and T 2 was similar to WT ( Figure 3a). However, the seed setting of mof1a plants was only 36.60%, which was significantly lower than WT (74.54%). To reveal the reproductive defects of mof1a, we investigated the mature embryo sac, pollen fertility, and pollen viability of mof1a and WT. The mature embryo sac fertility of mof1a mutants was 84.35%, which was no different compared to WT (94.16%). The pollen fertilities of mof1a and WT plants were 71.00% and 93.00%, respectively. Further staining the intravital pollen grains of mof1a by 2, 3, 5-Triphenyl-2H-tetrazolium chloride (TTC) solution displayed that only 21.36% pollen grains could be stained with light reddish color and 3.95% showed a deep red color, while this was 45.76% and 27.03% for WT plants, respectively ( Figure 3). These results indicate that low pollen fertility was the main reason for the low seed set of mof1a, and MOF1a might be involved in pollen development.

MOF1a defect Affects Tapetum Development and Chromosome Behavior at Meiosis Stage of Tetraploid Rice
In anther transverse sections, no obvious differences were observed between WT and mof1a-1 during S7 stage. The middle layer in WT displayed complete degeneration and became thick-lined

MOF1a Defect Affects Tapetum Development and Chromosome Behavior at Meiosis Stage of Tetraploid Rice
In anther transverse sections, no obvious differences were observed between WT and mof1a-1 during S7 stage. The middle layer in WT displayed complete degeneration and became thick-lined structure at S8a stage. However, the anther walls of mof1a-1 from S8 to S10 stage exhibited delayed the degeneration of the middle layer, which kept a cellular shape. Moreover, the tapetum cells of WT formed vacuoles and their stained areas reduced to form a cupped shape. At the S10 and S11 stage, tapetum cells degraded and differentiated. By contrast, the tapetum cells of mof1a-1 did not form a special structure, as WT did, and kept a large shape without reduction of stained areas during the S8 stage. During the S10 to S11 stage, the tapetum cells were lightly stained and displayed inward hyperplasia in mof1a-1; at this stage, the round shape stained materials surrounding microspore were not found in mof1a-1, which was different from WT. At the S12 stage, pollen grains of mof1a-1 still contained unfilled areas, while WT pollen grains were filled with lipid and starch (Table 2, Figure 4a). S8a, S8b, and S10 indicate anther development stages. S8a and S8b represents the meiosis process. S10 represents the single microspore stage.
During meiosis, H1 exhibited a lower number of abnormalities, which were similar to normal diploid rice. However, great changes were found in meiosis process of mof1a-1. At diakinesis, no difference was found in quantity or distribution of quadrivalents between H1 and mof1a, but the proportion of PMCs that were univalent or trivalent, was raised from 16.55% (145 observed WT PMCs) to 47.17% (104 observed mof1a-1 PMCs). From metaphase I to telophase II, more PMCs with abnormal chromosome behavior were observed in mof1a-1, including chromosome dragging at metaphase I (87.86% of 173 PMCs) and metaphase II (77.94% of 136 PMCs); chromosome lagging at anaphase I (59.32% of 59 PMCs) and anaphase II (50.00% of 48 PMCs); micronuclei at telophase I (61.60% of 125 PMCs) and telophase II (58.88% of 107 PMCs); and asynchrony of the chromosome during meiosis II (6.53% of 291 PMCs) (Figure 4b). The frequency of PMCs with chromosome behavior abnormalities in mof1a-1 was higher than H1 at all observed stages ( Table 2). These results indicate that MOF1a affects chromosome behavior during the meiosis process in neo-tetraploid rice. During meiosis, H1 exhibited a lower number of abnormalities, which were similar to normal diploid rice. However, great changes were found in meiosis process of mof1a-1. At diakinesis, no difference was found in quantity or distribution of quadrivalents between H1 and mof1a, but the proportion of PMCs that were univalent or trivalent, was raised from 16.55% (145 observed WT PMCs) to 47.17% (104 observed mof1a-1 PMCs). From metaphase I to telophase II, more PMCs with abnormal chromosome behavior were observed in mof1a-1, including chromosome dragging at metaphase I (87.86% of 173 PMCs) and metaphase II (77.94% of 136 PMCs); chromosome lagging at anaphase I (59.32% of 59 PMCs) and anaphase II (50.00% of 48 PMCs); micronuclei at telophase I (61.60% of 125 PMCs) and telophase II (58.88% of 107 PMCs); and asynchrony of the chromosome during meiosis II (6.53% of 291 PMCs) (Figure 4b). The frequency of PMCs with chromosome behavior abnormalities in mof1a-1 was higher than H1 at all observed stages (Table 2). These results indicate that MOF1a affects chromosome behavior during the meiosis process in neo-tetraploid rice.

MOF1a Defects Alter Expression Levels of Important Genes Involved in Tapetum and Meiosis in Neo-tetraploid Rice
To reveal genes regulated by MOF1a, RNA-seq of anthers during meiosis from a single plant of mof1a and WT was executed. A high degree of Pearson correlation (>0.905) was detected among three biological replicates of each material (Figure 5a). We identified 793 DEGs in mof1a, including 676 up-and 117 down-regulated genes (Figure 5b, Table S9). GO analyses revealed that these DEGs were enriched in key processes associated with anther development, such as sporopollenin biosynthetic process (GO: 0080110), pollen exine formation (GO: 0010584), and suberin biosynthetic process (GO: 0010345). KEGG pathway analysis demonstrated that these DEGs were involved in biosynthesis of indole alkaloid, O-glycan, flavonoid, cutin, suberin, and wax (Figure 5d,e). The GO and KEGG analyses showed that MOF1a was associated with important biosynthetic pathways during meiosis and pollen exine formation.

MOF1a Defects Alter Expression Levels of Important Genes Involved in Tapetum and Meiosis in Neo-Tetraploid Rice
To reveal genes regulated by MOF1a, RNA-seq of anthers during meiosis from a single plant of mof1a and WT was executed. A high degree of Pearson correlation (>0.905) was detected among three biological replicates of each material (Figure 5a). We identified 793 DEGs in mof1a, including 676 upand 117 down-regulated genes (Figure 5b, Table S9). GO analyses revealed that these DEGs were enriched in key processes associated with anther development, such as sporopollenin biosynthetic process (GO: 0080110), pollen exine formation (GO: 0010584), and suberin biosynthetic process (GO: 0010345). KEGG pathway analysis demonstrated that these DEGs were involved in biosynthesis of indole alkaloid, O-glycan, flavonoid, cutin, suberin, and wax (Figure 5d,e). The GO and KEGG analyses showed that MOF1a was associated with important biosynthetic pathways during meiosis and pollen exine formation.

MOF1a Plays a Key Role in the Pollen Development of Tetraploid Rice
Myeloblastosis (MYB) transcription factors are characterized by a highly conserved MYB DNA-binding domain. The MYB gene family contains at least 155 members in rice [42]. At least five MYB transcription factors play important roles during reproductive process in rice, including OsGAmyb, which modulates anther development in rice by regulating gibberellin [43], OsMYB103 (OsMYB80), which regulates anther development by regulating the tapetum development and pollen wall formation [44,45], OsMYB106, of which the down-regulation causes pollen sterility and shortens plant height [46], CSA, which promotes rice pollen development and affects its sugar distribution via brassinosteroids [47,48], and OsTDF1, which is involved in tapetal development [49].
Recently, MOF1 was found to encode a nucleoprotein with a MYB domain and two typical ethylene response factor-associated amphiphilic repression (EAR) motifs, and to regulate organ identity as well as spikelet determinacy in rice [36]. However, information about reproductive roles of MOF1 are limited. In this study, we identified a new haplotype of MOF1 from neo-tetraploid rice, in which some important SNPs were detected; as a result, it was named MOF1a. Importantly, MOF1a is a key gene that regulates pollen development of neo-tetraploid rice. The loss function of MOF1a caused abnormal tapetum degradation and abnormal meiosis process. Tapetal development and meiotic process have been well-investigated by many researchers [31,50]. Tapetal cells undergo degradation triggered by programmed cell death from meiosis to microspore stage, which is crucial to microspore development and pollen wall formation [51,52]. Thus, MOF1a played important role in pollen development of neo-tetraploid rice, especially affecting tapetum development and the meiosis process.

MOF1a May Affect the Gene Regulatory Network Associated with Tapetum Development or the Meiosis Process in Tetraploid Rice
In this study, loss function of MOF1a caused defective tapetum development with remarkable changes in the expression levels of 14 tapetal related-genes during meiotic anthers. Function of these 14 genes have been well understanding during anther development. CYP703A3 and CYP704B2 regulate 7-hydroxylated lauric acid generation and ω-hydroxylation of fatty acids [53,54]. OsACOS12 involves in sporopollenin synthesis [55]. DPW participates in fatty alcohol synthesis for anther cuticle and sporopollenin biosynthesis [56]. OsAP37 and OsCP1 involves in programed cell death of tapetum [52,57]. Mutants of OsPKS2, OsTKPR1, PTC1, and OsMYB80 lead to abnormal tapetum degeneration and most microspores lacked exines in mature anthers [45,[57][58][59]. OsABCG15, OsABCG26, OsABCG3, and OsC6 respond to transport lipidic molecules from tapetal cells [60][61][62]. These DEGs were involved in the biosynthesis pathway of important materials such as sporopollenin, the degeneration of tapetum, and the transportation of biosynthesis production from tapetum to microspore, which are essential for microspore subsequent development. The transcripts of these genes began to accumulate at early meiosis stages and accumulated a relative high level during meiosis, including CYP703A3 (S8) [53], CYB704B4 (S8) [54], PTC1 (S8) [57], OsABCG26 (S8) [62], and WDA1 (S8) [63]. These results indicate that meiosis is an important stage for the initiation of tapetal biosynthesis, and suggest that MOF1a might regulate the expression levels of important tapetal related-genes to affect tapetum development. The expression pattern of tapetal genes are strictly controlled, and premature expression against their normal pattern may cause abnormal tapetum development. For example, a series of important tapetal gens abnormally up-regulated are expressed in DTD/EAT1 mutant (dtd), including WDA1, CYP703A3, OsAP37, OsGAmyb, OsCP1, PTC1, UDT1, and TDR. As a result, dtd cannot normally start the programmed cell death of tapetum [64]. Similarly, 13 tapetal genes were up-regulated in meiotic anthers of mof1a. Interestingly, the expression levels of all differentially expressed tapetal genes gradually increased from S7 to S10 stage, including CYP703A3, CYP704B2, OsGAmyb, OsC6, PTC1, OsCP1 [64], OsABCG26 [62], DPW [56], OsACOS12 [55], OsABCG3 [60], OsPKS2 [59], OsAP37 [52], and OsTKPR1 [58]. However, the expression level of MOF1a showed a peak in an early stage (S7) and gradually decreased from S7 to S10 stage in WT plants, which is opposite to the expression pattern of these tapetal genes. Recently, it was found that MOF1 protein has two EAR motifs with a strong repression effect on downstream genes [36]. Thus, MOF1a regulates the development of tapetum by repressing expression levels of tapetal biosynthesis related-genes in an early stage of anther development.
In addition, loss function of MOF1a also increased the frequency of PMCs with abnormal chromosome behaviors in seven stage of meiosis, and the most common abnormality was the straggle chromosome. In a previous study, many mutants of meiosis-related genes, such as ZIP4, OsREC8, MRE11, OsSGO1, HEI10, OsDMC1, and OsMND1, caused similar straggle chromosomes during meiosis process because of abnormal homologous synapsis or double strand break [65][66][67][68][69][70][71]. Disorder of microtube distribution is also considered a key reason of chromosome behaviors [9,72]. However, no known meiotic-related genes differentially expressed in mof1a anthers with defective meiosis, suggesting that MOF1a might regulate the meiosis process via another pathway, that most well-known meiotic genes might function upstream to regulate meiotic function of MOF1a, or that MOF1a mainly affect factors involved in the post-transcriptional regulation of these meiotic genes. As a first conjecture, scientists have performed abundant researches to identified meiosis-related genes, and we found 335 mDEGs in mof1a. In particular, there are nine mDEGs that were differentially expressed in more than three group samples with different meiotic process, suggesting that they may be associated with meiotic chromosome behaviors [4,10]. Among these nine genes, LOC_Os08g28820 belongs to the Skp1 family, which is involved in Skp1-Cul1-F-box-protein (SCF) complex formation and regulates plant male meiosis process in Arabidopsis [73] and rice [74]; LOC_Os01g04409 encodes an OsWAK receptor-like cytoplasmic kinase, which is a homolog of OsWAK91/OsDEES1 regulating female gametophyte development [75]; LOC_Os11g47809 (OsMT1a) encodes an metallothionein to regulate reactive oxygen species (ROS) levels, while another metallothionein OsMT2b controls ROS levels to effect pollen development [76]. All of these nine genes were up-regulated expressed in mof1a, suggesting that MOF1a might alter expression levels of a set of important meiotic related genes, particularly these nine genes, to regulate meiosis process in neo-tetraploid rice. Knowledge about these regulations is still limited, and more explorations are required to verify them in future studies.

Plant Material
One neo-tetraploid line with high fertility Huaduo1 (H1) and two recombinant inbred lines, i.e., HF with high fertility and LF with a low fertility, were used for RNA-seq. HF and LF were developed from the progenies of T428 (autotetraploid line with low fertility) × H1. T428 was developed from the chromosome doubling of diploid rice (Oryza sativa L. ssp. japonica), Linglun, by colchicine. H1 was also used as the receptor during CRISPR/Cas9 transgenes to generate mof1a mutants.

Evaluation of Seed Setting and Pollen Fertility
At least five plants of each material were harvested to measure their seed setting at maturity [1]. Pollen fertility was evaluated by 1% iodine-potassium iodide solution staining experiment [77]. The live pollen grains were subjected to 1% 2, 3, 5-Triphenyl-2H-tetrazolium chloride (TTC) solution at 31 • C for 1 h to evaluate pollen viability. In this case, the viable pollen grains would turn red while the devitalized one would keep primary color. Three spikelets from each plant (three plants for each material) were collected for pollen fertility or pollen viability experiments. Five pictures of each sample were taken under a microscope (Motic BA200) for the estimation of pollen fertility and viability.

RNA-seq Analysis
To determine that the anthers (pollen mother cells) are at meiosis stage, we measured the spikelet length. Our results showed that most of the pollen mother cells were at a meiosis stage when spikelet length was between 4.00 to 6.00 mm in H1, HF, and LF. Thus, the anthers were isolated from 4.00 mm to 6.00 mm spikelets. All the samples were kept in liquid nitrogen and stored at −80 • C for RNA isolation. Each sample was collected in three biological replications. The total RNA extraction, RNA-seq, and the evaluation of differentially expressed gene (DEG) were done as described by Guo et al. [1]. The DEGs were validated by qRT-PCR as described by Wu et al. [20]. The Ubiquitin gene was used as an internal control and all primers for qRT-PCR were designed by Primer Premier 5.0 and Primer-Blast software in NCBI (Table S13). The relative expression of qRT-PCR was calculated by the 2 −∆∆CT method [78].

Bioinformatics Analysis
Candidate genes were annotated according to GFF3 file of MSU7 genome [37]. The whole gene expression profile was predicted by using the Rice Expression Profile Database [39], Rice eFP expression profile analysis website [38], and GEO Profiles in NCBI. GO analysis of candidate genes was performed by using AgriGO tool [79]. Heat map diagram was drawn by TBtools [80].

Construction of CRISPR/Cas9 Vectors
We constructed the CRISPR/Cas9 vectors as described by Ma et al. [81].
For each candidate gene, two sgRNA expression cassettes (U6a and U6b promoters) were designed to target specific sites at the coding sequence to cause frameshift mutations (Figure 2a). The sgRNA sequence = target sequence + "GTTTTAGAGCTAGAAATAGCAAGTTAAAATAAGGC TAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGC." These two sgRNA expression cassettes were cloned into pYLCRISPR/Cas9Pubi-H by Golden Gate ligation. After transformation into Escherichia coli DH5-alpha, the recombinants were selected and sequenced. BioRun (Wuhan, China) was employed to complete the tissue culture of H1 and transform each vector into calli of H1 using Agrobacterium-mediated method. Transgenic seedlings were examined under natural field condition at the experimental station of South China Agriculture University, Guangzhou, China. In T 0 generation, two primers (gRT14 and U6bT17) were used as molecular markers to identify T-DNA, which were specific to two sgRNA expression cassettes (Table S13). In T 1 and T 2 generation, two primers (Hyg-F and Hyg-R) were used as molecular markers for the hygromycin resistance gene to identify T-DNA (Table S13). Sanger sequencing was used for the phenotypic analysis of transgenic lines. T1MOF1 and t2MOF2 primers were used to detect the polymorphisms of target1 loci and target2 loci of MOF1a (Table S13).

Cytological Observation
Young panicles were collected from rice plants with 0.00 to 3.00 cm between their flag leaf cushion and the second-to-last leaf cushion, fixed in Carnoy solution (alcohol: acetic acid, 3:1 v/v) for more than 24 h, and then kept in 70% alcohol. Chromosome behavior and configuration were observed and analyzed as described by Wu et al. [4]. Meiosis stages were classified and explained according to He et al. [9].
A semi-thin sectioning was employed to characterize the anther wall development in HF and LF. The isolated anthers from spikelets in 70% alcohol were sequentially dehydrated in 85%, 90%, and 95% alcohol for 30 min. The dehydrated samples were added to penetrant (Basic resin: Activator = 50 mL: 0.5 g, Leica 7022 Histeresin Embedding Kit) with the same volume of 95% alcohol and stored at 4 • C for 10 h. The samples were transferred to pure penetrant and kept for 3 h. Then, 500 µL embedding medium (penetrant: Hardener= 15:1 v/v, Leica 7022 Histeresin Embedding Kit) and one sample were placed into a mold together. The solution with sample would condense after 30 min. Cross-sections (3 µm) were cut and stained with 1% Toluidine Blue O for several seconds and washed under the water. Then, the samples were observed under a microscope (Motic BA200).
A whole-mount eosin B-staining confocal laser scanning microscopy (WE-CLSM) was used to investigate the variations in mature embryo sac fertility. The spikelets of pre-flowering were collected and fixed in a Formaldehyde-Acetic acid-Alcohol solution (70% alcohol: acetic acid: formaldehyde = 18:1:1, v/v) for at least 24 h. Then, the samples were stored in 70% alcohol at 4 • C. The isolated ovaries were hydrated consecutively in 50%, 30%, and 10% alcohol and distilled water for 30 min. After an eosin B (10 mg/L in 4% sucrose solution) staining procedure for 10 h, the samples were dehydrated sequentially in 10%, 30%, 50%, 70%, 90%, and 100% (three times) alcohol for 30 min. The dehydrated samples were added to methyl salicylate with the same volume of alcohol and kept for 1 h. Finally, the samples were stored in pure methyl salicylate for 1 h and observed under WE-CLSM.
Supplementary Materials: The following are available online at http://www.mdpi.com/1422-0067/21/20/7489/s1, Figure S1: The Pearson correlation between three biological replicates of RNA-seq in HF, LF, and H1. Figure S2: Predicted expression pattern of MOF1. Figure S3: Expression pattern analysis of MOF1 by RNA-seq data. Figure S4. Transcription level of MOF1a in WT and mof1a plants by qRT-PCR. Table S1: Average RNA-seq data quality among three biological replicates in H1, HF, and LF. Table S2: Differentially expressed genes (DEGs) between H1 and LF. Table S3. Differentially expressed genes (DEGs) between HF and LF. Table S4. Common differentially expressed genes (coDEGs) between HF and H1. Table S5: Eight candidate genes related to the anther development. Table S6: The SNPs and InDels of MOF1a in H1. Table S7: Detailed information about the expression pattern analysis of MOF1 by RNA-seq data ( Figure S5). Table S8: The seed setting and mutant information of WT and mof1a. Table S9: Differentially expressed genes (DEGs) between mof1a and WT. Table S10: Meiosis-related gene set. Table S11: Differentially expressed meiosis-related genes in mof1a. Table S12: Differentially expressed meiosis-related genes (mDEGs) overlapped in mof1a/H1 and LF/H1-HF. Table S13: Primers used in this study.