ZmFAR1 and ZmABCG26 Regulated by microRNA Are Essential for Lipid Metabolism in Maize Anther

The function and regulation of lipid metabolic genes are essential for plant male reproduction. However, expression regulation of lipid metabolic genic male sterility (GMS) genes by noncoding RNAs is largely unclear. Here, we systematically predicted the microRNA regulators of 34 maize white brown complex members in ATP-binding cassette transporter G subfamily (WBC/ABCG) genes using transcriptome analysis. Results indicate that the ZmABCG26 transcript was predicted to be targeted by zma-miR164h-5p, and their expression levels were negatively correlated in maize B73 and Oh43 genetic backgrounds based on both transcriptome data and qRT-PCR experiments. CRISPR/Cas9-induced gene mutagenesis was performed on ZmABCG26 and another lipid metabolic gene, ZmFAR1. DNA sequencing, phenotypic, and cytological observations demonstrated that both ZmABCG26 and ZmFAR1 are GMS genes in maize. Notably, ZmABCG26 proteins are localized in the endoplasmic reticulum (ER), chloroplast/plastid, and plasma membrane. Furthermore, ZmFAR1 shows catalytic activities to three CoA substrates in vitro with the activity order of C12:0-CoA > C16:0-CoA > C18:0-CoA, and its four key amino acid sites were critical to its catalytic activities. Lipidomics analysis revealed decreased cutin amounts and increased wax contents in anthers of both zmabcg26 and zmfar1 GMS mutants. A more detailed analysis exhibited differential changes in 54 monomer contents between wild type and mutants, as well as between zmabcg26 and zmfar1. These findings will promote a deeper understanding of miRNA-regulated lipid metabolic genes and the functional diversity of lipid metabolic genes, contributing to lipid biosynthesis in maize anthers. Additionally, cosegregating molecular markers for ZmABCG26 and ZmFAR1 were developed to facilitate the breeding of male sterile lines.


Introduction
Anther and pollen development are critical for male reproduction of flowering plants. Plant microspore mother cells are surrounded by an anther wall consisting of four layers, from outer to inner including epidermis, endothecium, middle layer, and tapetum [1]. Specifically, developing pollen grains are protected from external abiotic stresses and biotic attacks by the anther cuticle that is an extracellular lipidic layer covering anther's that functioned in the production of fatty alcohols in plastids [33,34], and its expression was directly regulated by AtMYB103/AtMS188/AtMYB80 [35]. The gene AtFAR3/CER4 encoded an endoplasmic reticulum (ER)-localized FAR that was specifically involved in the production of C24 to C28 very long-chain primary alcohols [36]. In rice, Defective Pollen Wall (OsDPW), the ortholog of AtMs2, was similarly regulated by OsMYB80 and controlled the primary fatty alcohol synthesis for the anther cuticle and sporopollenin [37,38]. In bread wheat (Triticum aestivum L.), TaTAA1 encoded by Triticum aestivum anther 1 was the first reported alcohol-forming fatty acyl-CoA reductase responsible for male gametophyte development [39]. In maize, ZmMs6021/ZmMs25 [40,41] encoded a fatty acyl-ACP reductase and was directly activated and regulated by ZmMYB84, which is an MYB TF ortholog of AtMYB80 and OsMYB80 [41]. Even though the major catalytic steps of FARs in lipid biosynthesis have been revealed, our knowledge on the function and regulation of FAR family members involved in plant male reproduction is still limited.
The present study focuses on the potential posttranscriptional regulation of ZmABCG26 by zma-miR164h-5p during maize anther development and evaluation of ZmABCG26 and ZmFAR1 as genic male sterility (GMS) genes using the CRISPR/Cas9 mutagenesis to generate mutants for phenotypic and cytological observations. The subcellular localization of the two GMS genes encoding proteins was investigated, and analysis of enzymatic activity for ZmFAR1was conducted in vitro. Finally, a comparison was conducted for lipidomics data from anthers of ZmABCG26 and ZmFAR1 knockout lines. These findings will facilitate a greater understanding of the miRNA-regulated lipid metabolic network in maize anthers and may promote the application of lipid metabolic GMS genes in molecular breeding and hybrid seed production in maize.

