Huisheng Zhng,Huili Yng,Desheng Hu,Bing Li,Ynn Lin,Wen Yo,Zhnyong Guo,Hochun Li,Dong Ding,Zhnhui Zhng,Ynmin Hu,Ydong Xue,*,Jihu Tng,*

a National Key Laboratory of Wheat and Maize Crop Science,College of Agronomy,Henan Agricultural University,Zhengzhou 450002,Henan,China

b College of Life Sciences,Henan Agricultural University,Zhengzhou 450002,Henan,China

Keywords:

A B S T R A C T Normal microsporogenesis is determined by both nuclear and mitochondrial genes.In maize C-type cytoplasmic male sterility,it is unclear how the development of meiocytes and microspores is affected by the mitochondrial sterility gene and the nuclear restorer gene.In this study,we sequenced the transcriptomes of single meiocytes(tetrad stage)and early mononucleate microspores from sterile and restorer lines.The numbers of expressed genes varied in individual cells and fewer than half of the expressed genes were common to the same cell types.Four comparisons revealed 3379 differentially expressed genes(DEGs),with 277 putatively associated with mitochondria,226 encoding transcription factors,and 467 possibly targeted by RF4.KEGG analysis indicated that the DEGs in the two lines at the tetrad stage were involved predominantly in carbon metabolism and in amino acid biosynthesis and metabolism,whereas the DEGs during the transition from the tetrad stage to the early mononucleate stage were associated mostly with regulation of protein metabolism,fatty acid metabolism,and anatomical structure morphogenesis.Thus,meiocyte and microspore development was affected by the surrounding cells and the restorer gene,and the restorer gene helped restore the redox homeostasis of microspores and the normal cellular reconstruction during the transition.

1.Introduction

Heterosis,which underlies hybrid breeding,is a phenomenon in which multiple traits in the first filial generation are superior to those in the two parents or the means of the traits[1,2].Heterosis has been widely exploited for the crossbreeding of livestock and crops.However,large-scale application of heterosis in crops was hindered by practical difficulties associated with manual emasculation until male-sterile germplasm was introduced,in particular materials with cytoplasmically inherited male-sterility[3].The discovery and generation of plant materials showing cytoplasmic male sterility(CMS)further broadened the application of heterosis and promoted hybrid seed production in crops.All parts of CMS lines grow normally,with the exception of the anthers.

CMS genes,which are chimeric open reading frames(ORFs)located in the mitochondrial genome,have been identified and characterized in many plant species[4–6].The spatiotemporalspecific expression of CMS genes disrupts essential mitochondrial functions,including those involved in the tricarboxylic acid cycle and ATP generation,ultimately resulting in the production of nonfunctional pollen(conferring male sterility).The male fertility of CMS lines can be restored by restorer-of-fertility(Rf)genes via diverse mechanisms[6].In the CMS/Rfsystem,most clonedRfgenes belong to the pentatricopeptide repeat(PPR)family.The PPRRfgene products mitigate the detrimental effects of CMS genes via splicing,editing,cleavage,and/or polyadenylation at the posttranscriptional level.In rice,RF1A and RF5 are required to cleaveatp6-orf79dicistronic transcripts in CMS-BT andatp6-orfH79transcripts in CMS-HL,respectively[7,8].Some PPR family members restore the fertility of CMS lines at the translational or posttranslational level.The accumulation of WA352 protein is inhibited by RF3 in CMS-WA rice[9].

Besides the mRNA and protein levels,fertility restoration can occur at other levels.The first clonedRfgene,Rf2,restores the fertility of CMS-T plants at the metabolic level[10],whereasFrin the CMS-Sprite common bean mediates the excision of a specific CMSassociated mitochondrial genomic sequence.This latter example represents the first example of the fertility restoration of CMS lines at the genomic level[11].In most of these CMS cases,the levels of reactive oxygen species(ROS),including the superoxide anion(O2-),hydrogen peroxide(H2O2),and the hydroxyl radical(˙OH),are changed by the CMS factors[5].ROS can serve as signals to regulate plant development,and a breakdown in ROS balance leads to slowdown,pausing,or cessation of cell development,or even cell death[12,13].

The three major CMS types in maize have been designated as T,C,and S[14,15].The CMS-T and CMS-C types are sporophytic,whereas the CMS-S type is gametophytic.The expression of the mitochondrial geneT-urf13,which is toxic to tapetal cells,causes pollen abortion[16].In contrast,theRf2gene restores the fertility of CMS-T plants by detoxification[10].In CMS-S maize,the chimeric transcriptorf355-orf77affects microspore development[17]and is cleaved by RF3[18].The expression of mitochondrialorf355activates retrograde signaling to induce the expression of a nuclear gene,thereby further promotingorf355expression[19].The deleterious effects in CMS-C maize are alleviated by two duplicate genes[20,21].TheRf4gene,encoding a bHLH transcription factor(TF),is associated with a gain-of-function mutation in CMS-C maize[22],whereas the null allele is the underlying gene for male sterility(ms23)[23].

