Impact of Autophagy on Gene Expression and Tapetal Programmed Cell Death During Pollen Development in Rice

Autophagy has recently been shown to be required for tapetal programmed cell death (PCD) and pollen maturation in rice. A transcriptional regulatory network is also known to play a key role in the progression of tapetal PCD. However, the relationship between the gene regulatory network and autophagy in rice anther development is mostly unknown. Here, we comprehensively analyzed the effect of autophagy disruption on gene expression profile during the tapetal PCD in rice anther development using high-throughput RNA sequencing. Expression of thousands of genes, including specific transcription factors and several proteases required for tapetal degradation, fluctuated synchronously at specific stages during tapetal PCD progression in the wild-type anthers, while this fluctuation showed significant delay in the autophagy-deficient mutant Osatg7-1. Moreover, gene ontology enrichment analysis in combination with self-organizing map clustering as well as pathway analysis revealed that the expression patterns of a variety of organelle-related genes as well as genes involved in carbohydrate/lipid metabolism were affected in the Osatg7-1 mutant during pollen maturation. These results suggest that autophagy is required for proper regulation of gene expression and quality control of organelles and timely progression of tapetal PCD during rice pollen development.


INTRODUCTION
Reproductive development, both in animals and plants, requires drastic metabolic changes for nutrient supply to gametes. In flowering plants, anthers exhibit a four-layered structure composed of epidermis, endothecium, middle layer, and tapetum. The tapetum, the innermost of the four sporophytic layers of the anther wall, directly contacts with the developing microspores and provides metabolites and nutrients to pollen grains and pollen coat during their development (Ariizumi and Toriyama, 2011). Therefore, proper metabolic regulation is critical for proper pollen maturation and grain yield.
The tapetum acts as a nutritional source for the developing microspores by undergoing degeneration triggered by programmed cell death (PCD) from stage 7 (ST7) to ST11 (Ariizumi and Toriyama, 2011;Zhang and Yang, 2014). After two rounds of cell division during meiosis, tetrads covered by callose are formed at ST8. After callose degradation, at ST9, haploid microspores are released into the lobe, and the accumulation of starch and lipids are associated with the pollen maturation and pollen wall development (Zhang et al., 2008;Zhang et al., 2011).
Transcriptional regulatory network plays key roles in the progression of anther/pollen development in plants. Transcriptomic analyses of developing anthers and pollen grains have been described in many species (Huang et al., 2011;Rutley and Twell, 2015) and enriched our knowledge of the repertoire of genes expressed during anther/pollen development. Transcriptomic profiling at different developmental stages has been reported for Arabidopsis (Honys and Twell, 2004), rice (Wei et al., 2010), and tobacco (Bokvaj et al., 2015).
Rice is not only a widely-grown crop to feed almost half of the global population but also an ideal monocot model plant facilitating plant science research, e.g. the discovery of mechanisms underlying male gametophyte development (Kim and Zhang, 2018), and thus is of significance for both basic and applied plant research. Several bHLH-TFs including undeveloped tapetum1 (UDT1) (Jung et al., 2005), tapetum degeneration retardation (TDR) (Li et al., 2006), TDR-interacting protein2 (TIP2) (Fu et al., 2014;Ko et al., 2014), and eternal tapetum1 (EAT1) (Niu et al., 2013;Ono et al., 2018), a MYB transcription factor GAMYB (Aya et al., 2009), and a plant homeodomain (PHD)-finger motif DNA binding protein persistent tapetal cell1 (PTC1)  have been reported to be required for tapetal PCD and identified as key genes involved in a transcriptional network during anther development as well as tapetal PCD progression. Comparative studies with Arabidopsis developing anthers including tapetum have revealed complexed and different expression pattern of genes during reproductive periods (Li et al., 2006;Deveshwar et al., 2011).
Autophagy is a major intracellular mechanism that degrades organelles, proteins, and metabolites (Mizushima and Komatsu, 2011). The double membrane structure called autophagosomes is formed in the cytosol, which fuses with the lytic compartments, the vacuoles or lysosomes, where the components enclosed in the autophagosomes are degraded (Kurusu et al., 2016). More than 30 autophagy-related genes (ATGs) required for autophagy are well conserved basically in all eukaryotes including animals and plants (Avin-Wittenberg et al., 2012;Yoshimoto and Ohsumi, 2018).
Autophagy plays essential roles in growth, development, and survival in eukaryotic cells (Chung et al., 2009;Li and Vierstra, 2012;Yoshimoto, 2012). Under normal growth conditions, autophagy plays critical roles in nutrient recycling, cell waste management and quality control of organelles as a housekeeping mechanism. When cells are faced with nutrient deprivation or stress, remobilization of nutrients by autophagy becomes crucial (Thompson et al., 2005;Shibata et al., 2013;Yoshimoto et al., 2014).
In plants, autophagy has been suggested to play roles in the recycling of proteins and metabolites, including lipids, at the whole-plant level, and is involved in numerous physiological processes (Ghiglione et al., 2008;Ishida et al., 2008;Chung et al., 2009;Yoshimoto et al., 2009;Kwon et al., 2010;Izumi et al., 2015;Izumi et al., 2017). Transcriptomics and metabolomics in autophagy-deficient mutants of Arabidopsis (atatg5-1) and maize (Zmatg12) suggest the importance of autophagy in cell homeostasis and stress responses. These multi-omics also provides comprehensive data sets for the identification of proteins, protein complexes, organelles and processes directly or indirectly under autophagic control.
In rice, autophagy contributes to degradation of chloroplasts/ plastids  and efficient nitrogen remobilization and biomass production at the whole-plant level by facilitating protein degradation for nitrogen recycling  as well as starch metabolism during seed development (Sera et al., 2019). Rice mutants defective in autophagy, Osatg7-1, Osatg7-2, and Osatg9, show sporophytic severe male sterility under normal growth conditions . Autophagy is involved in phytohormone and lipid metabolism in anthers and is crucial for sexual reproductive development Kurusu and Kuchitsu, 2017). Pollens from the Osatg7-1 mutant are premature due to significant defects in anthers during pollen maturation , and autophagy is induced inside the tapetum at the uninucleate stage (ST9-10) . Moreover, transmission electron microscopy (TEM) analysis has shown that the morphology of tapetum at ST8-9 is normal, but the tapetal collapse is significantly delayed and intracellular structures including mitochondria and plastids remained at mature pollen stage in the autophagy-deficient mutant Osatg7-1 Kurusu et al., 2014), suggesting that autophagy contributes to tapetal degradation and PCD in rice. However, the relationship between gene regulatory network and autophagy in plant anther development and tapetal PCD progression is mostly unknown.
In this report, to reveal the role of autophagy and how autophagy affects the gene regulatory network during anther/ pollen development and tapetal PCD, we performed RNAsequencing (RNA-seq)-based transcriptome analyses in combination with a quantitative PCR (qPCR) of anthers from ST8 to ST11 of the wild-type (WT) and the autophagy-deficient mutant Osatg7-1 plants. Role of autophagy in the synchronized progression of tapetal PCD, metabolisms of carbohydrates and lipids and the quality control of organelles during rice anther/ pollen development are discussed.