The Potential Posttranscriptional Regulation of ZmABCG26 by zma-miR164h-5p during Maize Anther Development
ZmABCG26 has been predicted as one of the targets of zma-miR164h-5p by bioinformatics analysis based on transcriptomes of maize anthers [42]. ZmABCG26 belongs to the WBC/ABCG subfamily that contains 34 members in the maize genome ( Figure 1A). To further investigate the potential posttranscriptional regulation of ZmABCG26 by zma-miR164h-5p during maize anther development, we predicted the miRNA regulators of all maize WBC/ABCG family genes with additional transcriptome datasets. By a combined analysis of RNA-Seq and miRNA-Seq data from three maize WT lines (W23, Oh43, and B73), as well as three maize GMS mutants (mac1, ocl4, and ms23) [12,42,43], 29 miRNA-ABCG gene pairs were predicted that contained 13 maize WBC/ABCG family members and a total of 2 known and 23 potential new miRNAs ( Figure 1A, Table S1). Among the 25 miRNAs, 13 have detectable expression values in the anther transcriptome dataset ( Figure 1B). zma-miR164h-5p was highly expressed in B73 and Oh43 during anther developmental stages 8 to 9. miRN2025 was also highly expressed during anther developmental stages 3 to 7 in W23 and mutant ocl4. In addition, 8 of the 13 expressed known or potential miRNAs were negatively correlated with their target genes in expression levels, including zma-miR164h-5p, which was predicted to bind to the ZmABCG26 transcript at the region corresponding to the fifth exon ( Figure 1C,D and Figure S1). In B73 anther transcriptomes, zma-miR164h-5p was increasingly expressed during anther developmental stages 9 to 10, and the expression level of ZmABCG26 was decreased during the corresponding stages (Figure 1(D1)). Similarly, the negatively correlated expression pattern between zma-miR164h-5p and ZmABCG26 was detected in the Oh43 anther transcriptomes during the same developmental stages (Figure 1(D2)). Furthermore, expression patterns between zma-miR164h-5p and ZmABCG26 were confirmed by quantitative reverse transcription PCR (qRT-PCR) assay in B73 and Oh43 anthers during a broader and finely divided developmental period that includes stages 9, 9-10, 10 and 11 ( Figure 1E). The qRT-PCR results showed the expression of zma-miR164h-5p was continually maintained at a high level in B73 or increased into a higher level in Oh43 at anther developmental stages 10 and 11 (Figure 1(E1,E2)), while the expression of ZmABCG26 was hardly detected in the two lines at stages 10 and 11. Therefore, by the expression anticorrelation result, it can be speculated that ZmABCG26 is most probably regulated by zma-miR164h-5p at the posttranscriptional level, and the miR164h-ABCG26 module may be functionally important in maize anther development. into a higher level in Oh43 at anther developmental stages 10 and 11 ( Figure 1(E1,E2)), while the expression of ZmABCG26 was hardly detected in the two lines at stages 10 and 11. Therefore, by the expression anticorrelation result, it can be speculated that ZmABCG26 is most probably regulated by zma-miR164h-5p at the posttranscriptional level, and the miR164h-ABCG26 module may be functionally important in maize anther development. transcriptome analyses of expression patterns of predicted miRNA regulators of maize WBC/ABCG family genes in anthers of WT inbred lines (W23, B73, and Oh43) and GMS mutant lines (ocl4, mac1, ms23). The names of predicted miRNA regulators having negatively correlated expression levels with their ABCG target genes (Pearson correlation test) were marked in blue; (C) the sequence alignment of zma-miR64h-5p and its target site on the ZmABCG26 transcript. *, base mismatch; (D,E) the negatively correlated expression patterns between ZmABCG26 and zma-miR164h-5p in B73 (D1,E1) and Oh43 (D2,E2) inbred lines at anther developmental stages 9 and 10 revealed by transcriptome analyses (D) and at stages 9 to 11 confirmed by qRT-PCR analyses (E). r, Pearson correlation coefficient.

CRISPR/Cas9-Induced Loss Function of ZmABCG26 Leads to Complete Male Sterility in Maize
The phylogenetic analysis predicted that ZmABCG26 (Zm00001d046537) was an orthologous gene of GMS genes OsABCG15 and AtABCG26 (Figure 2A). Using the Figure 1. Posttranscriptional regulation of ZmABCG26 by zma-miR164h-5p during maize anther development: (A) the phylogenetic tree of maize ABCG family genes and their predicted miRNA regulators during anther development; (B) transcriptome analyses of expression patterns of predicted miRNA regulators of maize WBC/ABCG family genes in anthers of WT inbred lines (W23, B73, and Oh43) and GMS mutant lines (ocl4, mac1, ms23). The names of predicted miRNA regulators having negatively correlated expression levels with their ABCG target genes (Pearson correlation test) were marked in blue; (C) the sequence alignment of zma-miR64h-5p and its target site on the ZmABCG26 transcript. *, base mismatch; (D,E) the negatively correlated expression patterns between ZmABCG26 and zma-miR164h-5p in B73 (D1,E1) and Oh43 (D2,E2) inbred lines at anther developmental stages 9 and 10 revealed by transcriptome analyses (D) and at stages 9 to 11 confirmed by qRT-PCR analyses (E). r, Pearson correlation coefficient.