Because the proteins encoded by restorer genes may not be TFs that interact directly with CMS genes,there is substantial interest in the mechanism underlying the relationship between CMS-C andRfgenes.Several attempts involving bulk RNA sequencing(RNAseq)have been made to elucidate the complexity of the nuclearcytoplasmic interaction in CMS-C maize[24,25].The results clearly indicate that the protein encoded by the restorer gene does not directly interact with the CMS gene,but it remains unclear how theRf4gene product affects pollen development.

The ability to analyze the transcriptome of individual cells using single-cell RNA sequencing(scRNA-seq)technology has revolutionized investigations of cell development and differentiation.The profiling of RNA at the single-cell resolution has revealed the continuous developmental trajectory of highly heterogeneous cells in animals and plants,led to the discovery of new markers for specific cell types,and resolved intermediate or novel cell types[26–31].scRNAseq analysis of germinal cells from the archesporial cell stage to the prezygotene stage has clarified the premeiotic developmental programs of maize germinal cells[32].However,because of the presence of cell walls in plants,scRNA-seq technology has been used less extensively for plant research than for animal research.To date,sequenced plant samples have been limited to roots and early-stage germinal cells,from which protoplasts are readily isolated.

The two types of CMS systems,gametophytic and sporophytic,differ with respect to the anther cell part(sperm or somatic)affected by the CMS gene.In sporophytic CMS,the CMS gene directly interferes with tapetal cell development and indirectly induces microspore abortion.Previous functional characterizations and RNA-seq analyses focused on whole anthers,but the mechanism by which microspores are affected by the defective tapetal cells and the restorer gene remains unknown.

The objectives of this study were to investigate by scRNA-seq the transcriptomic changes in individual microspore cells to identify genes differentially expressed in the maize CMS-C and fertility restorer lines and to compare these changes with those in mitochondrial proteins expressed in the same developmental stages.

2.Materials and methods

2.1.Plant materials

The maize lines were identical to those used in our previous studies[25,33].The sterile lines(C,ES cytoplasm of CIII subgroup)and the restorer lines(R,ES cytoplasm of CIII subgroup)were grown at the experimental farm of Henan Agricultural University[Zhengzhou,Henan province,China(34°16′–34°58′N,112°42′–114°13′E)]from the spring to the summer of 2019.For each cultivar,more than 100 maize seedlings were grown to produce enough material for the experiments.Standard local field management practices(e.g.,irrigation,fertilization,weeding,and pest and disease control)were applied.

2.2.Microscopic observation and staging

The upper three anthers in a floret were collected;One anther was examined microscopically to determine the developmental stage and the other two were used for single-cell sampling(S8b and S9:tetrad stage and early mononucleate stage,respectively).Anthers were stained with acetocarmine and 4′,6-diamidino-2-phe nylindole to identify the precise anther developmental stage.Neighboring florets were then sampled and fixed in a formaldehyde:acetic acid:ethanol solution under vacuum.Fixed florets were dehydrated in an ethanol series and embedded in wax.Sections(4 μm thick)were stained and examined under an Olympus IX73 microscope(Olympus,Tokyo,Japan)at 20×magnification.

2.3.Tissue dissociation and cell isolation

Floret tissues were dissociated using a slightly modified version of the method described by Nelm and Walbot[32].Upper anthers at a known developmental stage within a floret were excised at 0.3–0.6 mm intervals and placed in a microcentrifuge tube containing 200 μL protoplast isolation buffer.Suspensions were incubated at 31 °C for 1 h with gentle inversion every 20 min.

Isolated protoplasts were gradually diluted by addition of 20 μL protoplasts to 2 mL 27% sorbitol solution until only a few protoplasts were detectable in a 10×visual field.A single protoplast was selected using a capillary glass tube(Fig.S1A and B),transferred to an 8-strip PCR tube(Axygen Low Profile 8-Strip PCR Tubes;ThermoFisher Scientific(China),Shanghai,CN)containing RNA lysis solution,and then flash-frozen in liquid nitrogen.Frozen cells were stored at-80°C.For each sample,16 cells per developmental stage were prepared in case of low-quality RNA isolation.

2.4.Single-cell RNA-seq library construction and sequencing