Plant Materials and Sample Preparation of Anthers
Surface-sterilized seeds of transgenic rice lines (Oryza sativa L. cv. Nipponbare (NB)) were germinated on MS medium (Murashige and Skoog, 1962) containing 0.8% agar and grown for 10 days in a growth chamber in long-day conditions (16-h light/8-h darkness, 28°C). Seedlings were transplanted into soil and grown in a greenhouse or paddy field at the National Institute of Genetics (NIG) in Mishima (Japan).

Sample Preparation for RNA Extraction
Anther samples at different developmental stages were separated based on the length and color of the anthers (Table 1), immediately frozen with liquid nitrogen in microtubes, and stored at −80°C. More than 180 anthers from 30 flowers were used at each stage sample.
Total RNA was extracted from stage 8 to stage 11 rice anthers, with three biological replicates each, using TRIzol reagent in accordance with the manufacturer's instructions (Life Technologies) and treated with DNase I (TaKaRa). Quality and integrity of RNA were checked using a spectrophotometer and Agilent 2100 bioanalyzer (Agilent).

RNA-sequencing and Gene Expression Profiling
RNA-seq libraries were prepared using the Illumina TruSeq ® Stranded RNA LT kit (Illumina), according to manufacturer instructions. To find differentially expressed genes (DEGs) throughout developmental stages in WT as well as Osatg7-1 mutant anthers, 30 libraries were prepared and sequenced using the NextSeq500 sequencing platform (Illumina) in accordance with the manufacturer's instructions. Approximately 20 million raw reads giving more than 3 Gb sequence data for each sample were obtained by single-end sequencing of 75 bp length. The data were deposited to the DDBJ Sequencing Read Archive database (Accession number DRA008977). The obtained reads were mapped to the reference rice genome (IRGSP-1.0) by TopHat2 (Kim et al., 2013), and the htseq-counts script in the HTSeq library was used to count the reads (Anders et al., 2015). Count data were subjected to trimmed mean of M-values (TMM) normalization in EdgeR (Robinson et al., 2010;McCarthy et al., 2012). Multidimensional scaling was performed via calculating log-fold changes between WT and using DEGs to compute distances in EdgeR with the "plotMDS" function (Supplementary Figure S1).
Transcript expression profiles and DEGs were defined using EdgeR generalized linear models (GLMs) (Robinson et al., 2010). Differential expression was calculated via fitting a GLM at the gene level using either developmental stages (from ST8 to ST11) or autophagy dependence as factors. The threshold for DEGs was a false discovery rate (FDR) of < 0.01; this yielded 20,391 genes. As a result, RAP-ID accessions were assigned in 19,748 genes of all DEGs, and these 20,391 genes were designated as antherspecific DEGs (ASDs). Moreover, for defining all DEGs under the control of autophagy activation during anther development, differential expression was calculated via fitting a GLM at the gene level using autophagy dependence as a factor. As a result, 2,359 genes were designated as autophagy-dependent DEGs (ADDs).