CRISPR/Cas9-Induced Loss Function of ZmABCG26 Leads to Complete Male Sterility in Maize
The phylogenetic analysis predicted that ZmABCG26 (Zm00001d046537) was an orthologous gene of GMS genes OsABCG15 and AtABCG26 (Figure 2A). Using the CRISPR/Cas9 technology zmabcg26 mutants were generated to evaluate the functional importance of ZmABCG26 in maize male reproduction. Firstly, a pCas9-ZmABCG26 vector with two targets in the first exon of ZmABCG26 was constructed ( Figure 2B) and then transformed into maize. PCR products containing the target site fragment of primary T 0 transgenic plants were sequenced and three homozygous loss-of-function mutants of ZmABCG26 with 7-bp indel, 26-bp deletion, and 6-bp indel, respectively, were obtained ( Figure 2C,D). The three zmabcg26 homozygous mutants exhibited normal vegetative organ growth and female development but showed shrunken anthers that failed to produce viable pollen grain, as compared to WT plants ( Figure 2E). Scanning electron microscopy (SEM) assay indicated that the outer and inner surfaces of zmabcg26 mutant anther wall were smooth and glossy without knitting cuticle and Ubisch bodies ( Figure 2F). These results indicated that ZmABCG26 plays an essential role in maize anther and pollen development. CRISPR/Cas9 technology zmabcg26 mutants were generated to evaluate the functional importance of ZmABCG26 in maize male reproduction. Firstly, a pCas9-ZmABCG26 vector with two targets in the first exon of ZmABCG26 was constructed ( Figure 2B) and then transformed into maize. PCR products containing the target site fragment of primary T0 transgenic plants were sequenced and three homozygous loss-of-function mutants of ZmABCG26 with 7-bp indel, 26-bp deletion, and 6-bp indel, respectively, were obtained ( Figure 2C,D). The three zmabcg26 homozygous mutants exhibited normal vegetative organ growth and female development but showed shrunken anthers that failed to produce viable pollen grain, as compared to WT plants ( Figure 2E). Scanning electron microscopy (SEM) assay indicated that the outer and inner surfaces of zmabcg26 mutant anther wall were smooth and glossy without knitting cuticle and Ubisch bodies ( Figure 2F). These results indicated that ZmABCG26 plays an essential role in maize anther and pollen development.  T 0 plants did not produce seeds by self-pollination, and maize inbred line Zheng58 was used as a male parent to pollinate plants. To avoid the complications from the continuing action of the Cas9/gRNA transgene element in Zheng58 alleles, Cas9-negative F 1 plants were selected to produce F 2 seeds. To detect genotypes of the zmabcg26 mutant F 2 plants, the cosegregating molecular markers with the primer pairs covering the corresponding mutations in ZmABCG26 were developed, by which the mutant genotypes were identified using PCR amplification, together with agarose gel electrophoresis or polyacrylamide gel electrophoresis (PAGE) ( Figure 2G, Table S2).

Spatiotemporal Expression Analysis of ZmABCG26 Gene and Subcellular Localization of ZmABCG26 Protein
ZmABCG26 was expressed in anthers from stages 8b to 9-10, with the expression peak at stage 9-10, while its expression was not detected in the vegetative tissues (root, stem, and leaf) and female organs (immature cob and silk) ( Figure 3A), indicating that ZmABCG26 is an anther-specific expression GMS gene. The expression pattern of ZmABCG26 is consistent with its functions in maize anther and pollen development.
plants, the cosegregating molecular markers with the primer pairs covering the corresponding mutations in ZmABCG26 were developed, by which the mutant genotypes were identified using PCR amplification, together with agarose gel electrophoresis or polyacrylamide gel electrophoresis (PAGE) ( Figure 2G, Table S2).

Spatiotemporal Expression Analysis of ZmABCG26 Gene and Subcellular Localization of ZmABCG26 Protein
ZmABCG26 was expressed in anthers from stages 8b to 9-10, with the expression peak at stage 9-10, while its expression was not detected in the vegetative tissues (root, stem, and leaf) and female organs (immature cob and silk) ( Figure 3A), indicating that ZmABCG26 is an anther-specific expression GMS gene. The expression pattern of ZmABCG26 is consistent with its functions in maize anther and pollen development.
Transiently expressed ZmABCG26-GFP constructs in maize protoplasts and tobacco leaves, respectively, were used to investigate the subcellular localization of ZmABCG26. In maize protoplasts, the ZmABCG26-GFP fusion protein was colocalized with both autofluorescence of chlorophyll in plastids and the ER marker protein mCherry-HDEL (Figure 3(B1)). In tobacco leaf epidermis, ZmABCG26-GFP fusion protein was colocalized with both the ER marker and plasma membrane marker proteins (Figure 3(B2)). Therefore, ZmABCG26 is primarily localized in ER, plasma membrane, and chloroplasts/plastids, implying that the functions of ZmABCG26 and its orthologs may be slightly different between plant species.  Transiently expressed ZmABCG26-GFP constructs in maize protoplasts and tobacco leaves, respectively, were used to investigate the subcellular localization of ZmABCG26. In maize protoplasts, the ZmABCG26-GFP fusion protein was colocalized with both autofluorescence of chlorophyll in plastids and the ER marker protein mCherry-HDEL (Figure 3(B1)). In tobacco leaf epidermis, ZmABCG26-GFP fusion protein was colocalized with both the ER marker and plasma membrane marker proteins (Figure 3(B2)). Therefore, ZmABCG26 is primarily localized in ER, plasma membrane, and chloroplasts/plastids, implying that the functions of ZmABCG26 and its orthologs may be slightly different between plant species.