The reverse transcription of RNA and subsequent amplification,library construction,and sequencing were completed using the Smart-seq2 kit(Illumina,San Diego,CA,USA)as described by Picelli et al.[34].The library(300–700 bp fragments)was sequenced by the Wan Cheng Biotech Company(Shanghai,China)using the Illumina MiSeq platform to produce 150-bp paired-end reads.At least one gigabase of sequences was generated for each individual cell,corresponding to more than 3.3 million 150-basepaired-end reads.After mapping,most cells showed at least 20%sequence coverage of the maize genome.

Clean reads were obtained from the raw data for each cell by filtering out adapter sequences and the low-quality sequences(quality score＜30,length＜50 bp,and N frequency greater than 10%)using SOAPnuke software(version 2.1.0)(https://github.com/BGIflexlab/SOAPnuke;BGI,Shengzhen,China).Sequences were mapped to the maize reference genome B73v4 and reads were counted using HISAT2(version 2.2.0)(http://daehwankimlab.github.io/hisat2/)and StringTie(version 2.1.3)(http://www.ccb.jhu.edu/software/stringtie/),respectively.Differentially expressed genes were identified using the DESeq2 package(version1.20.0)(https://github.com/mikelove/DESeq2).A principal component analysis was performed using the top 500 genes ranked in terms of read count variance.For DEG analysis,each cell was treated as one biological replicate.The DEGs were defined as genes with a false discovery rate≤0.01 and|log2(Fold change)|≥1.Gene Ontology(GO)functional classification was performed using the AgriGO online toolkit(version 2.0)with the GO dataset(version 2016–8)(http://systemsbiology.cau.edu.cn/agriGOv2/).The Kyoto Encyclopedia of Genes and Genomes(KEGG)database(release 95.2)(http://www.genome.jp/kegg/)was used to identify the enriched pathways among the DEGs.The Benjamini-Hochberg method[35]was used to control the false discovery rate.KEGG pathways withP＜0.05 were assigned as significantly enriched.

2.5.RNA isolation and quantitative real-time(qRT)PCR

Anthers in S8b were sampled and staged as described in section 2.2.Total RNA was isolated,cDNA was synthesized,and qRT-PCR analysis was completed following the corresponding manufacturer’s manuals as previously described[33].Three independent biological and three technical replicates were used for qPCR.Primers for detecting DEGs were designed using Primer3Plus software(http://primer3plus.com).MaizeActin7(Zm00001d032480)was selected as an internal control gene.Relative gene expression levels were calculated using the 2-ΔΔCTalgorithm[36].Details of primers used for amplifying DEGs and the control gene are listed in Table S1.

3.Results

3.1.Cytological differences during anther development

The sterile lines completely lacked mature pollen grains,whereas the restorer lines produced viable pollen grains that were normal in size.Microspore development was normal in both lines until late stage 9(Fig.1A).Although microspores remained intact at S8a and early S9,we previously determined that changes in gene expression caused by the CMS gene were first detectable at metaphase I and peaked at S8b[25].By transmission electron microscopy(TEM),the tapetal cells of the restorer lines showed clear cell structures with many small vacuoles at S8b and began to break down at S9(Fig.S1C b and d).In contrast,many large vacuoles formed in the sterile line at S8b and the tapetal cells began to lose defined organelles(Fig.S1C a and c).Generally,it is presumed that this change occurs mainly in the tapetal cells in sporophytic CMS,but it is unclear how genes in the microspores are influenced and how the restorer geneRf4rescues pollen development.Accordingly,single tetrads(S8b)and early mononucleate(S9)microspores were selected using a capillary glass tube and transferred to a PCR tube.A total of 16 cells were selected for each stage.The transcriptome of 12 cells with high-quality RNA was profiled at each stage by scRNA-seq(Fig.1B;Table S2).

Fig.1.Cytological comparison of anther development and schematic diagram of single-cell RNA sequencing.(A)Cytological comparison of anther development in CR87-1(a–e)and Cr87-1(f–j).S8a,dyad stage;S8b,tetrad stage;S9,mononuclear stage;S10,late mononucleate stage;S13,mature pollen stage;Dy,dyad cell;E,epidermis;En,endothecium;M,middle layer;MP,mature pollen;Msp,microspores;T,tapetal layer;Tds,tetrads.Scale bars,50 μm.(B)The process of single-cell RNA sequencing,including single-cell isolation,dispersing cells,single-cell RNA extraction,and sequencing.

3.2.Individual cell expression and cell clustering analysis

The sequencing of 48 libraries(one library per cell and 12 replicates per stage)generated＞2.05 billion paired-end reads.After removal of adapter and low-quality sequences,the remaining＞170 million high-quality reads were aligned to the genome(Table S2).Expressed genes were defined as yielding more than one fragment per kilobase per million reads in three more cells.The number of expressed genes per cell was determined(Fig.2A).Cells varied with respect to numbers of expressed genes,but the number of commonly expressed genes was higher between cells in the same group than between cells in different groups.The mean numbers of expressed genes in RQ and rQ(R,restorer line;r,sterile line;Q,S8b)were 12,588 and 12,121,respectively.The number of expressed genes decreased to 3561 in RM(M,early mononucleate stage)and 2132 in rM.Only 132 genes were expressed in all cells.In contrast to tagged pool single-cell sequencing,individual cells were at specific developmental stages.To test if the selected cells were uniform and at the specified stages,a principal component analysis was performed to cluster cells.As expected,the 48 cells were grouped into four major clusters(stages×materials),and the cells of each plant material at a specific stage were clustered(Fig.2B).Tetrad cells from two plant materials were more clearly separated than early mononucleate cells.

3.3.Major patterns of gene activity and associated pathways

Fig.2.Gene counting and PCA analysis for cells.(A)The number of expressed genes for each cell.Twelve cells were sequenced for each stage by genetic group.Blue bar:commonly expressed genes in the same group;Red bar:distinct expressed genes among individual cells in the groups;Horizontal black line:genes commonly expressed among the four groups.(B)Principal component analysis of expressed genes across all cells.Each dot represents a single cell.R,restorer line;r,sterile line;Q,S8b;M,early S9.

Fig.3.Major expression patterns exhibited by all expressed genes and enriched pathway sub-categories.Grouping of expressed genes was performed based on K-means clustering and pathway enrichment was tested with the hypergeometric test(P≤0.05).The top 10 pathways are shown for groups containing more than 10 pathways.All pathway sub-categories were significantly enriched after multiple correction with adjusted P-value(FDR＜0.05).The Y-axis of the left chart represents the relative expression level of each gene after normalization in each cell.The X-axis of the right chart shows the enrichment factor.R,restorer line;r,sterile line;Q,S8b;M,early S9.

To identify co-regulated gene clusters,K-means clustering was used to reveal major expression patterns associated with developmental progression.Six clusters(clusters 1–6)were extracted from all expressed genes(Fig.3;Table S3).Cluster 1 comprised 337 genes that were highly expressed in the tetrads of the restorer lines.These genes were associated with regulation of cell wall structure,cell redox state,and post-translational modifications.A total of 22 pathway sub-categories were represented in this cluster,with carbohydrate derivative metabolic process as the most common pathway,followed by organophosphate metabolic process,nucleotide metabolic process,and ribose phosphate metabolic process.Cluster 2 contained the most genes,which were involved in 75 pathways involved mainly in protein translation,proton transfer,and ATP metabolism.The expression levels of the genes in this cluster were high in tetrads and considerably lower in mononuclear microscopes in the restorer lines.In contrast,the expression of most of the genes in this cluster was misregulated at the same stages in the sterile lines(Fig.3).Cluster 3 included 445 genes that were expressed mainly the tetrads of the restorer lines.Some genes were also expressed at the mononuclear stage.The expression of these genes was almost undetectable in the tetrads of the sterile lines.A few genes contributed to the following pathways regulating protein degradation:ubiquitindependent protein catabolic process,proteolysis involved in cellular protein catabolic process,and cellular macromolecule catabolic process.Cluster 4 comprised 394 genes that regulate microspore development,including genes involved in protein localization and hormone synthesis.Clusters 5 and 6 consisted of 445 and 270 genes,respectively,but no enriched pathways were detected for these genes.Six gene clusters that showed cell and stagespecific expression were also identified(Fig.S2).Cluster 1 comprised 72 genes expressed specifically in rM cells,whereas cluster 2 comprised 28 genes expressed specifically in RM cells.The three cluster 3 genes were expressed specifically in mononuclear microspores.The 327 genes in cluster 4 were expressed exclusively in rQ cells.Clusters 1–4 lacked enriched pathways.Cluster 5 contained 839 genes that were expressed only in RQ cells.Many of these genes were involved in proteolysis,the small molecular biosynthesis process,and the ribose phosphate metabolic process.Cluster 6 consisted of 952 genes that were expressed in tetrad cells and were associated predominantly with proton transmembrane transport,the purine ribonucleoside triphosphate metabolic process,and translation.

3.4.Deduced functions of DEGs and validation by qRT-PCR

To investigate genes that were differentially expressed in meiocytes and microspores,individual cells that deviated from the main clusters(RM3;RQ1,3,4,and 5)(Fig.2B)or showed abnormal numbers of expressed genes(rQ9,10,11,and 12)(Fig.2A)were excluded from subsequent DEG analysis.DESeq2 was used to search for DEGs between the developmental stages and between the sterile and restorer lines.The following threshold was used to determine the significance of gene expression level differences:adjustedP≤0.01 and|log2Fold change|≥1.A total of 3,379 significant DEGs were identified in the four comparisons.Of these genes,2,113(974 up-regulated and 1,139 down-regulated)were detected in the rQ_vs._RQ comparison,whereas respectively 599(290/309),758(400/358),and 195(126/69)were detected in the rM_vs._RM,RM_vs._RQ,and rM_vs._rQ comparisons(Fig.4A and B).

Fig.4.Differentially expressed gene(DEG)analysis and validation of DEGs by qRT-PCR.(A)Venn diagram indicating the overlap of genes identified in comparisons.(B)Upregulated and down-regulated genes identified in comparisons.(C)KEGG enrichment analysis in the comparison between the sterile and restorer lines.(D)Relative expression of 18 genes measured by FPKMs in single cell RNA-seq results.(E)Relative expression of 18 genes measured by qRT-PCR.R,restorer line;r,sterile line;Q,S8b;M,early S9.

The functions and pathways that varied between the sterile and restorer lines and from S8b to early S9 were identified as follow.The DEGs detected in the rQ_vs._RQ,rM_vs._RM,RM_vs._RQ,and rM_vs._rQ comparisons were annotated(P＜0.01)with 539,188,0,and 13 GO terms,respectively.In the rQ_vs._RQ comparison,471,3,and 65 GO terms assigned to DEGs were from the biological process(BP),molecular function(MF),and cellular component(CC)categories,respectively(Table S4).For the rM_vs._RM comparison,the DEGs were annotated with 160,26,and 2 GO terms from the BP,CC,and MF categories,respectively.The DEGs revealed by the rM_vs._rQ comparisons were annotated with 13 GO terms from the CC category.For the rQ_vs._RQ comparison,‘cytoplasm,’‘cytoplasmic part,’and‘plasma membrane’were the most significantly overrepresented GO terms in the CC category;terms associated with responses to diverse stresses,such as metal ion and abiotic stimuli,were the most represented GO terms in the BP category;and‘pigment binding’,‘GTPase activity’,and‘copper ion binding’were the only significantly enriched GO terms in the MF category.In the rM_vs._RM comparison,‘electron carrier activity’and‘oxidoreductase activity’were the only overrepresented GO terms in the MF category;the most enriched GO terms in the BP category were almost identical to those revealed in the rQ_vs._RQ comparison;and‘a(chǎn)poplast’was the most enriched GO term in the CC category,followed by‘extracellular region’and‘cytoplasmic part’.For the rM_vs._rQ comparison,‘cytoplasmic part’,‘membranebounded organelle’,and‘intracellular membrane-bound organelle’were the most enriched GO terms in the CC category.

The DEGs detected in the four comparisons were assigned to 30 significant(P＜0.05)KEGG pathways,of which‘carbon metabolism’(n=48),‘biosynthesis of amino acids’(n=45),‘plantpathogen interaction’(n=34),‘a(chǎn)mino sugar and nucleotide sugar metabolism’(n=30),and‘MAPK signaling pathway plant’(n=27)were highly enriched(Fig.4A).The KEGG analysis also revealed that‘glutathione metabolism’,‘citrate cycle(TCA cycle)’,‘glycolysis/gluconeogenesis’,and‘fatty acid biosynthesis’were pathways enriched in more than one comparison(Fig.4C;Table S4).

Eighteen DEGs were randomly selected to validate the transcriptome data in S8b by qRT-PCR.The selected DEGs were involved predominantly in transcription,the cell cycle,RNA processing,cell wall modeling,transport,pollen development,and other pathways.The observed changes of 11 genes in relative expression levels were consistent between the scRNA-seq and qRT-PCR results,reflecting the accuracy of the scRNA-seq analysis(Fig.4D and E).However,seven genes displayed trends differing from those of RNA-seq because of the inclusion of the outer fourlayer tissues around meiocytes for RNA isolation.Some genes showed expression patterns differing between meiocytes and the tapetal layer.

3.5.Pollen development-associated TFs

Changes in TF abundances will lead to substantial changes in gene expression during the transition from tetrads to mononucleate microspores and determine subsequent development.We examined the types of TFs among the DEGs.In normally developing microspores(in the restorer lines),the expression levels of 26 and 22 TF-encoding DEGs were lower and higher,respectively,at the early S9 than at S8b(Fig.5A and B).In contrast,only seven down-regulated and four up-regulated TFencoding DEGs were detected in the same transition in the sterile lines.Only one TF gene(Zm00001d041395)was common in the restorer and sterile lines in this transition process.We also investigated the differentially expressed TF genes at the same stages in these two materials.In total,152 TF genes were detected at the S8b,of which 101 were down-regulated and 51 were up-regulated.At the early S9,4 and 22 TF genes were down-regulated and up-regulated,respectively,in the comparison between the microspores of the sterile and restorer lines.These TF-encoding DEGs included members of the HSF,MADS,MYB,NAC,ARF,bZIP,and other TF families.Most of the TF genes from these families,including the restorer geneRf4,which encodes a bHLH transcription factor[22],were down-regulated during the transition from tetrads to mononucleate microspores.These genes also showed lower expression in the sterile lines than in the restorer lines at S8b.

3.6.Cell phase transition-associated DEGs

From tetrad meiocytes to mononucleate microspores,cells undergo a series of morphological and internal structural changes that prepare them for the subsequent mitosis.We examined the expression changes of 456 genes involved in the process of cell phase transition[32].We identified 16 DEGs associated with cell cycles,with eight up-regulated and eight down-regulated genes in the restorer and normal lines.Most of the up-regulated DEGs were expressed in the G2 phase of the cell cycle,whereas most of the down-regulated genes were expressed in the G1 and S phases(Table 1).These results indicated that the selected microspores had finished synthesizing DNA and were preparing for mitosis.However,at S8b in the sterile lines,compared with the restorer lines,most of the G2 phase gene expression levels were down-regulated,unlike the up-regulated expression of the G1 and S phase genes,reflecting the delayed development of the meiocytes in the sterile lines(Table S5).

Table 1Cell cycle-related DEGs in cell phase transition.

Table 2Number of DEGs targeted to mitochondria.

Fig.5.Types and numbers of transcription factors in DEGs and WGCNA results for the Rf4-co-expressed gene.(A)Numbers of up-and down-regulated transcription factors in the different comparisons.(B)Distribution of the types of transcription factors.(C)The network of the co-expressed gene with the Rf4 gene.(D)The enriched pathways of genes in the Rf4-regulated network.R,restorer line;r,sterile line;Q,S8b;M,early S9;rQ_vs._RQ,comparison between sterile and restorer lines at S8b;WGCNA,weighted gene co expression network analysis.

3.7.Mitochondrial function-related DEGs

Genes involved in mitochondrial activities are vital for cell function and development.In a previous study[33],we identified 1840 mitochondrial proteins in the anthers of sterile and restorer lines at S8b using iTRAQ technology.In the present study,we determined whether the corresponding genes for these proteins were differentially expressed.A total of 219,45,36,and 15 DEGs in the rQ_vs._RQ,rM_vs._RM,RM_vs._RQ,and rM_vs._rQ comparisons,respectively,were detected in this mitochondrial gene set.In the normal/restorer lines,the expression levels of most DEGs associated with mitochondria were up-regulated(29 up-regulated vs.7 down-regulated)(Tables 2,S6).The up-regulated DEGs were involved in lipid metabolism and energy metabolism pathways.At S8b,the expression levels of most of the mitochondrial DEGs in the sterile lines were up-regulated(159 up-regulated vs.60 down-regulated),implying that the internal homeostasis in mitochondria had been disrupted by the defective tapetal cells.The expression levels of two DEGs encoding cytochrome P450 were 10 and 35 times higher in the normal/restorer lines than in the sterile lines.One of these genes was thems10male sterility gene.

3.8.Differential expression of the Rf4 gene and its targets

The allele of theRf4restorer gene in CMS-C maize is also annotated asMs23.TheRf4/Ms23gene is preferentially expressed in the tapetal cells,with peak expression at stage 10[22,23].An analysis of the microspores revealed the differential expression pattern of this gene.TheRf4(Ms23)expression level was significantly higher at early S9 than at S8b,indicating thatRf4expression was downregulated earlier than in tapetal cells(Table 3).Moreover,rf4expression at S8b was 3.6-times higher in the sterile lines thanRf4in the restorer lines.TheRf4gene belongs to the bHLH TF family,whose members typically bind to E-box sequences(CANNTG),such as the G-box(CACGTG),in the promoter region of target genes[37].We accordingly searched the promoter region of DEGs for Ebox sequences(at least two binding sequences)to identify the RF4-targeted genes,whose expression profiles were then determined.Sets of 346 DEGs from the rQ_vs._RQ comparison and 133 DEGs from the RM_vs._RQ comparison contained E-box sequences in their promoter regions(Table S7).After combining the non-DEGs,we used RF4 as the central controller to construct a regulatory network with only one layer.The transcriptome network contained 120 genes with three network clusters(Fig.5C).Functional enrichment analysis revealed that the network genes were associated primarily with protein-DNA complexes,chromatin,and DNA packaging complexes(Fig.5D).The RF4 protein is not localized to mitochondria,but one of its target genes must be imported into mitochondria to restore the fertility of CMS-C maize.Thus,we further narrowed the set of secondary restorer candidates to nine genes in the rQ_vs._RQ comparison(Table 3).Specifically,

Table 3Expression of Rf4 and its nine putative target genes.

Zm00001d010805andZm00001d030053encode proteins with a DNA-binding motif,suggesting that they may regulate the expression of mitochondrial genes.TheZm00001d051442andZm00001d039505genes encode a methyltransferase and mitochondrial-processing peptidase subunit alpha-2,respectively,and are probably involved in protein processing.In contrast,Zm00001d024595encodes a monocopper oxidase-like protein that may participate in cell wall expansion,whereasZm00001d048198encodes a yippee-like protein and is required for normal cell proliferation.The protein encoded byZm00001d043660regulates intracellular auxin homeostasis to promote pollen development.

Fig.6.The putative regulation pathway of the Rf4 gene on the tapetum layer and in microspores.The targeted genes were selected based on co-expression analysis and the recognition motif of bHLH-type transcription factor.S8b,tetrad stage;S9,mononucleate stage.

4.Discussion

The restorer gene for CMS-C maize,Rf4,encodes a bHLH-like TF[22]and its null allele is responsible for male sterility in maize[23].TheMs23gene specifies early tapetal layer development,and in thems23mutant,meiocytes are suspended at prophase I.However,in CMS-C maize,tapetal cells initially develop normally and meiocytes progress to the mononuclear stage(Fig.1A).The expression ofRf4ensures that meiocytes continue to develop into mature pollen grains,which has different effect with thems23allele.Based on transcriptome data from scRNA-seq,the 48 cells were classified into four groups with no marked difference between RM and rM.This is probably because genes related to meiosis were not expressed when cells entered microspores and few genes of microspores are affected by the restorer gene(Fig.4A and B).We also investigated howRf4contributes to meiotic processes.Our findings suggest that RF4 may regulate the expression of several male sterility-associated genes until S8b to ensure that the tapetum layer functions normally and microspores continue to develop(Fig.6).Compared with the RNA-seq data from the whole-anther,the majority of DEGs showed different expression trends,indicating that the regulation of expression of genes in the whole anther(tapetum)differed from that in meiocytes and microspores.

4.1.A master regulator encoded by the restorer gene influences meiocyte transcriptome at tetrad stage

Earlier studies of anthers confirmed that gene expression levels change only slightly at the pollen mother cell stage and during early meiosis,but subsequently undergo substantial changes at S8b and S9[25,38].These gene expression changes occur in the outer anther layer and in the inner meiocytes or microspores,possibly interfering with investigation of the mechanisms associated with sterility or fertility restoration.However,sequencing involving homogeneous cells can eliminate the problems associated with analyses of heterogeneous cells,and the scRNA-seq method is superior to the bulk RNA-seq method for detecting DEGs.Using a 2-fold change in expression as the threshold,we identified the most DEGs(2113)in the comparison rQ_vs._RQ(Fig.4A).Of these DEGs,152 encoded TFs(Fig.5A).Given that near-isogenic lines were studied,these gene expression changes were likely associated with a single fertility restorer gene.Absence of expression of this gene reportedly[23]alters the expression of several bHLH genes at the pre-meiosis stage.The downstream genes were involved mainly in energy-related pathways,including carbon metabolism and amino acid biosynthesis.In CMS-C maize,the expression levels of genes associated with mitochondrial energy metabolism and protein maturation are up regulated in the anthers of restorer lines[33].In the present study,the expression of DEGs involved in mitochondrial protein folding,ATP synthase,and mitochondrial structure was down-regulated in CMS-C lines,whereas the expression of DEGs encoding long chain acyl-CoA synthetase1,receptor-like protein kinase 1,and ascorbate peroxidase 2 was up-regulated.On the basis of the nature of the C cytoplasm in our maize materials,we deduced that the restorer gene rescued microspore development by promoting the expression of genes involved in mitochondrial functions.

4.2.Cell-cycle transition is probably remodeled by the restorer gene

Meiosis-related cell-cycle processes,especially early meiotic events,are critical to sexual reproduction.The sequencing of the maize transcriptomes of clustered meiocytes at specific developmental stages has resulted in the identification of many cell cycle-associated genes homologous toArabidopsisgenes[39,40].Further transcriptome analyses of multiple individual meiocytes confirmed the specificity of the expression of these cell-cyclephase genes,and identified marker genes for the entry into meiosis[32].However,the transition,exit from meiosis,and re-entry into mitosis have not been thoroughly investigated.Comparison of microspore development between S8b and early S9 showed different expressions of many cell-cycle-associated genes,suggesting that these changes were in preparation for entry into mitosis(Tables 1,S5).The expression of ribosomal genes and membraneassociated genes also changed significantly,possibly indicative of the remodeling of the cytoplasm from meiosis to mitosis.This expression switch and reprogramming has been observed in the early stages of the entry into meiosis in maize[32],in microspores at the five major stages in rice[41]and in the anther/tetrad/microspore transcriptomes of other plant species[42,43].Thus,in this late developmental stage,the expression ofRf4(Ms23)may ensure that remodeling occurs by regulating downstream cell-cycle genes.

4.3.TFs function synergistically with the restorer gene

Changes in TF gene expression or the structure of the encoded proteins can markedly affect plant developmental transformations[44].Many TFs mediate anther development and pollen maturation.In an earlier study in sorghum,five of seven DEG groups were enriched with TF genes,most of which were members of the B3,SBP,NAC,and MYB families[45].A mutation toOsMYB80causes premature tapetal cell death,loss of Ubisch bodies,and microspore degeneration[46].TheSPL8andBIM1genes,the latter of which encodes a bHLH TF,function redundantly to control pollen sac development and male fertility inArabidopsis[47,48].Transcriptome analyses of anthers from other CMS plants,including tobacco[49]and pigeon pea[50],indicated that MYB,bHLH,NAC,and other TF families are involved in the sterility and fertility restoration of CMS lines.In the present study,225 TF genes were included among the DEGs.Although these TF genes represented many TF families,the most common genes were from the MYB,bZIP,NAC,EREB,and GRA families(Fig.5B).Nine DEGs,including theRf4/Ms23gene,which encodes a bHLH TF,were detected in the tetrads and microspores of CMS-C maize plants.The rice bHLH protein,TIP2,is critical for meristemoid transition and differentiation during the early anther developmental stages[51].InArabidopsis,DYT1,a bHLH TF,regulates tapetum function and pollen development[52].Thus,TFs,including bHLH family members,are crucial for normal tapetum development and meiosis.Several bHLH TFs may function cooperatively in a sequential manner or as complexes to coordinate transcription in developing anthers[53].For example,theRf4gene product interacts with other bHLH TFs(e.g.,bHLH51)[23].Furthermore,by searching for the recognition motif in DEG promoters,two of nine bHLH-type DEGs were identified as putative target genes of RF4.However,the rQ_vs._RQ comparison indicated the expression levels of these two bHLH TFencoding genes andRf4showed opposite trends.Thus,the functions of the proteins encoded by these genes may differ between the tetrads and the tapetum.

4.4.Redox homeostasis in CMS-C maize is recovered by the restorer gene

Redox systems in mitochondria are central hubs in plants that affect cellular metabolism,signaling,plant development,and stress responses[13].CMS factors cause oxidative stress in mitochondrial inner membrane space and matrix,inducing ROS accumulation or burst[5].In contrast in the restored lines,the ROS level was reversed to the former status or reached a new balanced level.The glutathione metabolism pathway was enriched multiple times in the three-comparison analyses(Table S4),indicating that the ROS status was changed significantly.In the putative mitochondria-targeted DEGs,many genes were annotated as maintaining redox statein vivo,such as thiol protease,thioredoxin,and superoxide dismutase(Table S6).These results show that the restorer gene re-balanced the ROS state in maize CMS-C.

CRediT authorship contribution statement

Dong Ding,Yadong Xue,and Jihua Tang:Conceptualization.Huaisheng Zhang,Huili Yang,Haochuan Li,Zhanhui Zhang,and Yanmin Hu:Data curation.Wen Yao,Desheng Hu,and Zhanyong Guo:Formal analysis.Dong Ding,Yanmin Hu,and Jihua Tang:Investigation.Huaisheng Zhang,Huili Yang,Bing Li,and Yanan Lin:Funding acquisition.Huili Yang and Yadong Xue:Methodology.Dong Ding and Yadong Xue:Project administration.Yanmin Hu and Jihua Tang:Resources.Yadong Xue:Writing–original draft.Haochuan Li,Dong Ding,Zhanhui Zhang,and Yadong Xue:Writing–review&editing.All authors have read and agreed to the published version of the manuscript.

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

This work was supported by the National Natural Science Foundation of China(31571745 and 31971893),the Key Technology Research and Development Program of Henan Province(202102110164 and 212102110061),the Zhengzhou Major Science and Technology Innovation Project(188PCXZX803),and the Open Funds of the State Key Laboratory of Crop Genetics and Germplasm Enhancement(ZW202001).

Appendix A.Supplementary data

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