Bioinformatic Analyses: Principal Components Analysis with Self-Organizing Map Clustering and Gene Ontology Analysis
We applied a gene-expression clustering method (Chitwood et al., 2013) for both ASDs and ADDs defined using EdgeR. Scaled expression values were used for multilevel 3 × 3 and 4 × 4 rectangular self-organizing map (SOM) clusters (Figures 2 and 7) (Kohonen, 1982;Wehrens and Buydens, 2007). One hundred training interactions were used during clustering. Gene clusters were based on the final assignment of genes to winning units. In order to focus on gene-expression patterns instead of expression magnitude and to identify genes that vary in expression patterns between the WT and Osatg7-1, expression values were meancentered and variance-scaled separately between the WT and Osatg7-1 in a 3 × 3 rectangular SOM.  (Bastian et al., 2009) were used to visualize-as a directed network-the assignment of genes from different accessions to separate clusters. Direction of arrows indicates gene assignment to clusters, from the WT to Osatg7-1, with arrow size proportional to the gene number represented. Clustered and displaced gene sets among clusters were subjected to gene ontology (GO) analysis using agriGO v2.0 (http:// systemsbiology.cau.edu.cn/agriGOv2/).
The MapMan program allows the grouping of genes into different functional categories and visualization of data through various diagrams (Jung and An, 2012). To obtain functional classifications, we uploaded RAP locus IDs for 1,556 ADDs [early phase of stage 10 (ST10E); FDR < 0.05] to the MapMan tool kit. We then investigated the overviews of both regulation and metabolism pathways containing photosynthesis pathways, carbohydrate metabolism, N-dependent pathways such as amino acid, cell wall, lipid, and secondary metabolisms.

Quantification of mRNA by Real-Time PCR
First-strand complementary DNA (cDNA) was synthesized from 500 ng of total RNA with ReverTra Ace ® qPCR RT Master Mix with gDNA Remover (TOYOBO). Real-time PCR was performed using a Bio-Rad CFX Connect Real-Time System (Bio-Rad) with the THUNDERBIRD ® SYBR qPCR Mix (TOYOBO) and the specific primers (Supplementary Table S1). Relative mRNA levels were calculated using the 2 -DDCt method and normalized to corresponding an OsUbiquitin5 gene (Os01g0328400) level. The relative level of each gene in the WT anthers at tetrad stage (ST8) was standardized as 1.

Transcriptomic Analyses by RNA Sequencing
Autophagy has been shown to be induced at the uninucleate stage throughout the tapetal cells during anther development Hanamata et al., 2019). To identify genes affected by the activation of autophagy throughout the reproductive period, we conducted RNA-seq experiments using whole anther samples and compared the data between WT and Osatg7-1 plants. We obtained data from five different stages: tetrad (ST8), early uninucleate stage (ST9), late uninucleate stages (early and late phases of stage 10, ST10E and ST10L, respectively), and bicellular stage (ST11), each with three biological replicates (Table 1). Overall, 210,249,723 reads from the WT and 212,842,580 reads from the Osatg7-1 anthers were mapped to the rice genome. The reliability between each sample was confirmed using multidimensional scaling (MDS) plot analysis (Supplementary Figure S1). To confirm the trends in gene expression levels of anther samples from different stages between the WT and the Osatg7-1 mutant, we used correlation matrix (Figure 1). Similar expression patterns mainly depended on the developmental stage rather than on autophagy disruption. This trend was also confirmed by MDS plot analysis (Supplementary Figure S1).