CRISPR/Cas9-Based Confirmation of ZmFAR1 Function in Controlling Maize Male Sterility
ZmFAR1 (GRMZM2G120987/Zm00001d048337) is the maize ortholog of GMS genes OsDPW and AtMs2 encoding fatty acyl reductases (FARs) [34,37]. One CRISPR/Cas9 vector of ZmFAR1 was constructed and transformed into maize ( Figure 4A). Three homozygous zmfar1 mutants (ZmFAR1-Cas9-1, -2 and -3) with 198-bp, 27-bp, and 3-bp deletions in the first exon, respectively ( Figure 4B,C), were produced for the evaluation. zmfar1 mutant F 2 plants exhibited normal vegetative organ growth and female fertility but showed complete male sterility with smaller anthers lacking pollen grains, compared with WT ( Figure 4D). SEM assay revealed zmfar1 mutant anthers lacked the knitting cuticle and Ubisch bodies ( Figure 4E). These results confirmed that ZmFAR1 is indispensable to male fertility in maize. Cosegregating molecular markers were developed and successfully used to detect zmfar1 mutant genotypes via PCR amplification with gel electrophoresis ( Figure 4F, Table S2). also cotransformed with the CFP-AtROP10 as a plasma membrane marker in tobacco leaf cells. The 35S-GFP vector was used as a negative control. Scale bars, 8 µm (B1) and 40 µm (B2).

CRISPR/Cas9-Based Confirmation of ZmFAR1 Function in Controlling Maize Male Sterility
ZmFAR1 (GRMZM2G120987/Zm00001d048337) is the maize ortholog of GMS genes OsDPW and AtMs2 encoding fatty acyl reductases (FARs) [34,37]. One CRISPR/Cas9 vector of ZmFAR1 was constructed and transformed into maize ( Figure 4A). Three homozygous zmfar1 mutants (ZmFAR1-Cas9-1, -2 and -3) with 198-bp, 27-bp, and 3-bp deletions in the first exon, respectively ( Figure 4B,C), were produced for the evaluation. zmfar1 mutant F2 plants exhibited normal vegetative organ growth and female fertility but showed complete male sterility with smaller anthers lacking pollen grains, compared with WT ( Figure 4D). SEM assay revealed zmfar1 mutant anthers lacked the knitting cuticle and Ubisch bodies ( Figure 4E). These results confirmed that ZmFAR1 is indispensable to male fertility in maize. Cosegregating molecular markers were developed and successfully used to detect zmfar1 mutant genotypes via PCR amplification with gel electrophoresis ( Figure 4F, Table S2).

Spatiotemporal Expression Analysis of ZmFAR1 Gene and Subcellular Localization and Enzymatic Activity Analysis of ZmFAR1 Protein
qRT-PCR analysis showed ZmFAR1 expression was undetectable in vegetative organs (root, stem, and leaf) and female organs (immature cob and silk), but was specific in anthers from stages 8b-9 to 9-10 with the peak at stage 9-10 ( Figure 5A). To determine the subcellular localization of ZmFAR1 in maize, the ZmFAR1-GFP construct was transiently expressed in maize protoplasts, and ZmFAR1-GFP fusion protein was colocalized with autofluorescence of chlorophyll in plastids ( Figure 5B), indicating that ZmFAR1 functions as a plastid-localized fatty acyl reductase.

Spatiotemporal Expression Analysis of ZmFAR1 Gene and Subcellular Localization and Enzymatic Activity Analysis of ZmFAR1 Protein
qRT-PCR analysis showed ZmFAR1 expression was undetectable in vegetative organs (root, stem, and leaf) and female organs (immature cob and silk), but was specific in anthers from stages 8b-9 to 9-10 with the peak at stage 9-10 ( Figure 5A). To determine the subcellular localization of ZmFAR1 in maize, the ZmFAR1-GFP construct was transiently expressed in maize protoplasts, and ZmFAR1-GFP fusion protein was colocalized with autofluorescence of chlorophyll in plastids ( Figure 5B), indicating that ZmFAR1 functions as a plastid-localized fatty acyl reductase.  Enzymatic activity analysis of ZmFAR1 was based on the prokaryotic expression system and in vitro activity assay. ZmFAR1 contains a NAD(P)H-binding Rossmann-fold domain with two conserved motifs GGTGFLA and YVFTK at the N terminus and a FAR_C domain at the C terminus. Alignment of ZmFAR1 and its putative orthologs in plants revealed that there were four conserved residues, i.e., G101 and G104 in motif I (GGTGFLA) and Y327 and K331 in motif II (YVFTK) ( Figure 5(C1)). The MBP-tagged ZmFAR1 protein (ZmFAR1-MBP) was expressed in Escherichia coli and confirmed by Western blotting (Figure 5(C2)). Then, using three fatty acyl-CoAs (C12:0-, C16:0-and C18:0-CoAs) as substrates, the enzymatic activities of ZmFAR1 were investigated in vitro. Based on a spectrophotometric assay following the consumption of NADPH, ZmFAR1 catalyzed the reduction of the three fatty acyl-CoAs ( Figure 5(C3)) and showed significantly higher catalytic activity to C12:0-CoA than C16:0-and C18:0-CoAs ( Figure 5(C4)), suggesting that ZmFAR1 has enzyme activities to multiple fatty acyl-CoA substrates with various chain lengths (C12 to C18). In addition, the effects of four mutations (G101A, G104A, Y327F, and K331I) on substrate activities were also evaluated and compared with WT ZmFAR1. Almost all the four variant enzymes showed significantly decreased specific activities to the three fatty acyl-CoAs ( Figure 5D), indicating that all the substituted residues at positions 101, 104, 327, and 331 can significantly impair the reduction capability of ZmFAR1 toward different substrates.

