Ziwen Li,Toto Zhu,Shungshung Liu,Yilin Jing,Hoyun Liu,Yuwen Zhng,Ke Xie,,Jinping Li,Xueli An,,*,Xingyun Wn,,*

a Zhongzhi International Institute of Agricultural Biosciences,Biology and Agriculture Research Center,University of Science and Technology Beijing,Beijing 100024,China

b Beijing Engineering Laboratory of Main Crop Bio-Tech Breeding,Beijing International Science and Technology Cooperation Base of Bio-Tech Breeding,Beijing 100192,China

Keywords:

A B S T R A C T Anther development is a programmed biological process crucial to plant male reproduction.Genomewide analyses on the functions of transcriptional factor(TF)genes and their microRNA(miRNA)regulators contributing to anther development have not been comprehensively performed in maize.Here,using published RNA-Seq and small RNA-Seq(sRNA-Seq)data from maize anthers at ten developmental stages in three genic male-sterility(GMS)mutants(ocl4,mac1,and ms23)and wild type W23,as well as newly sequenced maize anther transcriptomes of ms7-6007 and lob30 GMS mutants and their WT lines,we analyzed and found 1079 stage-differentially expressed(stage-DE)TF genes that can be grouped into six(premeiotic,meiotic,postmeiotic,premeiotic-meiotic,premeiotic-postmeiotic,and meiotic-postmeiotic clusters)expression clusters.Functional enrichment combined with cytological and physiological analyses revealed specific functions of genes in each expression cluster.In addition,118 stage-DE miRNAs and 99 miRNA-TF gene pairs were identified in maize anthers.Further analyses revealed the regulatory roles of zma-miR319 and zma-miR159 as well as ZmMs7 and ZmLOB30 on ZmGAMYB expression.Moreover,ZmGAMYB and its paralog ZmGAMYB-2 were demonstrated as novel maize GMS genes by CRISPR/Cas9 knockout analysis.These results extend our understanding on the functions of miRNA-TF gene regulatory pairs and GMS TF genes contributing to male fertility in plants.

1.Introduction

In plants,anther development is important for male reproduction during which four somatic cell layers(epidermis,endothecium,middle layer,and tapetum)of anther wall are differentiated and developed to support the development of inner microspore mother cells into mature pollen grains.Various biological activities are involved in this process,including mitotic and meiotic cell divisions,lipid metabolism,reactive oxygen species(ROS)homeostasis,starch turnover and photosynthesis in endothecium chloroplasts[1–3],and so on.Genes that have a decisive influence on male fertility are called GMS(genic malesterility)genes,and more than a hundred of GMS genes have been cloned and characterized mainly inArabidopsis thaliana(At),rice(Oryza sativa,Os)and maize(Zea mays,Zm)[3,4].Besides CMS(cytoplasmic male-sterility)genes,GMS genes are an important genetic basis of plant male sterile lines by which crop heterosis can be applied to increase crop yield.Maize has the highest total grain yield in the world by using maize hybrid.Discovering GMS genes and investigating their molecular regulatory mechanisms are necessary to elucidate plant male reproductive processes and promote their application in crop breeding and hybrid seed production.

Transcriptional factors(TFs)are major regulators in gene regulatory networks controlling developmental and stress response processes,including plant anther development and male fertility.In maize,a total of 17 GMS genes have been cloned,and six of them encode TFs,includingOUTER CELL LAYER4(OCL4)belonging to homeobox-leucine zipper(HD-ZIP)gene family[5],MALE STERILITY 32(ZmMs32)andZmMs23belonging to basic helix-loop-helix(bHLH)gene family[6,7],ZmMs9belonging to myeloblastosis(MYB)gene family[8],ZmMs7encoding a plant homeodomain(PHD)finger TF[9,10],andINDETERMINATE GAMETOPHYTE 1(IG1)belonging to lateral organ boundaries domain(LBD)gene family[11,12].InArabidopsisand rice,more than 26 GMS TF genes have been identified and functionally studied[4],and a core regulatory pathway controlling tapetum development and pollen formation has been established in the two species[13–15].Specifically,theArabidopsiscore regulatory pathway is composed of five TF genes includingDYSFUNCTIONAL TAPETUM 1(DYT1)[16],TAPETAL DEVELOPMENT AND FUNCTION 1(TDF1)[17],ABORTED MICROSPORE(AtAMS)[18],AtMYB103[19],andAtMS1[20].In rice,the orthologous genes of the fiveArabidodpsisTF genes have been cloned,includingUNDEVELOPED TAPETUM 1(UDT1)[21],OsTDF[22],TAPETUM DEGENERATION RETARDATION(TDR)[23],OsMYB80[24],andPERSISTANT TAPETAL CELL 1(PTC1)[25],and the loss-of-function mutant of each gene exhibited male sterility.Furthermore,the transcriptional regulatory pathway reported inArabidodpsisfor tapetum development is relatively conserved in rice[15],and may be also conserved in maize[4].In addition to the above five GMS TF genes,other GMS TF genes are also important regulators controlling anther development.For example,rice TF geneTDR INTERACTING PROTEIN2(TIP2)andArabidopsis bHLH010,bHLH089,andbHLH091are important nodes in the GMS gene regulatory network[26,27],and maizeZmMs23is their ortholog[7].Therefore,GMS TF genes set up the framework of complex regulatory networks regulating anther development and male fertility in plants.

miRNAs are small non-coding RNA regulators binding and/or cleaving target transcripts to post-transcriptionally regulate the gene expression levels or translation processes of target genes.Function deficiencies of some miRNAs may lead to plant male sterility due to disordered regulations on the target TF genes.For example,the maizezma-miR172loss-of-function mutant obviously lacked male characteristics and displayed male sterility[28].InArabidopsis,function deficiency ofath-miR319aresulted in abnormal stamen due to disrupted regulation onTCP4gene encoding a teosinte-branched 1/cycloidea/proliferating(TCP)TF[29].

miR160-ARF,miR166-PHB(a HD-ZIP TF gene),miR319-TCP24,andmiR396-GRF(a growth-regulation factor TF gene)were crucial regulatory modules duringArabidopsisanther development[30–33].In rice,miR159targetedOsGAMYBregulating anther development,and themiR159-MYBregulatory module was also confirmed inArabidopsis[34,35].Meanwhile,themiR164-NACmodule was important in rice anther development and male fertility[36].In tomato(Solanum lycopersicum,Sl),anthers with silencedmiR171displayed disordered tapetum development and reduced callose deposition around the tetrads that were induced by the deregulated expression ofmiR171target gene,SlGRAS24(a GRAS TF gene)[37].These results indicate that miRNAs are of great importance in anther development by their post-transcriptional regulation roles on TF genes that most probably function as GMS genes in plants.However,when compared to the obtained results inArabidopsisand rice,the regulatory roles of miRNAs on anther development and male fertility are not well unveiled in maize.Therefore,more exploring work on the computational identification and experimental confirmation of TF genes and their miRNA regulators potentially associated with maize anther development are required.