Anther-Specific Differentially Expressed Genes and Induction of Autophagy-Related Genes During Pollen Maturation
To define all differentially-expressed genes (DEGs) throughout developmental stages in the WT and the Osatg7-1 mutant anthers, we used EdgeR package in R software. Of 38,311 genes located on the rice genome, 20,391 genes were defined as anther-specific DEGs (ASDs) based on a GLM at the gene level using either developmental stages (from ST8 to ST11) or a u t o p h a g y d e p e n d e n c e a s f a c t o r s ( F D R < 0 . 0 1 ) (Supplementary Dataset 1). The lists of ASDs that were significantly up-or down-regulated in each developmental stage are shown in Supplementary Dataset 1.
We found that the expression patterns of a variety of ATG genes including ATG8s defined as ASDs were up-regulated from ST10E to ST11 in the WT throughout anther development (Supplementary Figure S2). Ectopic overexpression of OsATG8s, such as ATG8a and ATG8b, has been shown to result in enhancement of autophagic activity in rice Zhen et al., 2019). Since autophagy is dramatically induced at the uninucleate stages (ST9-10) throughout the tapetal cells during pollen maturation , sustainedautophagic activity in rice tapetum may be regulated by transcription of ATG genes including ATG8s at least in part.
Visualization and Assessment of Self-Organizing Map Clustering Using the Anther-Specific Differentially Expressed Genes Profiles Self-organizing map (SOM) analysis allows us to identify a subset of genes with similar expression profiles (Nakayama et al., 2018). We next performed SOM clustering to further understand the differences in expression patterns. Gene expression values from the WT and the Osatg7-1 mutant were mean-centered and variance-scaled separately, allowing a focus on the differences in the expression patterns instead of expression magnitude. As a result, ASDs were assigned to clusters, and 16 clusters were obtained successfully (Figures 2A, B and Supplementary Dataset 2), based on box and line plots showing genes in each cluster with distinct, nonredundant expression patterns (Supplementary Figure S3). For instance, in cluster 9, the expression levels in ST8 and ST9 were lower than the average of the expression levels at each stage, and then the expression levels increased to ST10L. After that, their levels decreased at ST11. As shown in Figure 2C, the numbers of ASDs in the WT anthers were enriched in clusters 16 (3,189 DEGs; 16% of all ASDs in WT), 8 (1,843 DEGs; 9%), 9 (1,725 DEGs; 8%), and 4 (1,623 DEGs; 8%). In cluster 16, which has the most abundant of ASDs in the WT, the expression levels were lower than the average levels until ST10L, and then the expression levels drastically increased to ST11. In contrast, in the Osatg7-1 mutant, the numbers of ASDs in cluster 9 were significantly reduced (1,104 DEGs; 5% of all ASDs in the Osatg7-1 mutant). Conversely, the genes in clusters 15 (1,897 DEGs; 9%), 13 (1,751 DEGs; 9%), and 14 (1,533 DEGs; 8%) were increased compared with those in the WT ( Figure 2C). The common feature of these clusters (13, 14 and 15) was that the induction of gene occurs at ST10E.
We then focused on displaced gene sets between the WT and Osatg7-1 mutant among clusters based on SOM clustering results from ASD data ( Figure 3A). We extracted and visualized genes with different SOM cluster numbers in both the WT and Osatg7-1 from SOM clustering data excluding genes that were distant from the typical cluster pattern (distance < 0.8) (Supplementary Dataset 3). The 7,813 displaced genes exhibited certain tendencies ( Figure 3B; all directions from WT to Osatg7-1). As shown in Figure 3D, the numbers of displaced ASDs were enriched in clusters 4 and 9 of the WT, and these displaced gene sets were significantly enriched in displacements 9 ! 13 (718 genes; 9% of displaced genes) and 4 ! 8 (507 genes; 6%), suggesting that pre-and post-displacement differences in expression patterns occurs apparently at ST10E in anthers ( Figure 3C). This trend was also confirmed by the data in Figure 2. Moreover, the numbers of displaced ASDs were also enriched in clusters 1, 5 and 8 of the WT, and these displaced gene sets were significantly enriched in displacements 8 ! 3 (431 genes; 6%), 5 ! 9 (394 genes; 5%), and 1 ! 5 (355 genes; 5%) ( Figure 3D and Supplementary Figure S4). These results suggest that the state of gene expression patterns of anthers is delayed in the Osatg7-1 mutant compared with those of the WT.
For further characterization of these displaced gene sets between the WT and Osatg7-1 mutant, we performed gene ontology (GO) enrichment analysis using the agriGO v2.0 program. Among these, the GO term "Mitochondrial component" as cellular components was apparently enriched in the categories of displacements 9 ! 13, 4 ! 8, 8 ! 3, 1 ! 5, and 5 ! 9 (Figure 4). In contrast, the GO term "Plastid/ Chloroplast component" was enriched in the category of displacements 1 ! 5, 5 ! 9, and 9 ! 13. "Microtubuleassociated complex" was also enriched in the category of displacements 4 ! 8 (Figure 4). In the biological process group, the GO terms "Response to chemical stimulus" and  "Secondary metabolism" were enriched in the category of displacement 9 ! 13. "Cellular protein catabolic process containing ubiquitin-dependent protein" and "Cellular protein metabolic process" were enriched in the category of the displacement 4 ! 8. "Cell cycle process" was enriched in the category of the displacement 8 ! 3. "Cell wall organization" was enriched in the category of the displacements 9 ! 13 and 8 ! 3 (Figure 4). Moreover, in the molecular function, the GO term "Oxidoreductase activity" was ranked as the top category in displacement 1 ! 5 (Figure 4).
Overall, the present results indicate that the state of gene expression patterns in anthers is delayed in the autophagy-