Functional Comparison of ZmABCG26 and ZmFAR1 Contributing to Lipid Metabolism in Maize Anthers
The cytological defects of anther cuticles (Figure 2FandFigure 4E) indicated that zmabcg26 and zmfar1 mutations may alter lipid biosynthesis or transport for anther cuticle formation. Anthers of WT and two mutants collected before anthesis were used to detect the components of cutin, wax, and internal lipid by using gas chromatography-mass spectrometry (GC-MS).
Lipidomics analysis showed that the total anther cutin contents in zmfar1 (0.345 µg/mm 2 ) and zmabcg26 (0.095 µg/mm 2 ) were significantly lower (p < 0.001) than that in WT (0.793 µg/mm 2 ), and the total cutin content in zmabcg26 anther was significantly lower (p < 0.01) than that in zmfar1 anther ( Figure 6A(1)), indicating zmabcg26 mutation had a stronger effect on the reduction of anther cutin content. The differences of total cutin content among WT, zmfar1, and zmabcg26 anthers were mainly due to the alterations of 21 cutin monomers, including seven reduced (91.35% of the total cutin), seven increased (1.29% of the total cutin), and seven oppositely or differently changed (7.36% of the total cutin) ( Figure 6A(2), Figure S2). Notably, the content of C18:2 FA (5.57% of the total cutin) was significantly increased in zmfar1 anther but decreased in zmabcg26 anther. different alterations between zmfar1 and zmabcg26 anthers. For example, the amounts of C29, C32, and C34 ALKs (18.44% of the total wax) and C20 FA content (0.16% of the total wax) were significantly increased only in zmabcg26 anther, while C33 ALK amount (2.56% of the total wax) was significantly reduced only in zmfar1 anther (Figures 6(B2) and S3).  Figure 6(A1,B1,C1), ** and *** indicate the significant levels of 1% and 1‰ (Student's t test, n = 3), respectively. In Figure 6(A2,B2,C2), the red or blue lines represent the increase or decrease of the cutin, wax, and internal lipid monomers in mutant anthers, respectively. The thickness of the lines represents the change magnitude of the increase or decrease. The total wax contents in zmfar1 and zmabcg26 anthers (0.067 and 0.077 µg/mm 2 , respectively) were significantly higher (p < 0.001) than that (0.048 µg /mm 2 ) in WT anther (Figure 6(B1)). Although the total wax content of zmabcg26 was seemingly higher than that of zmfar1, the difference was not statistically significant. Eight wax constituents, C21 and C24-28 alkanes (ALKs) (19.64% of the total wax) and C24 and C28 FAs (undetectable in WT), were greatly increased in both zmfar1 and zmabcg26 anthers when compared to those in WT anther, which resulted in significantly higher contents of total wax in the two mutant anthers (Figure 6(B2), Figure S3). Interestingly, some of the wax monomers displayed different alterations between zmfar1 and zmabcg26 anthers. For example, the amounts of C29, C32, and C34 ALKs (18.44% of the total wax) and C20 FA content (0.16% of the total wax) were significantly increased only in zmabcg26 anther, while C33 ALK amount (2.56% of the total wax) was significantly reduced only in zmfar1 anther (Figure 6(B2) and Figure S3).
In addition, no significant difference in the total internal lipid contents was found between zmfar1 (16.44 µg/mg DW) and WT anthers (16.62 µg/mg DW), while the total internal lipid in zmabcg26 anther (4.81 µg/mg DW) was greatly decreased when compared with those in WT and zmfar1 anthers (Figure 6(C1)), suggesting that ZmFAR1 and ZmABCG26 play significantly different roles in the metabolism of internal lipid in maize anthers. Notably, the contents of four internal lipid constituents (C16:1, C18:1, C18:2, and C20 FAs, 23.26% of the total internal lipid) were significantly increased in zmfar1 anthers but significantly decreased in zmabcg26 anthers. Inversely, the level of C14 FA (0.18% of the total internal lipid) was significantly decreased in zmfar1 anthers but significantly increased in zmabcg26 anthers ( Figure 6(C2) and Figure S4).
Therefore, lipidomics analysis indicates that both ZmABCG26 and ZmFAR1 participate in anther lipid metabolism that is necessary for anther and pollen development in maize, but they play significantly different roles in the metabolism of cutin, wax, and internal lipid in maize anthers.

ZmFAR1 and ZmABCG26 May Be Regulated by a Complex TF Regulatory Network during Anther Development
In addition to miRNAs, the expressions of lipid metabolic GMS genes are usually regulated by TF GMS genes, and also by feedback regulation of other lipid metabolic GMS genes. The expression pattern changes of ZmFAR1 and ZmABCG26 in anther development were analyzed by using seven sets of maize anther transcriptomes, including three TF GMS mutant lines (ms23, ms7-6007, and lob30) [12,42,44,45], two TF overexpression GMS lines (p5126-ZmMs7 and p5126-ZmLOB30) [44], and two lipid metabolic GMS mutant lines (ms30 and ms33) [7,[46][47][48] (Figure S5A,B, Table S3). In ms23 mutant, ZmFAR1 expression was not significantly changed, while ZmABCG26 expression was significantly downregulated at stages 6 and 7 ( Figure S5(A1), indicating TF ZmMs23 may positively regulate ZmABCG26 expression. The expression of ZmFAR1 displayed downregulation in ms7-6007 mutant and upregulation in p5126-ZmMs7 overexpression line, while that of ZmABCG26 was downregulated in ms7-6007 mutant and its overexpression line (p5126-ZmMs7) ( Figure S5(A2),(A3)), reflecting that TF ZmMs7 may positively regulate ZmFAR1 and can affect ZmABCG26 expression. The expressions of ZmFAR1 and ZmABCG26 were upregulated in lob30 mutant, but downregulated in p5126-ZmLOB30 overexpression line ( Figure S5(A4),(A5)), suggesting that TF ZmLOB30 may negatively regulate ZmFAR1 and ZmABCG26. ZmFAR1 and ZmABCG26 were not differentially expressed between WT and ms30 or ms33 mutants, showing that their expressions are not affected by lipid metabolic genes ZmMs30 and ZmMs33 ( Figure S5B). These different or consistent expression pattern changes indicate that ZmFAR1 and ZmABCG26 may be regulated by a complex TF regulatory network.
The expression level of ZmABCG26 was significantly changed in zmfar1 anthers, compared with WT anthers at stages 8b and 9 ( Figure S5C), implying the lipid biosynthetic gene ZmFAR1 may affect the expression and function of the lipid transport gene ZmABCG26 during anther development.

miRNA-Mediated Lipid Metabolism Is Most Probably Critical for Maize Anther Development
In maize, though a certain number of miRNA regulatory modules and networks have been revealed by bioinformatics analyses [12,42], further experimental studies are required to confirm that the important roles of miRNAs have on regulating the expression of lipid metabolic genes during anther development. The negatively correlated expression pattern between zma-miR164h-5p and ZmABCG26, a lipid transporter gene, revealed in the present study, indicates the potential regulatory role of zma-miR164h-5p on ZmABCG26 expression. ZmABCG26 was demonstrated as a maize GMS gene by CRISPR/Cas9-mediated mutagenesis (Figure 2). These results suggest that the predicted miR164h-ZmABCG26 module is most probably necessary for maize male reproductive process by regulating lipid transport, which supports the proposed miRNA-mediated regulatory module in anther lipid metabolism underlying plant male fertility. In addition to the known miRNA, several miRNA candidates were discovered by analyzing anther transcriptome data from different genetic backgrounds or mutant lines, and some of the miRNA candidates (e.g., miRN2025 in Figure 1B) have high expression levels and stage-specific expression patterns similar to those of the known miRNAs during anther development. Importantly, the miRNA candidates were predicted as potential expression regulators of ABCG family genes in maize ( Figure 1A). This result implies that the newly predicted miRNA candidates may also contribute to anther development and pollen formation in maize. Nevertheless, the knockout/down and overexpression experiments of these miRNA loci as performed in some previous studies [11,13,14] should be required in future investigation to obtain more direct lines of evidence supporting the critical roles of the miRNA-ABCG gene modules for male reproduction in maize and other plant species.

ZmFAR1 and ZmABCG26 Display the Functional Differences in Anther Lipid Metabolism
Two lipidic barriers, anther cuticle and pollen exine covering the outer surfaces of anther and pollen grain, respectively, play important roles in protecting male reproductive development in flowering plants. The formation of anther cuticle and pollen exine shares a part of lipid metabolic pathways in the anther tapetum [4]. Accordingly, the tapetum is considered as the main place to produce cuticle and sporopollenin precursors that are transported to the anther or pollen surfaces for polymerization. Many tapetum-expressed lipid metabolic GMS genes have been reported to be involved in anther cuticle and pollen exine development, such as ZmMs10/APV1 [49], ZmMs20/IPE1 [50,51], ZmMs26 [52], ZmMs30 [46], and ZmMs33 [47,48] in maize, and OsCYP703A3 [53], OsCYP704B2 [54], OsDPW [37], OsABCG15 [26], and OsABCG26 [27] in rice. In these GMS mutants, the lipidic contents and compositions of the anther cuticular layer are significantly altered, while comparing the impact of different lipid metabolic GMS genes on metabolic profiles of cutin and wax has been rarely reported.
Anther cuticle mainly consists of cuticular cutin and wax in which the lipidic polyester cutin lays on the anther surface, while the cuticular wax is embedded in the cutin matrix [55]. Cutin mainly consists of FA oxygenated derivatives with chain lengths of C16 and C18. In the present study, cutin amounts in zmfar1 and zmabcg26 anthers were significantly decreased ( Figure 6A and Figure S2), which was similar to results reported in the mutants of their orthologous genes, OsDPW [37] and OsABCG15 [26], indicating conserved roles in anther cutin synthesis and transport in monocot plants. Among the major cutin monomers, hydroxy FAs (C16 ωHFA and C18-9,10 DHDFA) and unhydroxylated FAs (C16 and C18 FAs) were significantly decreased in contents in zmfar1 and zmabcg26 anthers, and the change trends of these components were also observed in previously reported lipid metabolic mutants, such as zmms30 [46], zmms10/apv1 [49], and zmms20/ipe1 [50,51], suggesting the functional conservation of these GMS genes in cutin synthesis and transport. Notably, the contents of five FAs (C18, C20, C22, C18:1, and C18:2 FAs) were increased only in zmfar1 anther ( Figure 6A and Figure S2), indicating functional divergences among lipid metabolic genes in controlling anther cuticle formation. In plant anthers, de novo FA (up to C18) synthesis occurs in plastids where AtMs2 and OsDPW, the orthologs of ZmFAR1, function as FA reductases in the production of fatty alcohols. The content increases of C18, C18:1, and C18:2 FAs in zmfar1 anthers may directly result from the deficient function of ZmFAR1, which leads to the inhibition of FA reduction in plastids. Cuticular wax is a complex mixture of very long-chain FA derivatives [56], including ALKs, FAs, and fatty alcohols. The total wax contents in zmfar1 and zmabcg26 anthers were significantly increased ( Figure 6B and Figure S3), which were predominantly due to significant increases in the contents of unsaturated alkanes, such as C24-C28, suggesting that the loss function of ZmFAR1 and ZmABCG26 may contribute positive feedback regulation of wax biosynthesis and transport. However, the amounts of C29, C32, and C34 alkanes were increased only in zmabcg26 anthers, suggesting different positive feedback regulation mechanisms in zmfar1 and zmabcg26. Similarly, a slight increase in the total amount of wax was also detected in rice osdpw anthers [37], reflecting the functional conservation.
ZmFAR1 and its orthologs, AtMs2 and OsDPW, are consistently localized in plastids and required for lipid biosynthesis in plant anthers. However, ZmABCG26 essential for lipid transport is primarily localized in ER, plasma membrane, and plastids ( Figure 3B), which is not completely consistent with the subcellular localization of its orthologs, AtABCG26 and OsABCG15, which are localized on the plasma membrane [18,19,26], implying that the functions of ZmABCG26 and its orthologs may be slightly different between plant species.
The GMS mutants zmabcg26 and zmfar1 can be used for developing biotechnologybased male sterility (BMS) systems and creating male-sterile lines, which are critical for maize hybrid breeding and seed production. The cosegregating molecular markers developed based on CRISPR/Cas9-induced mutation sites or fragments in mutants zmabcg26 and zmfar1 (Figures 2G and 4F) are beneficial to determine the mutant genotypes of crossing or backcrossing plants before flowering and pollination so that cross or backcross breeding can be efficiently performed to create male-sterile lines in different genetic backgrounds.

Plant Materials and Growth Conditions
Maize inbred lines B73, Oh43, and Zheng58 and hybrid Hi II were used in this study.

Characterization of Mutant Phenotypes
Images of the tassels and anthers were captured with a Canon EOS 700D digital camera (Canon, Tokyo, Japan) and a SZX2-ILLB stereomicroscope (Olympus, Tokyo, Japan) respectively. Routine analysis of I 2 -KI staining and SEM (scanning electron microscopy) were performed as described previously [44]. SEM images of anther and pollen grain were detected with a HITACHI S-3400N scanning electron microscope (HITACHI, Tokyo, Japan).

Plasmid Construction and Maize Transformation
For generating CRISPR/Cas9-mediated mutants, the CRISPR/Cas9 plasmids were constructed based on the pBUE411 vector as described previously [66]. The CRISPR-P 2.0 (http://crispr.hzau.edu.cn/CRISPR2/, accessed on 13 May 2019) was used to choose specific gRNAs targeting the selected genes, and then the off-target analysis of gRNAs sequences was carried out on the website (http://www.rgenome.net/cas-offinder/) (accessed on 13 May 2019). A CRISPR/Cas9 plasmid contains two gRNAs. For assembly of the two gRNAs, a PCR fragment was amplified based on pCBC-MT1T2 with the primer pair specific for each gene (Table S2) and then purified PCR fragments were cloned into BsaI-digested pBUE411 vector.
Maize hybrid Hi II was used as the transformation receptor. The CRISPR/Cas9 vector was introduced in Agrobacterium tumefaciens strain EHA105 and transformed into immature embryos as described previously [67]. The Bar gene was used as a selectable marker and positive transformants were then selected by PCR amplification using primers OGF41 and OGF42. The primers were listed in Table S2.

Genotyping Maize Plants
Genomic DNA was extracted from leaves of maize seedlings for the T 0 , F 1, and F 2 generations using the cetyltrimethylammonium bromide (CTAB) method [68]. PCR amplifications of relevant regions were performed using the specific primers flanking the target sites listed in Table S2, and then the purified PCR products were cloned into the pEASY-T1 vector (TransGen Biotech, Beijing, China) for Sanger sequencing. The sequencing chromatograms were carefully examined for exact patterns that might indicate the mutation types, such as homozygous, monoallelic, or diallelic mutations. For genotyping plants from the F 2 generation, cosegregating molecular markers were designed according to mutations of target genes (Table S2) as described previously [51], and the PCR products were analyzed by agarose gel electrophoresis or polyacrylamide gel electrophoresis.

Quantitative Real-Time PCR (qRT-PCR) Analysis
Total RNA was isolated using TRIzol reagent (Invitrogen, Waltham, MA, USA), and genomic DNA was removed with DNase I (Promega, Madison, WI, USA). The cDNA was synthesized using 5X All-In-One RT MasterMix (abm, Richmond, Canada) following the manufacturer's protocol. The reverse transcription of microRNA was performed using PrimeScript RT reagent Kit (TaKaRa, Kusatsu, Japan) according to Chen et al. [69]. qRT-PCR was conducted with the corresponding primer set (Table S2) on a QuantStudio 5 Real-Time PCR system (ABI, Waltham, MA, USA) using TB Green Premix EX Tag (TakaRa, Kusatsu, Japan). Both ZmUbi2 and ZmCyanase were used as two internal controls. U6 snRNA was used as the housekeeping RNA in the expression analysis of microRNA. Each sample had three technical replicates and three biological replicates. The amplification data were calculated by the 2 −∆∆Ct method [70], and quantitative results were shown as means ± SD.

Phylogenetic Analysis
WBC/ABCG family members in maize, rice, and Arabidopsis thaliana were obtained from previous studies [16,71,72]. The phylogenetic tree was reconstructed in MEGA7.0 using the neighbor-joining method, with Poisson correction, pairwise deletion, and 1000 bootstrap replicates.

Subcellular Localization of ZmABCG26 and ZmFAR1
The corresponding full-length cDNAs of ZmABCG26 and ZmFAR1 were amplified using the primer pairs listed in Table S2 and fused in-frame with N terminal of GFP and then inserted to the expression vectors pUC19 and pJG186 via infusion method, respectively. The ER marker mCherry-HDEL was created by adding the sequence encoding HDEL to the 3 end of mCherry [73]. The vector of plasma membrane fluorescent marker CFP-AtROP10 was kindly provided by Yule Liu's Lab [74]. Constructs were transiently expressed in maize protoplasts or tobacco leaves. To detect the GFP-or mCherry-tagged proteins, an LSM780 confocal laser scanning microscope (Zeiss, Jena, Germany) was used.

Enzyme Activity Assay of ZmFAR1 and Its Mutants
Using a basic local alignment searching tool-protein to protein (BLASTP) search on the GRAMENE website (http://www.gramene.org/, accessed on 28 August 2020), the sequences of three orthologs of ZmFAR1 were obtained, and then sequence alignment of these orthologous proteins was performed to predict the conserved motifs (I and II) and active sites [75]. To generate the four mutants of ZmFAR1, the overlap-extension PCR method was performed with specific primers by substitutions of G101A, G104A, Y327F, and K331I. To construct the expression vectors, each of the full-length cDNAs of ZmFAR1 and the point mutation sequences was fused to the maltose-binding protein (MBP) gene and cloned into the vector pMAL-C2X (New England Biolabs, Ipswich, MA, USA). The fusion proteins of ZmFAR1-MBP, G101A-MBP, G104A-MBP, Y327F-MBP, and K331I-MBP were expressed in E. coli strain BL21, purified by using Ni affinity column, and then subjected by SDS-PAGE analysis with conformation by Western blot using MBP antibody (New England Biolabs, Ipswich, MA, USA). The primer sequences used for PCR amplification were listed in Table S2.

Analysis of Anther Wax, Cutin, and Internal Fatty Acids
Stage 13 WT, zmfar1, and zmabcg26 anthers were collected and freeze-dried using an α 2-4 LD plus freeze dryer (Christ, Osterode, Germany). Micrographs were collected to determine the surface area of anthers, assuming a cylindrical body for maize anthers as described previously [37]. Wax, cutin, and internal fatty acids were analyzed as described previously [46].