High-throughput RNA-Seq and sRNA-Seq technologies and the comparative transcriptomics analysis can be effectively combined to investigate the functions and regulations of TF genes and miRNAs.Previous genome-wide studies have revealed part of the potential functions of miRNAs and other regulatory factors on maize anther development[38,39],whereas transcriptome data used in these studies was relatively limited.Therefore,the comparative transcriptomics analysis on TF genes and miRNAs has not been performed across the anther whole developmental process and different genetic backgrounds in maize.Fortunately,the available data set including RNA-Seq and sRNA-Seq data from maize anthers sampled at ten stages[40],largely representing the anther whole developmental period,provides a good opportunity to address this question.

Here,we sequenced miRNA transcriptome of maizems7-6007mutant(Oh43 genetic background)anthers,as well as mRNA and miRNA transcriptomes oflob30mutant(B73 genetic background)anthers.Both mutants displayed complete male sterility without pollen grain.By transcriptome analyses on published and our newly sequenced data in maize anthers and experimental verification,we found a number of stage-DE TF genes with annotated functions well corresponding to the cytological and physiological alterations during anther development,and detected 118 miRNAs and 99 miRNA-TF gene pairs associated with maize anther development.Bothzma-miR319andzma-miR159,as well as ZmMs7andZmLOB30 were found to regulate or influence expression ofZmGAMYB,which was further verified to be critical for maize male fertility by gene knockout analysis based on the CRISPR/Cas9 system.The obtained results expand and deepen our understanding on the functions and regulations of TF genes and their miRNA regulators controlling plant male reproductive development.

2.Materials and methods

2.1.Plant materials and phenotypic characterization

All maize plants were grown in the experimental stations of USTB in Beijing and Hainan province.The T0transgenic plants were greenhouse-grown in Beijing.Morphological images of maize tassels were captured by a Canon EOS 700D digital camera(Canon,Tokyo,Japan).Spikelets were photographed with a SZX2-ILLB stereomicroscope(Olympus,Tokyo,Japan).Pollen grains were stained with 1%I2-KI solution and imaged by a BX-53F microscope(Olympus,Tokyo,Japan).

2.2.Cytological observation and microscopy

Fresh anther samples for semi-thin transverse section and scanning electron microscopy(SEM)analyses were fixed in FAA solution(Coolaber,Beijing,China),the subsequent procedures were performed as described previously[41,42].

2.3.Chloroplast development analysis

Chloroplast isolation and observation were performed as described previously[1].Briefly,fresh maize anthers at different developmental stages and the fresh maize leaves were collected,and then chloroplasts were isolated using a chloroplast isolation kit(product code CP-011;Invent).Chloroplasts were observed under a Leica TCS-SP8 confocal laser scanning microscope(Leica,Wetzlar,Hessen,Germany).

2.4.RNA-Seq and sRNA-Seq

Total RNA was extracted from anthers using an RNAprep Pure Plant Kit(Polysaccharides and Polyphenolics-rich)(TIANGEN,Beijing,China).The ribosomal RNA was removed by Epicentre Ribozero rRNA Removal Kit(Epicentre,Wisconsin,USA).RNA-Seq libraries were generated using the rRNA-depleted RNA by NEBNext UltraTM Directional RNA Library Prep Kit for Illumina(NEB,Massachusetts,USA).After adapter ligation and cDNA synthesis,the libraries were sequenced on an Illumina Hiseq 4000 platform to generate 150 bp paired-end reads.The sRNA-Seq libraries were constructed using NEBNext Multiplex Small RNA Library Prep Set for Illumina(NEB,USA).After adapter ligation and cDNA synthesis,140–160 bp DNA fragments(corresponding to 18–30 nt small RNA)were recovered by 8%PAGE.The libraries were sequenced on an Illumina Hiseq 2500 platform to generate 50 bp single-end reads.The summary information of RNA-Seq and sRNA-Seq data used here was listed in Table S1.

2.5.RNA-Seq data and TF gene analysis