Defining Differentially Expressed Genes Depending on Autophagy Activation in Anthers During Pollen Development
To identify DEGs affected by the activation of autophagy during pollen development, we used the R packages EdgeR. As a result, 2,359 DEGs throughout the developmental stages examined were defined as autophagy-dependent DEGs (ADDs) (FDR < 0.05) (Supplementary Dataset 4). In ST8 anthers, only 45 genes (1.9% in ADDs) were up-regulated and 108 (4.6%) were downregulated in the Osatg7-1 mutant compared with those in the WT. In ST9 anthers, 144 genes (6.1%) were up-regulated and 473 (20.1%) were down-regulated. In ST10E anthers, the effect of autophagy disruption showed the largest effect and 635 genes (26.9%) were up-regulated, while 921 (39.0%) were downregulated. Then in ST10L anthers, only 1 gene (0.04%) was upregulated and 9 (0.4%) were down-regulated, and in ST11 anther, 14 (0.6%) up-regulated and 215 (9.1%) down-regulated genes were identified ( Figure 5). These results indicate that the effect of autophagy disruption on gene expression throughout pollen development was largest at ST10E, exactly at the same stage when autophagy is induced in the tapetum Hanamata et al., 2019). Previous TEM analysis has shown that in Osatg7-1 mutant, morphology of tapetum at ST8-9 is normal, while the tapetal collapse is significantly delayed and remaining intracellular components including mitochondria and plastid are observed at ST11 . This is consistent with the present results that disruption of autophagy affects gene expression very little at ST8, while most severely at ST10E during pollen maturation ( Figure 5), suggesting that autophagy is crucial for timely progression of tapetal PCD, not for its initiation.

Differences in the Transcriptome Profile between the Wild-Type and Osatg7-1 Mutant Anthers During Pollen Maturation
To compare the expression profiles throughout pollen developmental stages between the WT and Osatg7-1 mutant anthers, we performed principal component analysis (PCA) using autophagy-dependent DEGs (ADDs; 2,359 genes). Major sources of variance in the transcriptome were investigated with a PCA that considered ADDs between the WT and Osatg7-1. As shown in Figure 6A, the first component (PC1) explained 53.3% of the variation and discriminated clearly between the WT and Osatg7-1 at ST10E. The second component (PC2) explained 20.3% of the variation and discriminated between the WT and Osatg7-1 at ST9 (Figures 6A, B; Supplementary Figure S5). These results of PCA clearly showed that the gene expression patterns were quite similar between the WT and Osatg7-1 at ST8, while quite different at ST9 and ST10E. Then the different patterns between the WT and Osatg7-1 became smaller thereafter. Of note, the expression patterns in the WT at ST8 were quite similar to those of Osatg7-1 at ST9, and the patterns of WT at ST9 were again similar to those of Osatg7-1 at ST10E.
Using heat maps, we illustrated the significantly up-or down-regulated transcripts (Osatg7-1 vs. WT; 2,359 genes in Supplementary Dataset 4). Severe changes in the mRNA profiles throughout the reproductive period were evident, with most transcripts being consistently affected ( Figure 6C).

Effects of Autophagy Disruption on the Metabolism of Carbohydrates, Lipids and Phytohormones as well as Quality Control of Organelles During Pollen Maturation
To identify genes under the control of autophagy during anther development precisely, we next performed SOM clustering to further understand the difference in expression patterns depending on autophagy disruption. We constructed a SOM to extract genes depending on the count data of ADDs. Gene expression values from Osatg7-1 were mean-centered and variance-scaled using the expression data from the WT (Sum of cpm through all stages in WT and Osatg7-1 > 10; 2,149 genes; Supplementary Dataset 5). As a result, all ADDs were assigned to clusters, and 9 clusters were successfully obtained ( Figures  7A, B; Supplementary Figure S6). The numbers of ADDs in the Osatg7-1 mutant were enriched in clusters 2, 3, 7, and 8, which showed peak in gene expression fluctuation is ST10E ( Figures  7A, B; Supplementary Figure S6), indicating that autophagy disruption affects gene expression levels throughout anther development and most severely at ST10E.
For further characterization of each cluster, we performed a GO enrichment analysis with the 9 clustered gene sets ( Figure 8). In the biological process groups, GO terms "Aromatic amino acid family metabolic process", "Carbohydrate catabolic process", "Cellular amino acid metabolic process", and "Catabolic process containing carbohydrate" in cluster 1 (q < 0.05), GO term "Cell cycle process" in cluster 8, and "Cell wall organization" and "Carbohydrate metabolic process" in cluster 9 were enriched (Figure 8). Of note, GO terms "Mitochondrial component" and "Plastid/Chloroplast component" as cellular components were also enriched in almost all clusters (Figure 8).
We then analyzed the overview in the MapMan tool to investigate the functions of genes differentially expressed between the WT and Osatg7-1 mutant during tapetal PCD process, especially in ST10E. Of the 1,556 gene that significantly changed (FDR < 0.05) in ST10E of ADDs, 1,500 genes were annotated and subjected to MapMan pathway analysis. As a result, 166 genes were mapped to the FIGURE 5 | Histogram showing the number of up-regulated and down-regulated ADDs during anther developmental stages of Osatg7-1 compared with WT. The number of genes two-fold up-regulated or down-regulated at FDR < 0.05 are plotted on the red or blue bars in the graph, respectively. metabolism overview, and 351 genes were mapped to the regulation overview, respectively (Figure 9). Almost all metabolism-related genes are categorized as 24 genes for cell wall metabolism, 24 for lipid metabolism, 40 for secondary metabolism, 21 for amino acid metabolism, 23 for carbohydrate metabolism, 10 for nucleotide metabolism and 5 for photosystems including light reaction and Calvin-Benson cycle ( Figure 9A).
Regulation-related genes were categorized as 102 genes for TFs, 48 genes for protein modification, 111 for protein degradation, 18 for hormone metabolism, and 57 signaling-related genes ( Figure 9B).
Hormone metabolism was further classified as 4 genes associated with IAA, 2 with abscisic acid (ABA), 9 with ethylene, 1 with cytokinin, 1 with jasmonate (JA), and 2 with gibberellin (GA) ( Figure 9B). These results suggest that many genes affected by autophagy were associated with cell wall, lipid and carbohydrate metabolic processes as well as some phytohormones.
Hormone profiling analyses at the flowering stage have shown that endogenous levels of active-forms of GAs, were significantly lower in the Osatg7-1, than in the WT . Treatment with GA 4 partially rescued the phenotype of pollen maturation in Osatg7-1 . In addition, the expression levels of GA20-oxidases, which have been reported to be required for GA biosynthesis in rice (Sasaki et al., 2002), were clearly decreased in Osatg7-1 at ST10E in comparison with that of WT ( Figure 9B and Supplementary Figure S7), indicating that autophagy affects GA biosynthesis during pollen maturation in rice. Precursors of bioactive GAs are synthesized in plastids (Lange, 1998;Hedden and Phillips, 2000). Indeed, the GO term "Plastid/ Chloroplast component" was enriched in almost all clusters, and photosynthesis-related genes including those involved in light reaction and Calvin-Benson cycle were decreased in the Osatg7-1 mutant at ST10E (Figures 8 and 9A). Since autophagy plays an important role in the quality control of plastids including elimination and turnover of photodamaged chloroplasts Izumi et al., 2017), the lower activity of GA biosynthesis in anther may be attributed to the accumulation of abnormal plastids in Osatg7-1 due to the disruption of autophagy. The relationship between autophagy and quality control of organelles throughout the progression of tapetal PCD should be an important future research topic.
Transcriptome analyses in the developmental process of rice anthers have shown that the genes involved in lipid and secondary metabolisms are expressed in tapetal cells and microspores at the uninucleate stages (ST9-10) (Hobo et al., 2008;Deveshwar et al., 2011). Autophagy has been shown to contribute to amino acid metabolism, lipid turnover, and maintenance of secondary metabolic pathways in nutrientstarved conditions in Arabidopsis and maize leaves (Masclaux-Daubresse et al., 2014;McLoughlin et al., 2018). In the present transcriptome analyses followed by the SOM and GO analyses of ADDs, several genes associated with both amino acid and carbohydrate metabolic/catabolic processes were enriched in clusters 1 and 9 ( Figure 8). Moreover, in the SOM and GO analyses of ASDs, the expression patterns of these genes were enriched in cluster 4 in the WT anthers (Figure 4), while moved to cluster 8 in the Osatg7-1 anthers ( Figure 4). TEM analysis has shown that accumulation of both starch granules and lipid bodies in pollen grains are impaired in Osatg7-1 at the mature stage . Defects in pollen maturation in Osatg7-1 anthers may be associated with the down/delayed-regulation of carbohydrate metabolism-related genes.
Lipidic exine synthesis is an important component of the pollen wall in rice and Arabidopsis (Yang et al., 2007). At the flowering stage, development of the coat structure of pollen grains and accumulation of lipid bodies in pollen grains are impaired in Osatg7-1 . Moreover, MapMan analysis clearly showed that autophagy disruption affects the expression of genes involved in lipid and cell wall metabolic processes at ST10E in anthers ( Figure 9). Therefore, the expression patterns of several marker genes related to pollen sporopollenin biosynthesis were checked by RNA-seq as well as qPCR analysis ( Figure 10). Induced expression of CYP704B2, which has been reported to be required for sporopollenin and pollen wall biosynthesis in rice (Li et al., 2010), was lower at ST10E-10L in the Osatg7-1 than WT ( Figure 10). Moreover, the expression of a rice homolog of MYB80, which is required for tapetal PCD and anther development in Arabidopsis (Phan et al., 2011), was also lower and delayed at ST9-10E ( Figure  10). Induced expression of PTC1 encoding a key TF regulator of sporopollenin biosynthesis, was also delayed in the Osatg7-1 at ST9 (Figure 11), suggesting the involvement of autophagy in sporopollenin biosynthesis and pollen wall formation during pollen maturation. Moreover, a lipid transfer protein of rice, anther specific protein 4 (C4), is specifically expressed in tapetal cells (Tsuchiya et al., 1994) and is able to transport lipid molecules such as fatty acids from tapetal cells to developing microspores. The induction of C4 gene was also suppressed in the Osatg7-1 compared with WT at ST10E (Figure 10). These results suggest that the immature pollen phenotype of Osatg7-1 is attributed to the defects in timely transport of lipids from the tapetum cells to developing microspores, which may lead to less