RNA-Seq data was analyzed according to the workflow described previously[38,43].Raw data were processed by using NGSQCToolkit[44]to filter out adaptor sequences and low quality reads,and the remaining clean reads were mapped to the maize reference genome(version AGPv4 downloaded from Ensembl release 37)by using TopHat2[45],gene expression levels were estimated by using Rsubread[46],and DE genes(DEGs)were identified using edgeR[47].Genes with fold changes more than two times and adjusted significant levels smaller than 0.05 were considered as DE genes.The genome-wide annotation of maize TF genes was downloaded from the Plant Transcription Factor Database(http://planttfdb.cbi.pku.edu.cn/).Co-expression analysis of stage-DE TF genes was performed by K-means clustering method measured by within-cluster sum of squared errors.Genes coexpressed with stage-DE TF genes were assigned to the 60 expression patterns(Fig.S1)according to their highest expression correlation coefficients related to the 60 types(Pearson correlationr＞0.9).Gene family enrichment analysis was performed based on the sizes of TF gene families listed in Table S2.GO enrichment analysis was performed by agriGO web-based tool[48].Pathway enrichment analysis was performed based on KEGG database[49].

2.6.sRNA-Seq data and miRNA analysis

sRNA-Seq data was analyzed according to the method described previously[38,43].Generally,short reads were processed and filtered by Infernal package and Rfam database[50,51],and the clean reads were mapped to the maize reference genome by Bowtie[52].The expression levels of known miRNAs and the DE miRNAs were estimated and identified by the same method used in the RNA-Seq analysis.The annotation information of maize known miRNAs was extracted from miRBase[53].

2.7.miRNA-TF gene pair analysis

The target TF genes of miRNAs were predicted by tapir and psRobot programs according to the approach described previously[38].Both predicted results of tapir and psRobot programs were used.Potential miRNA-TF regulatory pairs were identified by using two methods.In method one,miRNA-TF pairs with significant negative correlations(Pearson correlation,P＜0.05)in expression levels in any of the six transcriptome data(W23 line andocl4,ms23,ms7,lob30andmac1mutant lines)across their investigated stages(Table S1)were considered as potential miRNA-TF gene regulatory pairs.In method two,miRNA-TF pairs with a half of investigated stages showing opposite expression changes in transcriptomes for miRNAs and TF genes were considered as potential regulatory pairs.Two sets of investigated stages were used,including 15 stages forocl4,ms23andmac1mutants and 7 stages forms7andlob30mutants.

2.8.qRT-PCR analysis

qRT-PCR was performed as described previously[54].For miRNA quantification,stem-loop RT-PCR was used as described previously[55].Primers for each gene and stem-loop assay are listed in Table S3.

2.9.CRISPR/Cas9 mutagenesis and maize genetic transformation

For site-directed mutagenesis ofZmGAMYB(Zm00001d012544)and its paralogZmGAMYB-2(Zm00001d043131), the CRISPR/Cas9 plasmid was constructed based on the pBUE411 vector as described previously [56]. For assembly of two gRNAs, PCR fragment was amplified from pCBC-MT1T2 with primers sets(ZmGAMYB-1-target1-F and ZmGAMYB-2-target2-R) among which target sites (T1: GCCCCTGCGAGGCGACCTG, T2: AGGCCCAACCTCAAGAAGG) were incorporated, and then, the purified PCR fragment was digested withBsaI and cloned into the pBUE411 vector. The colonies were confirmed by sequencing and the corresponding vector was designated as pCas9-2 g-GAMYB (Fig. S2).

Immature zygotic embryos of the maize(Zea mays)Hi II hybrid genotype were used as genetic transformation materials.Agrobacterium-mediated genetic transformation was performed according to the previous procedures[57].The pCas9-2g-GAMYB transformants were screened by PCR amplification using primers OGF-41/42 specific toBargene.For the knockout transformants,sequencing analyses ofZmGAMYBandZmGAMYB-2genomic fragments were performed using primers ZmGAMYB-F/R and ZmGAMYB-2-F/R.

3.Results

3.1.Stage-differentially expressed TF genes associated with anther development in maize

Anther development beginning with anther primordia and archesporial cell specification and ending with pollen emission is a multistage and programmed biological process(Fig.1A),which is of great importance for plant male reproduction.In maize,the previously published transcriptome data are valuable for maize anther transcriptomics analysis[40].It contains anther samples at ten different anther lengths including 0.2,0.4,0.7,1.0,1.5,2.0,2.5,3.0,4.0 and 5.0 mm,representing the key developmental processes of W23 maize anthers such as cell fate specification,cell differentiation,meiosis,uninucleate microspore,binucleate microspore,gametogenesis and pollen mature based on the cytological observations(Fig.1A,B).Clustering analysis of expressed genes revealed that these transcriptomes of the ten stage anthers were grouped into three categories(Fig.1B),including early stages or premeiotic stages(with anther lengths of 0.2 to 1.0 mm responding stages 2 to 5),middle stages or meiotic stages(with anther lengths of 1.5 to 3.0 mm responding stages 6 to 9),and late stages or postmeiotic stages(with anther lengths of 4.0 and 5.0 mm responding stages 10 to 13).The three categories represented the major developmental framework of maize anther,in which anther somatic cells were divided into four cell layers(epidermis,endothecium,middle layer,and tapetum)composing anther wall and archesporial cells were formed and divided at early stages(S1–S5)(Fig.1A),microspore mother cells were divided to form microspores at middle stages(S6–S9)(Fig.1A),and microspores further developed to form mature pollen at late stages(S10–S14)(Fig.1A).These results indicate that the transcriptomic changes well correspond to the cytological alterations during the whole developmental process of maize anther.

Fig.1.Cytological observation of maize anthers and analysis of stage-differentially expressed TF genes during anther development.(A)Cytological observation of maize anther during developmental stages 1 to 14.Ar,Archesporial cell;BP,Bicellular pollen;C,Connective tissue;Dy,Dyad cell;E,Epidermis;En,Endothecium;L1,L2,and L3,the three cell layers in stamen primordia;MC,Meiotic cell;ML,Middle later;MMC,Microspore mother cell;MP,Mature pollen;Msp,Microspore parietal cell;PPL,Primary parietal layer;Sp,Sporogenous cell;SPL,Secondary parietal cell layer;Ta,Tapetum;Tds,Tetrads.Scale bars,50 μm.(B)Clustering result of maize anther RNA-Seq transcriptomes from anther samples of ten stages.The re-analyzed data and sample information were referred to[40].(C)Venn diagram of stage differentially expressed(DE)and expressed transcriptional factor(TF)genes,showing an overlap of 1087 stage DE TF genes.(D)Stage-DE TF genes were grouped into six expression clusters.Black lines represent the mean expression levels of included expression types in each cluster.(E)Gene family enrichment analysis of stage-DE TF genes in early(E),middle(M),late(L),early-middle(EM),early-late(EL)and middle-late(ML)clusters.FE,fold enrichment.The FE value of each TF family in the‘‘All”column was the ratio value(the number of expressed TF genes to TF family size)against the average value of these ratios in the 56 TF families.The FE value of each TF family in each of the other columns was the ratio value(the number of stage-DE TF genes to the number of expressed TF genes)against the average value of these ratios in the 56 TF families.

TFs play crucial roles in the regulation of biological processes.There are 39,005 protein-coding genes in the maize genome(AGP v4,Ensembl release 37),of which 2216 genes(5.7%)are predicted to encode TFs belonging to 56 gene families according to the annotation in plantTFDB(http://plantregmap.gao-lab.org/).Among the 56 TF gene families,42 had a family size with more than or equal to 10 members(Table S2).The top three TF gene families in the maize genome are bHLH,ethylene response factor(ERF)and MYB including 185, 173, and 167 members, respectively (Table S2). Here, we found 1301 of the 2216 TF genes had detect-able expression levels (read counts per kilobase per million mapped reads, RPKM ＞ 1) in at least one stage in anther transcrip-tome data of W23 line. The numbers of expressed TF genes in the top three TF families were 87, 80, and 81 for bHLH, ERF and MYB, respectively.In addition, the two smallest TF gene families were signal transducer and activator of transcription (STAT) and HRT-like (a novel DNA-binding protein first found in barley) with only one member of each family detected in anther transcriptomes. Furthermore, among the 39,005 protein-coding genes, 21,783 genes have detectable expression levels and 18,123 genes of them were identified as DEGs between any two of the ten investigated stages(stage-DEGs) in W23 anther transcriptomes, resulting in the identification of 1087 stage-differentially expressed (stage-DE) TF genes (Fig. 1C; Table S4). The ratio (6.0%, 1087/18123) of stage-DE TF genes among stage-DEGs is similar to that of TF genes among genome-wide protein-coding genes (5.7%) (Chi-square test,P= 0.1363), indicating it is a genome-wide transcriptome alteration of TF genes during anther development.

Through K-means clustering analysis,expression patterns of the 1087 stage-DE TF genes were grouped into 60 sub-clusters,59 of them were further integrated into six expression clusters including 1079 TF genes(Fig.S1).The six expression clusters included 294,187 and 44 stage-DE TF genes at early,middle and late stages,respectively,and 309,119 and 126 stage-DE TF genes at three combined stages containing early-middle,early-late and middle-late stages,respectively(Fig.1D,S1).Gene family enrichment analysis revealed that TF genes in the six expression clusters varied in enrichment patterns(Fig.1E;Table S2),possibly reflecting the functional preference of stage-DE TF genes at different stages during maize anther development.

3.2.Stage-DE TF genes and their co-expressed genes closely associated with the cytological and physiological processes during maize anther development

TFs and their regulated genes are usually co-expressed.To better understand the underlying molecular regulation of anther development,we first identified genes co-expressed with stage-DE TF genes.The expressed protein-coding genes other than stage-DE TF genes were assigned into the six expression clusters based on the highest expression correlations(Pearson’sr＞0.9)between the investigated genes and the 59 sub-clusters of expression patterns(Fig.S1).As a result,3886,1975,1106,2598,1115,and 944 protein-coding genes,were grouped into the six expression clusters(early,middle,late,early-middle,early-late and middle-late stages)as co-expressed genes of stage-DE TF genes,respectively.

To reveal the overall function of stage-DE TF genes and their coexpressed genes in each of the six expression clusters,we investigated their functions by gene ontology(GO)and KEGG pathway enrichment analyses(Figs.2A,S3),and further confirmed their roles by the cytological and physiological experiments,e.g.,observing the formation process of cutin and wax layer and Ubisch bodies on anther wall outer and inner surfaces,respectively,and of microspores and pollen grains in anther locule by SEM analysis(Fig.2B),observing the developmental process of anther endothecium chloroplasts by using confocal laser scanning microscope,and measuring photosynthetic activities of chloroplasts during anther development(Fig.2C).As a result,we found highly or preferentially expressed genes at early,early-middle,and early-late stages had overlapping enriched functions including regulation of gene expression/transcription and RNA biosynthesis(Fig.2A),indicating transcription activities are relatively high across the whole anther developmental process.For early stages(S1–S5),the preferentially expressed genes were exclusively enriched in biological processes including cell cycle checkpoint,cell proliferation,DNA replication,mitotic cell cycle and reproductive shoot system development(Fig.2A),and KEGG pathways of DNA replication,homologous recombination and base excision repair(Fig.S3).Meanwhile,the cytological observations during early stages suggested that anther primordial cells are divided and differentiated to form a four-layer anther wall and archesporial cells(Fig.1A).

The middle-stage DE genes had preferentially enriched GO functions in pollen exine formation and sporopollenin biosynthesis,and KEGG pathways of fatty acid(FA)biosynthesis and sugar/starch metabolism(Fig.2A,S3).Meanwhile,the SEM observations showed that a layer of cutin and wax was beginning to cover anther wall outer surface at stages 9 to 10,and Ubisch bodies and pollen exine started to form at the same stages(Fig.2B).The physiological assay displayed that anther endothecium chloroplasts progressively developed and their photosynthetic activities gradually increased from stages 7 to 9(Fig.2C).Additionally,the late-stage DE genes were functionally enriched in ion transport,nitrogen compound transport and transmembrane transport in GO analysis(Fig.2A),and several KEGG metabolic pathways such as FA degradation,amino acid biosynthesis,and carbon,starch and galactose metabolism(Fig.S3).By the SEM analysis,we found that the formation of a concentrated cutin and wax layer and the accumulation of starch in pollen grains occurred at late stages(Fig.2B).Notably,the middle-and late-stage DE genes had overlapping enrichment in oxidation-reduction process(Fig.2A).Similarly,the cytological observations showed the degradation of tapetal cells since stage 10(Fig.1A)that was induced by ROS accumulation and programmed cell death[1].Interestingly,the middle-late stage DE genes were preferentially enriched in photosynthetic processes in GO analysis(Fig.2A),and photosynthesisantenna proteins and carbon fixation in photosynthetic organisms in KEGG analysis(Fig.S3).We also found that anther endothecium chloroplasts were reached to the similar sizes of leaf chloroplasts at stage 10,and that anther chlorophyll fluorescence(Fv/Fmvalue)also reached to the similar level of leaf at stage 10(Fig.2C).From stages 11 to 13,anther endothecium chloroplasts kept the highly active photosynthesis in maize anthers(Fig.2C).Therefore,the enriched functions of stage-DE TF genes and their co-expressed genes could be largely confirmed by the cytological and physiological phenotypes of anthers during premeiotic,meiotic and postmeiotic developmental stages.

3.3.GMS TF genes control maize male sterility through regulating expression of a large number of TF genes during anther development

To further understand the roles and regulation of stage-DE TF genes, we sequenced and analyzed anther miRNA transcriptomes of WT andms7-6007GMS mutant with functional deficiency ofZmMs7gene (Zm00001d020680) that encodes a PHD TF [9,10]. AsZmMs7was grouped into the 49th expression pattern (Fig. S1,genes highly expressed at stages 8 and 9) belonging to the meiotic stages (Table S4), WT and mutant anthers at stages 8, 9, and 10 were sampled and sequenced (Fig. S4). Meanwhile, we also sequenced and analyzed mRNA and miRNA transcriptomes of another TF gene (ZmLOB30with the gene IDZm00001d036435)loss-of-function mutant (lob30) and its WT anthers (Fig. 3A–C,S5).lob30displayed no exerted anther and lacked pollen grain in the smaller and wilted anthers compared to its WT fertility plants(Fig. 3A). SEM analysis revealed smooth outer and inner surfaces of anther wall inlob30anther (Fig. 3B). Similar toZmMs7,ZmLOB30mainly expressed at stage 9 after the meiosis of microspore mother cell (Fig. 3D). Thus, WT andlob30anthers at developmental stages 8a, 8b, 9 and 10 were sampled and sequenced for RNA-Seq and sRNA-Seq analyses (Fig. 3C, S5). In addition, the published anther transcriptome dataset of two premeiotic GMS mutants,ocl4andms23,were also used for the following analysis.In total,anther transcriptomes of two premeiotic(ocl4andms23)[40]and two meiotic/postmeiotic(ms7-6007andlob30)TF gene mutants were used in this study.The expression patterns of the four GMS TF genes were investigated in WT anther during developmental stages 5 to 12 by qRT-PCR analysis(Fig.3D).

Fig.2.Cytological and physiological analyses of maize anthers and function annotations of stage-DE TF genes and their co-expression genes.(A)GO enrichment analysis of genes in the six expression clusters.FDR,false discovery rate.FE,fold enrichment.(B)SEM(scanning electron microscopy)analysis of outer and inner surface of anther wall as well as the microspore and pollen grains during maize anther developmental stages 6 to 13.The mesh structure on the outer surface of anther wall is cutin and wax layer(the representative structure at stage 12 was pointed with arrows).The dense dots on the inner surface of anther wall are Ubisch bodies(the representative structure at stage 11 was pointed with arrows).Scale bars,10 μm.(C)Chlorophyll fluorescence(C1)and chloroplast development(C2)analyses during maize anther developmental stages 6 to 13.Chlorophyll fluorescence and chloroplast development in mature leaf were used as control.Fv,variable fluorescence from dark-adapted tissue.Fm,maximal fluorescence from dark-adapted tissue.Fv/Fm reflects the potential maximum photosynthetic capacity of plant chloroplasts.Scale bars,1 mm in C1 and 1 μm in C2.

Fig.3.GMS-DE TF genes in maize anther transcriptomes of four GMS mutants(ocl4,ms23,ms7-6007 and lob30).(A)Phenotypic comparison of anthers between WT and lob30 mutant.(B)SEM(scanning electron microscopy)analysis of anthers between WT and lob30 mutant.(C)DEG(differentially expressed genes)analysis of anther transcriptomes between WT and lob30 mutant.(D)Expression patterns of four GMS TF(genic male-sterility transcriptional factor)genes(ZmOCL4,ZmMs23,ZmMs7 and ZmLOB30)during maize anther development.(E)Expression fold changes(FC)of GMS-DE TF genes between each of ocl4,ms23,ms7-6007 and lob30 mutant anthers and their corresponding WT line anthers.(F)Gene family enrichment analysis of GMS-DE TF genes between each of ocl4,ms23,ms7-6007 and lob30 mutant anthers and their corresponding WT line anthers.FE,fold enrichment.The FE value of each TF family in each column was the ratio value(the number of stage-DE TF genes to the number of expressed TF genes)against the average value of these ratios in the 56 TF families.

Compared to transcriptomes of their WT anthers,ocl4,ms23,ms7-6007, andlob30GMS mutant anthers had 6841, 8244, 3732,and 3663 DEGs (GMS-DEGs) across the investigated anther developmental stages, respectively. Among them, 439, 507, 269, and 295 DEGs were annotated as TF genes (GMS-DE TF genes) (Fig. 3E;Table S5). Notably,ms7-6007andlob30mutant anthers had higher proportions (7.2% and 8.1%) of GMS-DE TF genes among all GMSDEGs thanocl4andms23mutant anthers (6.4% and 6.2%, Chisquare test, allP＜ 0.01) which were similar to the proportion(6.0%, 1087/18123) of stage-DE TF genes among stage-DEGs in maize anther genome. It reflects that the regulated genes by GMS TFs ZmMs7 and ZmLOB30 may be more enriched in TF genes than those by GMS TFs ZmOCL4 and ZmMs23. Furthermore, among the GMS-DE TF genes inocl4,ms23,ms7-6007, andlob30mutant anthers, 90.1%, 92.0%, 93.9%, and 94.1% of them were stage-DE TF genes, which were significantly higher than the ratio of stage-DE TF genes among expressed TF genes during the W23 whole anther development (83.6%, 1087/1301) (Chi-square test, allP＜ 0.01)(Fig. 1C), possibly indicating the important roles of TF genes in maize anther development and male sterility.

To further investigate the functions and regulation of stage-DE TF genes,we analyzed the components and biological pathways of GMS-DE TF genes inocl4,ms23,ms7-6007andlob30mutant anthers by gene family and pathway enrichment analyses(Fig.3F,S6).The gene family enrichment results of GMS-DE TF genes varied in the four mutants(Fig.3F;Table S6).Notably,LBD TF family was enriched in all the four mutants,indicating LBD family genes may play critical roles in the gene regulatory networks controlling maize anther development.In KEGG analysis,we found that GMS-DE TF genes in the four mutants were all enriched in biosynthesis of secondary metabolites(Fig.S6).GMS-DE TF genes inocl4andms23mutants had overlapping enriched functions in plant hormone signal transduction,while those inms7-6007andlob30had overlapping enriched functions in phenylpropanoid and flavonoid biosynthesis,plant-pathogen interaction,and sugar metabolism.In addition,GMS-DE TF genes in each mutant specifically enriched in several pathways related to plant anther development(Fig.S6).These results suggested that stage-DE TF genes could participate in various pathways contributing to anther development in maize.

3.4.Stage-DE and GMS-DE miRNAs associated with anther development in maize

Many miRNAs are critical regulators mainly regulating TF genes at post-transcriptional level during anther development[58].To extend the gene regulatory networks in maize anther development,we analyzed small RNA transcriptomes of maize anthers and then identified potential miRNA-TF gene regulatory modules associated with anther development in maize.Among the 205 known miRNA groups in miRBase(miRNAs with the same mature sequences being integrated as one miRNA group in this study),118 groups expressed with CPM(read counts per million mapped reads)values＞1 in at least one sampled stage in W23 line anthers(Fig.4A,middle;Table S7).Among the 118 miRNA groups,20 were extremely highly expressed(at least one stage,the CPM value≥1000),27 were highly expressed(1000＞CPM≥100),29 were normally expressed(100＞CPM≥10),and 42 were low expressed(10＞CPM≥1)(Table S7).Among the extremely highly expressed miRNA groups, 10 groups includingzma-miR529-5p,zma-miR2118(a, b, c, d, and g),zma-miR166(k/j/n-3p, l/m-3p,and h/e/i/f/g/b/c/d-3p) andzma-miR319a/c/b/d-3p, eight groups includingzma-miR2275(b/c-3p and b-5p),zma-miR159(g/h/i-5p,b/j/k-5p, and a/b/f/g/k-3p),zma-miR408a/b-3pandzma-miR168(a/b-5p and a-3p), and two groups (zma-miR167g-3pandzmamiR396c/d) had expression peaks at early, middle, and late stages,respectively (Fig. 4A; Table S7). Similar to stage-DE TF genes, the 118 miRNA groups were clarified into three expression groups (I,early stage; II, middle stage; III, late stage) with 43, 43, and 32 miRNA groups, respectively (Fig. 4A, left). These stage-DE or highly expressed miRNAs may function in maize anther development by post-transcriptional regulation on their target genes.

As for miRNA expression in GMS mutant anthers, 24, 22, 12, and 31 miRNA groups were identified as GMS-DE miRNAs inocl4,ms23,ms7-6007andlob30mutant lines, respectively (Fig. 4B). Among all the above 67 independent GMS-DE miRNAs, 60 were stage-DE miRNA groups in W23 line anthers (marked with red cells in Fig. 4A, right). For the remaining seven GMS-DE miRNA groups,three (zma-miR395k-3p,zma-miR395e/h/j/p-5pandzma-miR395i-5p) were specifically expressed inocl4, two (zma-miR159c/d-3pandzma-miR393a-3p) in bothocl4andms23, one (zma-miR164f-5p)inms7-6007and one (zma-miR164h-5)inlob30mutant anther,respectively (Fig. 4B, marked in red). In summary, majority of GMSDE miRNAs were stage-DE miRNAs, which is similar to the high proportion of stage-DE TF genes among GMS-DE TF genes, implying potential functions of stage-DE miRNA groups during maize anther development.

3.5.Identification of miRNA-TF gene pairs potentially contributing to anther development in maize

miRNAs and their target genes are usually negatively correlated in expression levels[59].By combining RNA-Seq and sRNA-Seq analyses,the expression correlation coefficient was estimated for each miRNA-TF gene pair in W23 and the above four GMS mutant lines(ocl4,ms23,ms7-6007,andlob30)and another GMS mutant line(mac1)across their investigated stages.We identified 32 miRNA-TF gene pairs with negatively correlated expression in W23 line anthers(Pearson correlation,P＜0.05).Similarly,11,12,7,8 and 10 miRNA-TF gene pairs were detected inocl4,ms23,ms7-6007,lob30andmac1mutant lines,respectively(Pearson correlation,P＜0.05).In total,52 independent miRNA-TF gene pairs showed significant negative correlations in at least one investigated line of W23 and the five GMS mutant lines(Fig.S7;Table S8).The number of putative miRNA-TF gene regulatory pairs in W23 transcriptomes was larger than that in the five investigated GMS lines,mainly resulting from a larger sample size in W23 anther with 10 stages than that in the GMS line anthers with 3 to 5 stages.Notably,the previously reportedzma-miR172and its three target genes(ZmAP2,ZmSID1,andZmIDS1)were identified(Fig.S7,marked in red)[28,60],implying the remaining miRNATF gene pairs identified here may function in contributing to maize anther development.

Fig. 4. Expression analyses of stage-DE and GMS-DE miRNAs during maize anther development. (A) 118 stage-DE miRNAs in anthers of maize W23 line. Left, Z-score normalized expression tendencies displayed stage-DE (stage-differentially expressed) patterns; Middle, expression levels represented by CPM (read counts per million mapped reads) values; Right, differentially expressed miRNAs in ocl4, ms23, ms7-6007, and lob30 GMS mutants compared to their WT anthers, respectively. (B) Expression fold changes (FC, mutant/WT) of GMS-DE miRNAs between ocl4, ms23, ms7-6007, and lob30 mutant anthers and their corresponding WT line anthers. miRNAs with their names marked in red in Fig. 4B were absent from Fig. 4A.

To increase the detection sensitivity and identify more potential miRNA-TF gene pairs associated with anther development,we used an approach based on gene expression changes at each anther developmental stage between WT and the mutant line,and hypothesized that a miRNA-TF gene pair with down-regulated miRNA expression and up-regulated mRNA expression,or vice versa,in one GMS line at one sampled stage is an indicator of the real regulatory role.As a result,we found 56 miRNA-TF gene pairs showed opposite expression changes in more than half of the premeiotic/meiotic anther stages(S3 to S7)(i.e.more than seven of 15 stages in the three lines in total)inocl4,ms23andmac1lines,and 49 pairs showed opposite expression alterations in more than half of the meiotic/postmeiotic anther stages(S8 to S10)(i.e.more than three of seven stages in the two lines in total)inms7-6007 and lob30lines(Fig.5A;Table S9).The total 81 independent miRNATF gene pairs detected here were classified into three groups(Fig.5A),of which 24 pairs were detected in both premeiotic/meiotic and meiotic/postmeiotic transcriptomes(group I).When considering that miRNA-TF gene pairs with miRNAs in the same miRNA family and TF genes in the same TF family were classified into the same miRNA-TF module,the 81 miRNA-TF gene pairs were grouped into 25 miRNA-TF regulatory modules(Table S9).Among them,10 modules had 17 orthologous miRNA-TF gene pairs involved in anther development in other plant species(Fig.5B),includingmiR159-MYB,miR319-MYB,miR160-ARF,miR164-NAC,miR166-HD-ZIP,miR171-GRAS,miR172-AP2,miR319-TCP,andmiR396-GRF[29–34,36,37,60,61],as well asmiR390-ARFthat was mediated by trans-acting small interfering RNAs(tasiRNAs)[62].The miRNAs and their binding sites in TF transcripts were well matched(Fig.5C).In total,nearly a half(40%,10/25)of the detected miRNA-TF regulatory modules have been reported to function in contributing to plant anther development.Combining the 52 miRNA-TF gene pairs identified above by negative expression correlations(Fig.S7),a total of 99 independent miRNA-TF gene pairs were identified here,and might have higher probabilities to be involved in maize anther development.

3.6.ZmGAMYB,a maize TF homolog in miR319-and miR159-GAMYB modules,is essential for maize male fertility

In several plants,miR159-GAMYBandmiR139-GAMYBregulatory modules have been reported to play important roles in male reproductive process[34,35,63,64].The two miRNA-TF regulatory modules were also identified in this study(Fig.5A,B).In sRNA-Seq data of W23 anther,zma-miR319was highly expressed at premeiotic stages with the expression peak at stage 3 among the investigated stages(Fig.6A),whilezma-miR159was dominantly expressed during meiotic stages with the expression peaks at stages 7 to 9(Fig.6B).In contrary,in W23 anther RNA-Seq data,the expression ofZmGAMYB(Zm00001d012544)maintained at relatively high levels at stages 6 and 10(Fig.6A,B).The qRT-PCR analysis further confirmed the expression patterns ofzma-miR319,zma-miR159andZmGAMYBduring developmental stages 5 to 12 in WT anthers(Fig.6C).The negative correlations of expression levels betweenzma-miR319andZmGAMYB,as well as betweenzma-miR159andZmGAMYB,supported the combined regulatory roles of the two miRNAs onZmGAMYBat premeiotic and postmeiotic stages,respectively.By transcriptome analysis in GMS mutants,we found thatocl4mutant had little effect on the expression of the two miRNAs andZmGAMYB,ms23mutant lost the high expression ofmiR159at meiotic stages(S6 and S7),ms7-6007mutant resulted in down-regulated expression ofZmGAMYBat stage 10,andlob30mutant similarly repressedZmGAMYBexpression at stage 10 and increasedmiR319expression at stage 8b(Fig.6A,B).The gene expression alterations in GMS mutant lines were consistent with thatZmOCL4andZmMs23were preferentially expressed at premeiotic/meiotic stages,whileZmMs7andZmLOB30were mainly expressed during meiotic/postmeiotic stages,implying thatZmGAMYBmay be regulated by GMS TFs ZmMs7 and ZmLOB30 in addition to the two miRNA regulators(zma-miR319andzma-miR159).This inference was demonstrated by qRT-PCR analysis in whichZmGAMYBexpression was significantly down-regulated in bothlob30andms7-6007mutants at stage 10,respectively(Fig.6D).

ZmGAMYBexpression was associated with GMS TFs ZmMs7 and ZmLOB30 as well asmiR319andmiR159.Meanwhile,orthologous genes ofZmGAMYBin rice(OsGAMYB)andArabidopsis(AtMYB33andAtMYB65)were reported as GMS genes[63,65].We concluded thatZmGAMYBis important in maize anther development and male fertility.To verify this hypothesis,ZmGAMYBandZmGAMYB-2(Zm00001d043131,a paralog highly similar in sequence toZmGAMYB)double-gene loss-of-function mutant was created and obtained by targeted knockouts using the CRISPR/Cas9 system.The maizegamyb-1/2mutant displayed complete male sterility with shrunken and brown anthers and no pollen grain(Fig.6E,F,S2),demonstrating thatZmGAMYB,the predicted target gene ofzma-miR319andzma-miR159,andZmGAMYB-2are critical for maize anther development and male fertility.

4.Discussion

4.1.The expression pattern complexity of stage-differentially expressed TF genes and their regulatory roles in maize anther development

It is well known that tissue-or stage-specific expressed genes are important for the specific functions on tissue development and formation during the investigated developmental stages,including plant anther development[66].Among the reported GMS genes in plants,one third of them(31/101)encode TFs[4].Additionally,all the 17 cloned maize GMS genes show stage-DE manners during anther development[3,4,54,67–70].Therefore,stage-DE TF genes may play key regulatory roles in anther development.Here,by using a maize transcriptome dataset across ten anther developmental stages,we found that a half(1087/2216)of TF genes in maize genome showed stage-differential expression in maize anthers,and that the 1087 TF genes were mainly grouped into six expression clusters corresponding to the premeiotic,meiotic and postmeiotic development routes in maize anthers(Fig.1).Notably,majority of stage-DE TF genes had one expression peak during anther development,while there were hundreds of TF genes highly expressed with two or three peaks(Figs.1D,S1).Genes with high expression peaks at several discontinuous developmental stages may be involved in several biological functions or pathways during anther development.The different roles of the same gene at different developmental stages makes its functional complexity and should be deeply investigated to elucidate the molecular mechanism underlying anther development and plant male fertility in the future.

In functional enrichment analysis,we found that function annotations of stage-DE TF genes and their co-expressed or regulated genes well corresponded to the cytological and physiological observation results of maize anthers(Figs.1,2).For example,during anther development,archesporial cells and the four layers of anther wall were formed at premeiotic stages(S1–S5)(Fig.1A)[2],well corresponding to the enriched functions in biological processes including cell cycle checkpoint,cell proliferation,DNA replication,mitotic cell cycle and reproductive shoot system development in the premeiotic gene cluster(Figs.2A,S3).The cutin and wax layer and Ubisch bodies(Fig.2B)as well as tapetum degeneration(Fig.1A)were observed at stages 9 to 10,being consistent with the enriched functions in FA metabolism,sporopollenin biosynthesis and oxidation–reduction process for stage-DE genes in the meiotic cluster(Figs.2A,S3).At postmeiotic stages,pollen grains were filled with starch at stages 12 to 13(Figs.1A,2B),and the cutin and wax layer and Ubisch bodies were further concentrated and significantly enlarged at stages 12 to 13(Fig.2B),well corresponding to the enriched functions in active starch metabolism and transmembrane transport for genes in the postmeiotic cluster(Figs.2A,S3).From meiotic to postmeiotic stages,anther endothecium chloroplasts developed to a mature size and began to display photosynthesis activities as leaf chloroplasts since stage 10(Fig.2C)[1],which was also reflected by the enriched functions in photosynthesis and carbon fixation for genes in the middle-late cluster(Figs.2A,S3).These high associations between the anther cytological/physiological changes and the transcriptome alterations represented by stage-DE TF genes and their co-expressed genes during anther development,demonstrate that the stage-DE TF genes may play important roles in maize anther development,and that anther developmental processes may be precisely regulated by these stage-DEGs.

Fig.5.81 potential miRNA-TF gene pairs during maize anther development by expression alteration analysis in five GMS lines.(A)81 potential functional miRNA-TF gene pairs associated with maize anther development.Groups II,III and I indicate miRNA-TF gene pairs supported by premeiotic/meiotic GMS mutant lines(ocl4,mac1,and ms23),meiotic/postmeiotic GMS mutant lines(ms7-6007 and lob30),and both of them.(B)17 reported orthologous miRNA-TF gene pairs in maize and other plant species.(C)The miRNA binding sites on transcripts of target TF genes.

Fig.6.Expression and function analyses of ZmGAMYB gene and zma-miR159-ZmGAMYB and zma-miR319-ZmGAMYB regulatory modules.(A)Expression correlations between miR319 and ZmGAMYB in seven transcriptome data sets.(B)Expression correlations between zma-miR159 and ZmGAMYB in seven transcriptome data sets.(C)qRT-PCR analysis on expression patterns of zma-miR159,zma-miR319 and ZmGAMYB.(D)Expression pattern alterations of ZmGAMYB in ms7-6007 and lob30 mutant anthers in transcriptome data,respectively.(E)Gene knockout analysis of ZmGAMYB and ZmGAMYB-2 by using the CRISPR/Cas9 system.Phenotypic analyses of tassels,anthers and pollen grains stained with 1%I2-KI solution in WT and the double-gene knockout line gamyb-1/2.(F)Mutation site and fragment in the double-gene knockout line gamyb-1/2.

In addition, the function importance of stage-DE TF genes was also confirmed in the four maize GMS mutant (ocl4,ms23,ms7-6007, andlob30) anther transcriptome data containing a lot of DE TF genes, the majority of which were stage-specific expressed. As stage-DE GMS genes,ZmOCL4,ZmMs23,ZmMs7, andZmLOB30encoding TFs belong to HD-ZIP, MYB, PHD and LBD families,respectively, and play important roles in maize anther development and male fertility [5,7,9,10]. Their loss-of-function mutations influenced hundreds of GMS-DE TF genes (Fig. 3E, F) and their coexpressed genes, which are functionally enriched in biological processes and pathways closely related to anther development and male fertility (Fig. S6). Therefore, it can be speculated that multilayer and interwoven gene regulatory networks exist among stage-DE TF genes in anther development.

4.2.miRNAs contribute to maize anther development by regulating stage-differentially expressed TF genes

The molecular basis underlying complex traits in organisms is composed of key functional genes and their gene regulatory networks.Both transcriptional and post-transcriptional regulations are crucial to normal development and stress-responsive processes of cellular life.The miRNA-TF regulatory module is a well-studied post-transcriptional regulation pathway that also plays important roles in flower organ development and male fertility[3,58,71],such asmiR396-AtGRF5for the specification and formation of archesporial cells andmiR319-AtTCP24for secondary wall thickening in the anther endothecium cells inA.thaliana[32,33],miR164-OsCUC1for the establishment of the meristem boundary inO.sativa[36],andmiR171-SlGRAS24for tapetum ontogenesis and callose deposition around the tetrads inS.lycopersicum[37].Though several case studies and computational analyses have revealed many functionally important miRNA-target gene pairs regulating anther development in different plant species[38,72–74],the post-transcriptional regulation in anther development still remains elusive.

A comprehensive expression profile of miRNAs covering the anther whole development stages reported in our previous study[38]and new data supplemented here are a valuable data set for miRNA-TF regulatory analysis in maize anther development.The expression profile of miRNAs associated with anther development will elevate our understanding on the post-transcriptional regulatory mechanism of miRNAs to the GMS TF genes and other essential genes during anther development.In addition,comparative transcriptomics analysis has been commonly used to analyze and reveal the gene expression changes between WT and male-sterility mutants in plants,which can further identify functional genes involved in anther development[72–76].Here,by using RNA-Seq and sRNA-Seq data in the five GMS mutant lines(ocl4,ms23,ms7-6007,lob30,andmac1)and their corresponding WT lines,we identified 99 miRNA-TF gene pairs associated with maize anther development,and a large amount of them have been demonstrated to have important functions in male fertility in other plants(Figs.5A,S7).Asocl4,ms23,andmac1represent premeiotic/meiotic GMS lines andms7-6007andlob30represent meiotic/postmeiotic GMS mutants,the identified miRNA-TF gene pairs cover most of the anther whole developmental stages.Therefore,miRNA-TF regulatory modules identified here could provide many valuable candidates for further verification studies that will help us to better understand the regulatory roles of miRNAs and TF genes during anther development in maize and even other plant species.

4.3.The conserved functions of miR319-and miR159-GAMYB modules contribute to anther development and male fertility in plants

ZmGAMYBencodes a MYB TF and its orthologous genes in rice(OsGAMYB)andArabidopsis(AtMYB33andAtMYB65)were identified as plant GMS genes[34,35,63].In this study,ZmGAMYBwas grouped into the middle-late expression cluster(Table S4)and was demonstrated as a novel maize GMS gene by gene knockout assay(Fig.6E).Therefore,plantGAMYBgenes have relatively conserved functions in contributing to anther development and male fertility.Furthermore,asGAMYBwas preferentially regulated bymiR319in liverwort[77,78],mainly post-transitionally regulated bymiR159in rice andArabidopsis[34,63],and secondarily regulated bymiR319with a relative week effect compared to that ofmiR159inArabidopsis[61],it is most probable thatZmGAMYBis regulated by both or one of the two miRNAs in maize.By using transcriptome analysis and qRT-PCR assay,we foundZmGAMYBwas highly expressed at stages 5–6 and 10,while bothmiR319andmiR159were lowly expressed at these stages(Fig.6A–D).Correspondingly,whenZmGAMYBhad low expression levels before stage 5(premeiotic stages)and at stage 9,miR319andmiR159were highly expressed during the two periods,respectively(Fig.6A–D).In addition,bothzma-miR319andzma-miR159well matched binding sites onZmGAMYBtranscripts(Fig.5C).These results support thatZmGAMYBis regulated byzma-miR319andzma-miR159at premeiotic and meiotic/postmeiotic stages,respectively.Therefore,the functions ofmiR319-GAMYBandmiR159-GAMYBmodules controlling anther development are relatively conserved inArabidopsis,rice and maize,while the regulation patterns are likely different betweenmiR319-GAMYBandmiR159-GAMYBmodules in maize.

Additionally,ZmGAMYBexpression was down-regulated inms7-6007andlob30GMS mutant anthers at stage 10(Fig.6D),whilezma-miR319expression was up-regulated inlob30mutant anther at stage 8b but not largely changed inms7-6007mutant anther,and the expression patterns ofzma-miR159were also not largely altered in the two mutants(Fig.6A,B).Therefore,ZmGAMYBis likely regulated by both the two miRNAs(zma-miR319andzmamiR159)and the two GMS TFs(ZmMs7 and ZmLOB30).However,the detailed mechanism on how the two miRNA and two GMS TFs coordinately regulateZmGAMYBshould be further investigated.BesidesZmGAMYB,other stage-DE TF genes and miRNA-TF gene pairs identified in this study were candidates of gene knockout or overexpression analyses,by which their potential functions and regulation models in anther development can be further elucidated and verified.

CRediT authorship contribution statement

Ziwen Li:Methodology,Software,Formal analysis,Writingoriginal draft,Writing-review and editing.Taotao Zhu:Methodology,Investigation,Writing-original draft.Shuangshuang Liu:Formal analysis.Haoyun Liu:Formal analysis.Xueli An:Formal analysis,Investigation,Supervision.Yuwen Zhang:Investigation.Ke Xie:Investigation.Xiangyuan Wan:Writing-original draft,Writing-review and editing,project administration,Supervision.Jinping Li:Project administration.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We thank Virginia Walbot(Department of Biology,Stanford University)and Jixian Zhai(Department of Biology,Southern University of Science and Technology)for their help in the acquisition of transcriptome data.This research was funded by the National Natural Science Foundation of China(31771875,31971958,and 31871702),the Fundamental Research Funds for the Central Universities of China(2302019FRF-TP-19-013A1,06500136),and the National Key Research and Development Program of China(2017YFD0102001,2018YFD0100806,and 2017YFD0101201).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2021.03.016.