Effects of Autophagy Disruption on the Progression of Tapetal PCD During Pollen Maturation
Appropriate temporal regulation of tapetal PCD is vital for normal pollen development in plants. The signal initiating tapetal PCD has been suggested to be first produced during the tetrad stage (ST8) (Kawanabe et al., 2006). A transcriptional regulatory network as well as activation of proteases is also known to play a key role in the progression of tapetal PCD, and some key genes involved in the tapetal PCD include the PHD-finger transcription factors (PHD-TFs) and the bHLH transcription factors (bHLH-TFs) as well as aspartic proteases (APs) and cysteine proteases (CPs) (Niu et al., 2013;Ono et al., 2018;Yang et al., 2019).
To investigate the effects of autophagy disruption on the transcriptional regulatory networks during tapetal PCD process, we searched the expression profiles of PHD-TFs and bHLH-TFs as well as CPs and APs genes in rice (Chen et al., 2009;Wang et al., 2018). We extracted 5 PHD-TFs and 10 bHLH-TFs as well as 9 APs and 4 CPs (ADDs; FDR < 0.05). The expression level of PTC1, which has been reported to be required for rice tapetal PCD , was decreased and delayed from ST8 to ST9 in the Osatg7-1 compared with WT ( Figure 11). The induction of selected bHLH-TFs was also decreased in the Osatg7-1 anthers throughout tapetal PCD (Figure 11), and these characteristic expression patterns of some key genes, such as EAT1 involved in PCD progression were also confirmed by qPCR analysis (Figure 12). Moreover, the induction of GAMYB, which has been reported to be required for tapetal PCD and pollen development in rice (Aya et al., 2009), was also suppressed at the ST10-11 in the anthers of Osatg7-1, suggesting that autophagy disruption alters the expression patterns of the key regulatory TFs associated with PCD throughout the progression of tapetal PCD.
EAT1 is known to trigger tapetal PCD by regulating the expression of two aspartic proteases (AP25 and AP37) at the young microspore stage (Niu et al., 2013). The induction of EAT1 gene was suppressed at the ST10E in the anthers of Osatg7-1. Previous TEM analysis indicated the delay of tapetal PCD in Osatg7-1 anthers , which was correlated with the down-regulation of EAT1 expression at pollen developmental stages (Figures 11 and 12). Moreover, the expression levels of 4 CP genes were significantly down-regulated in Osatg7-1 anthers compared with those of the WT throughout tapetal PCD progression ( Figure 11). The induction of AP25, which has been reported to be required for rice tapetal PCD (Niu et al., 2013), was also delayed and decreased at ST9-10E in Osatg7-1 compared with the WT (Figures 11 and 12). Wild-type anthers exhibited strong induction of EAT1, AP25 as well as AP38 at ST9 to turn on timely progression of tapetal PCD (Figures 11 and  12). As shown in Figure 7, the expression patterns of almost all ADDs apparently fluctuated at ST10E. These trends were also confirmed using heatmap analysis, SOM clustering as well as qPCR analysis (Figure 12). Thus, the impairment and delay of FIGURE 10 | Verification of expression profiles of selected genes related to lipid metabolism by qPCR analysis. Quantitative expression levels of CYP704B2 (A) and MYB80 (B) and C4 (C) during anther developmental stages. The amount of each mRNA was calculated from the threshold point located in the log-linear range of RT-PCR. The relative level of each gene in the WT anthers at tetrad stage (ST8) was standardized as 1. Data are the means ± SE of three independent experiments. *P < 0.05, **P < 0.01; significantly different from the controls. tapetal PCD in the autophagy-deficient mutant may at least in part be attributed to the reduced expression of these genes involved in the progression of tapetal PCD at ST9-10 in anthers.

Future Perspectives
Autophagy-deficient mutants in Arabidopsis and rice are reported to be hypersensitive to oxidative stress (Xiong et al., 2007a;Xiong et al., 2007b;Shin et al., 2010), suggesting that autophagy plays an important role in oxidative stress responses during anther development. As an important protein quality control mechanism, autophagy has been suggested to play a critical role in the degradation of these misfolded/denatured and potentially highly toxic proteins or protein aggregates (Kraft et al., 2010;Johansen and Lamark, 2011;Shaid et al., 2013). Since autophagy is involved in the turnover of damaged mitochondria and chloroplasts Izumi et al., 2017), autophagy The heatmap shows log2-fold change in expression of ADDs involved in tapetal PCD during anther developmental stages of Osatg7-1 compared with WT. For SOM clustering of ADDs, both 4 × 4 (A) and 3 × 3 (B) rectangular topologies are used, and extracted genes are assigned to each SOM cluster, respectively.
FIGURE 12 | Verification of expression profiles of selected genes related to tapetal PCD progression and pollen maturation by qPCR analysis. Quantitative expression levels of TIP2 (A) and EAT1 (B) and GAMYB (C) and AP38 (D) and AP25 (E) during anther developmental stages. The amount of each mRNA was calculated from the threshold point located in the log-linear range of RT-PCR. The relative level of each gene in the WT anthers at tetrad stage (ST8) was standardized as 1. Data are the means ± SE of three independent experiments. *P < 0.05, ***P < 0.001; significantly different from the controls. disruption in anthers could cause over-accumulation of reactive oxygen species (ROS) under various environmental stress. Autophagy induced during tapetal PCD may play critical roles in the quality control of organelles and environmental stress responses including oxidative stress during pollen maturation.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the DRA008977.