Next Article in Journal
Myotonic Dystrophies: A Genetic Overview
Previous Article in Journal
Epigenetic and Physiological Responses to Varying Root-Zone Temperatures in Greenhouse Rocket
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Genes
Volume 13
Issue 2
10.3390/genes13020366
genes-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptomics-Based Identification of Genes Related to Tapetum Degradation and Microspore Development in Lily

by
Juanjuan Sui
1,
Wenjie Jia
2,3,
Yin Xin
2 and
Yuanyuan Zhang
4,*
1
Department of Biology, Biology and Food Engineering College, Fuyang Normal University, Fuyang 236037, China
2
Beijing Key Laboratory of Development and Quality Control of Ornamental Crops, Department of Ornamental Horticulture and Landscape Architecture, China Agricultural University, Beijing 100193, China
3
Flower Research Institute, Yunnan Academy of Agricultural Sciences, Kunming 650205, China
4
College of Biological Sciences, China Agricultural University, Beijing 100193, China
*
Author to whom correspondence should be addressed.
Genes 2022, 13(2), 366; https://doi.org/10.3390/genes13020366
Submission received: 24 January 2022 / Revised: 7 February 2022 / Accepted: 12 February 2022 / Published: 17 February 2022
(This article belongs to the Section Plant Genetics and Genomics)
Browse Figures
Review Reports Versions Notes

Abstract

:
Lily is a popular and economically ornamental crop around the world. However, its high production of pollen grains causes serious problems to consumers, including allergies and staining of clothes. During anther development, the tapetum is a crucial step for pollen formation and microspore release. Therefore, it is important to understand the mechanism of tapetum degradation and microspore development in lily where free pollen contamination occurs. Here, we used the cut lily cultivar ‘Siberia’ to characterize the process of tapetum degradation through the use of cytology and transcriptomic methods. The cytological observation indicated that, as the lily buds developed from 4 cm (Lo 4 cm) to 8 cm (Lo 8 cm), the tapetum completed the degradation process and the microspores matured. Furthermore, by comparing the transcriptome profiling among three developmental stages (Lo 4 cm, Lo 6 cm and Lo 8 cm), we identified 27 differentially expressed genes. These 27 genes were classed into 4 groups by function, namely, cell division and expansion, cell-wall morphogenesis, transcription factors, LRR-RLK (leucine-rich repeat receptor-like kinases), plant hormone biosynthesis and transduction. Quantitative real-time PCR was performed as validation of the transcriptome data. These selected genes are candidate genes for the tapetum degradation and microspore development of lily and our work provides a theoretical basis for breeding new lily cultivars without pollen.

1. Introduction

Lily is one of the most important cut flowers in the world due to its beautiful blossoms and pleasant fragrance [1,2]. However, the pollens released from anthers cause environmental pollution leading to consumer issues, such as pollen allergies and staining of clothes and petals [1,3,4]. Although doubled lily varieties have no stamens and grow additional petals instead, consumers still prefer single lilies, such as popular cultivars ‘Sorbonne’ and ‘Siberia’, which still maintain their anthers and stamens [1,4]. Therefore, a sterile single lily would be the ideal cultivar for consumers.
A mature anther is composed of four somatic layers [5,6]. From the outside to the inside, they are the epidermis, the endothecium, the middle layer and the tapetum [5,6]. The tapetum is the innermost layer to the anther chamber and it plays an essential role in microspore development. The tapetum provides nutrition for microspores, degrades callose and participates in pollen exine and coating formation [5,6,7]. Previous studies have shown that there are several genes involved in tapetum degradation and microspore development in model plants, e.g., Arabidopsis and rice. In the myb33/myb65 mutant, tapetum cells vacuolated at the pollen mother-cell stage, leading to microspore degradation before meiosis [8]. In the dyt1 (dysfunctional tapetum 1) mutant, although normal tapetum cells could form, the tapetum cells were non-functional [9]. At the 6th and 7th stages of anther development in Arabidopsis, the tapetal cells were abnormally hypertrophic and filled the chamber, squeezing the developing tetrad, which led to the early degradation of microspores [9]. Similarly, in the tdf1 (defective in meristem development and function 1) mutant, the tapetum showed severe vacuolization and irregular division, followed by abnormal hypertrophy, leading to microspore degeneration [10]. In the ams (aborted microspores) mutant, the microspore mother cells experienced normal meiosis, but, at a later stage, the tapetum vacuolated and hypertrophied abnormally after the microspore was released from the tetrad, which led to the early degradation of microspores [11,12,13]. In the ms188 (male sterile 188) mutant, pollen mother cells developed normally in the early stage; however, the tapetum could not develop into normal secretory cells at a later stage, which also led to the early degradation of microspores [6,14,15]. In the ms1 mutant, the tapetum cells vacuolated, which led to the dysfunction of the tapetum and the degradation of microspores [16,17]. Therefore, abnormal tapetum development can result in pollen abortion.
Male sterility reduces the work of artificial emasculation and solves the problem of pollen contamination and allergy in ornamental plants, such as lily [3]. Although much progress has been made in the field of anther development, the mechanism largely remains to be investigated [13]. More recently, RNA sequencing (RNA-Seq) has been used as a tool in ornamental plants to identify genes related to a variety of desirable traits [18].
To better understand the degradation of the tapetum and independent microspore releasing, we histologically analyzed lily anthers by paraffin sectioning. According to the sections analyzed, anthers at three different stages were selected for transcriptome sequencing, namely, Lo 4 cm (flower buds at 4 cm in length), Lo 6 cm and Lo 8 cm. The BGISEQ-500 sequencing platform was used to identify the differentially expressed genes (DEGs) related to tapetum degradation and independent microspore releasing in ‘Sibera’. Several important DEGs were identified, including some transcription factors (TFs). Our findings provide new insights into the potential molecular mechanism of tapetum degradation and microspore development. In addition, the candidate genes identified can be used for breeding new lily varieties without pollen contamination by genetic engineering.

2. Materials and Methods

2.1. Plant Materials and Treatments

Anthers collected from the Oriental Hybrid Lily cultivar ‘Siberia’ were used in the experiment. The lily was planted in greenhouses and grown under a photoperiod of 16 h of light and 8 h of darkness, at 22 °C, at Fuyang Normal University, Fuyang, Anhui Province, China. The anthers were collected from the flower buds in different sizes (1, 2, 3, 4, 5, 6, 7, 8, 9 and 10 cm in length). Each lily bud contained 6 anthers, one of which was fixed with FAA (Formalin-Aceto-Alcohol) solution and was used for histological analysis, while the remaining five were quickly frozen in liquid nitrogen and transferred to −80℃ to be stored for RNA sequencing and analysis of gene expression.

2.2. Lily Anther Histological Observation

Anthers embedded in paraffin were cut into 8 μm thick sections, corresponding to a size of about 0.5 mm × 0.5 mm × 0.5 mm. The processes were performed as previously described [4]. The sections were stained by safranine and fast green double dyeing. Section observation was performed under a light microscope (Nikon Eclipse E100; Tokyo, Japan).

2.3. RNA Extraction and Library Construction

The anthers were removed from the buds and total RNA was extracted from these samples using an RNAprep Pure Kit (DP441, Tiangen, Beijing, China). More than 20 μg of RNA was extracted from each group of samples and was used to construct cDNA libraries. The quality and quantity of the RNA samples were assessed using Qubit and Agilent 2100. A Clontech SMARTer PCR cDNA Synthesis Kit was used to construct the Iso-Seq library. The Iso-Seq library construction started with the mRNA purified with Oligo(dT) magnetic beads (NEB) and cDNA was synthesized by using the SMARTer PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA) according to the Iso-Seq protocol (Pacific Biosciences, Menlo Park, CA, USA). Then, cDNA was amplified using KAPA HiFi DNA Polymerase (Roche, Basel, Switzerland) for 12 cycles followed by purification and size selection. The BluePippin Size Selection System protocol was performed by following the instructions described by Pacific Biosciences (PN 100-092-800-03). The insert size of the library was detected using Agilent 2100 to ensure the quality of libraries. Short-read RNA-Seq libraries were prepared using TruSeq stranded RNA Library Preparation kits and the supplied protocol (Illumina, San Diego, CA, USA) and sequenced on a HiSeq2500 platform using v2 sequencing chemistry to generate 2 × 150 paired-end reads.

2.4. Error Correction of PacBio Iso-Seq Full-Length cDNA Reads

The samples’ transcriptomes were sequenced using the Illumina NextSeq 500 (Illumina, San Diego, CA, USA) and Pacbio SMRT Sequel (Pacific Biosciences, Menlo Park, CA, USA) platforms for RNA-Seq and Iso-Seq, respectively. The Iso-Seq subread raw reads were firstly obtained from the raw SMRT raw long reads and used to generate the corrected circular consensus sequences (CCS) by SMRT Link v8.0 software (Pacific Biosciences, Menlo Park, CA, USA) (minLength, 50; maxLength, 15,000; minPasses, 1). Then, the sequences were divided into full-length sequences and non-full-length sequences according to the existence of the 5′-primers, 3′ primers and polyA tails. To improve consensus accuracy, full-length non-chimera (FLNC) reads were polished by Illumina short reads using the arrow and consensus transcripts were obtained using the isoform-level clustering algorithm ICE (iterative clustering for error correction). Finally, RNA-Seq short reads from the nine samples were employed to correct sequencing errors in the consensus transcripts using the program LoRDEC [19].

2.5. Gene Prediction and Annotation

Cogent (https://github.com/Magdoll/Cogent, accessed on 23 June 2021; v 2.1) was used to reconstruct isoforms into one or several full-length unique transcript model(s) based on the k-mer clustering and De Bruijn graph methods. It generates UniTransModels as distinct genes. GMAP [20] software was used to compare the transcripts to UniTransModels and SUPPA software was used to identify different alternative splicing types. Transcripts with different splicing sites were considered as alternative splicing isomers. Finally, highly similar isoforms were clustered and removed using CD-HIT [21]; the non-redundant representative transcripts were kept for further analyses.
To obtain comprehensive gene function information, the non-redundant representative transcripts were aligned to the nucleotide and protein databases, i.e., Nr, Pfam, KOG/COG, Swiss-Prot, KEGG and GO. Open reading fragments (ORFs) for each transcript were detected by TransDecoder with a minimum length of 100 amino acids (AAs). Where multiple ORFs were found in a single transcript, we used the first ORF in each transcript as the ‘representative’. In each transcript locus, the isoform containing the longest ORF was the ‘locus representative’ for functional annotation. Protein sequences were downloaded from Swiss-Prot (ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz, accessed on 23 June 2021) and NCBI NR (ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/, accessed on 23 June 2021). We then removed the records for proteins described as hypothetical or predicted in these well-characterized protein databases for the ‘locus representative’ protein sequences using BLASTp (NCBI-BLAST v2.2.28+) with the following parameters: outfmt, 6; max_target_seqs, 1; evalue 1e-5. Meanwhile, the Pfam database of 16,712 protein families was downloaded from (ftp://ftp.ebi.ac.uk/pub/databases/Pfam; release Pfam31.0, accessed on 23 June 2021). Pfam Scan was used to identify occurrences of the known protein domains documented in the Pfam database with the parameters ‘E 1e-5’. Transcription factors were predicted using iTAK v 4.0 [22]. The GO terms and pathways in KEGG for each gene were obtained from KOBAS using Arabidopsis thaliana protein sequences as reference [23].

2.6. Quantification of Transcripts Using Illumina Short Reads

Paired-end Illumina RNA-Seq reads from each tissue sample were trimmed to remove the adaptor sequences and low-quality bases using Trimmomatic (v0.32) [24] with the explicit option settings ILLUMINACLIP: adapters.fa: 2:30:10:1:true LEADING:3 TRAILING:3 SLIDINGWINDOW: 4:20 LEADING:3 TRAILING:3 MINLEN:25. Trimmed Illumina reads were aligned against the obtained transcripts from the Iso-Seq analysis by Bowtie2 [25]. RSEM (v1.3.3) [26] was used to calculate the FPKMs of all transcripts in each sample. Isoforms FPKMs were computed by summing the FPKMs of transcripts in each gene group. FPKMs mean fragments per kilo-base of exon per million fragments mapped, calculated based on the length of the fragment and read count mapped to this fragment.
DEGs were identified using the DESeq2 package with the threshold of logarithm transformation of fold change larger than 1 and p-value < 0.05 [27].

2.7. The ‘eFP-Seq’ Browser Analysis of DEGs

The Arabidopsis ‘eFP-Seq Browser’ tool is a useful tool to display transcriptome data visually through electronic fluorescent pictographic (eFP) images, which can present the relative expression level of genes in different stages of development (http://bar.utoronto.ca/efp/cgi-bin/efpWeb.cgi, accessed on 23 June 2021) [28]. In this study, tissue-specific data, including pollen development, were analyzed.

2.8. Quantitative Real-Time PCR

To further verify the gene expression levels of DEGs identified by RNA sequencing, quantitative real-time PCR was used to detect gene expression to further confirm the transcriptome results. All primers of genes are listed in Table S1. The Lo18S rRNA was used as the internal reference for the quantitative expression analysis. Three biological replicates were performed. The 2−ΔΔCT method was used for the data analyses [29].

3. Results

3.1. Lily Anthers Cytology Observation

The lily flower contains six anthers, each with four chambers. The anther chamber is the location of microspore development and maturation. When the anther is mature, it cracks naturally and pollen grains scatter from the chamber (Figure 1A). To track the developmental stage for tapetum degradation and independent microspore formation, we collected anthers of different sizes of flower buds (from 1 to 10 cm) and took out one of the six anthers from each flower for paraffin embedding and microscopic analysis. The histological results show that the four layers’ structure (tapetum, middle layer, endothecium and epidermal layer) gradually formed from 1 cm to 3 cm stage buds. The tapetal cells became long with a large nucleus and dense cytoplasm in the 3 cm buds, which is in accordance with previous results [4]. At the Lo 4cm bud stage, several changes occurred in the anther, i.e., (1) binuclear and polynuclear cells appeared; (2) cytoplasm condensed with cellular vacuolization; and (3) the pollen mother cells developed into dyads and tetrads (Figure 1B,C). From the Lo 4 cm to Lo 6 cm stages, the tapetal cytoplasm was further concentrated and independent microspores were released from the tetrads, suggesting callose degradation during this transition (Figure 1B,D) [4,14]. Thus, the transition from the Lo 4 cm to Lo 6 cm stages was the key stage of callose degradation. At the Lo 8cm stage, the whole cell structure of the tapetum could not be observed except for parts of the tangential wall (Figure 1B,E). The tapetum provides protein, lipid and other nutrients for microspores development; consequentially, it was completely degraded before pollens were mature [5]. We observed microspore development from binuclear cells to mature pollen from the Lo 4 cm to Lo 8 cm stages. Taken together, alongside tapetum degradation, microspores developed and matured.

3.2. Overall Analysis of Transcriptome Sequencing and Gene Annotation

To obtain accurately expressed transcripts from the samples, we conducted Iso-Seq of pooled samples and RNA-Seq of nine samples from three development stages with three biological replicates. After transcript assembly and sequence correction, a total of 35,213 transcripts (17,065 genes) were identified and principal compound analysis (PCA) indicated high agreement among biological replicates (Figure 2A and Figure S1). The maximum length of unigenes was nearly 8.9 kb and the minimum length was 87 bp. N50 and N90 are important parameters reflecting the quality of transcript assembly [30]. The N50 and N90 of the transcripts were 2356 bp and 1237 bp, respectively, suggesting that the sequence quality was high and could be used for subsequent analyses. Furthermore, more than 95% of transcripts were less than 4 kb in length (Figure 2B). Gene annotation showed that there 97.07% (16,566) of genes were found at least in one database and 7553 transcripts were annotated by all the seven databases used for the analyses (NR, SwissProt, KEGG, KOG, GO, NT and Pfam) (Figure 2C). Finally, a total of 16,380, 14,287, 16,434, 12,984, 10,947, 11,929 and 14,287 genes were annotated in the NR, SwissProt, KEGG, NT, KOG, GO and Pfam databases, respectively (Figure 2C).

3.3. Screening of Differentially Expressed Genes (DEGs) and Co-Expression Analysis

To identify genes involved in tapetum degradation and independent microspore development, we analyzed the DEGs among three anther development periods, including Lo 4 cm, Lo 6 cm and Lo 8 cm. Here, DEGs were selected according to their possible functions and fold changes in the expression levels, as well as the screening threshold of padj < 0.05, which were based on the NR annotation database.
A total of 6896 DEGs was screened among the three stages (Figure 3A). In addition, a heat map representation of the expression pattern of the 6896 DEGs, including a hierarchical clustering analysis, can be seen in Figure 3B. The heat map analysis showed that most DEGs among the Lo 4 cm, Lo 6 cm and Lo 8 cm stages showed anti-correlated expression patterns (i.e., they were exactly opposite) (Figure 3B). Furthermore, in total, 4592 DEGs, 2595 DEGs and 4901 DEGs were identified in the comparisons Lo 6 cm vs. Lo 4 cm, Lo 8 cm vs. Lo 6 cm and Lo 8 cm vs. Lo 4 cm, respectively (Figure 3C). The number of up- and down-regulated genes was also counted in each group comparison (Figure 3C). The number of DEGs in the comparison Lo 6 cm vs. Lo 4 cm was almost twice that in Lo 8 cm vs. Lo 6 cm, indicating that the transition from the Lo 4 cm to Lo 6 cm stages was important for tapetum degradation and microspore development, which is consistent with the results of the cytology observations (Figure 1).
Next, we conducted a co-expression analysis of 6896 DEGs and identified 9 different clusters of gene expression patterns (Figure 3). The DEGs were divided into nine clusters and, within the same cluster, the DEGs had similar expression patterns (Figure 3D; Table S2). In cluster 1, the DEGs were highly expressed from Lo 4 cm to Lo 6 cm but were down-regulated from Lo 6 cm to Lo 8 cm, where the genes related to plant hormone signal transduction were enriched (Table S2). In cluster 2, the DEGs were only highly expressed at the Lo 6 cm stage, enriching metabolic genes related to fatty acid degradation, amino acid and energy metabolism, MAPK signal pathway, autophagy and peroxisome (Table S2). In cluster 3, the DEGs were down-regulated from Lo 4 cm to Lo 6 cm and had low gene expression levels from Lo 6 cm to Lo 8 cm, containing genes related to DNA replication, translation and photosynthesis (Table S2). In cluster 4, the DEGs were down-regulated at Lo 6 cm and enriched genes were involved in translation. In cluster 5, the DEGs were continuously down-regulated from Lo 4 cm to Lo 8 cm, especially genes that participated in photosynthesis, spliceosome, pentose phosphate pathway, phenylpropanoid biosynthesis and flavonoid biosynthesis (Table S2). In cluster 6, the DEGs had low expression levels from Lo 4 cm to Lo 6 cm but were highly expressed at Lo 8 cm and only one pathway was enriched, i.e., basal transcription factors (Table S2). In cluster 7, the DEGs were continuously up-regulated from Lo 4 cm to Lo 8 cm, especially genes involved in oxidative phosphorylation, lysosome, ubiquitin-mediated proteolysis and autophagy (Table S2). In cluster 8, the DEGs were down-regulated from Lo 4 cm to Lo 6 cm but were up-regulated from Lo 6 cm to Lo 8 cm and genes specific to the regulation of phagosome were enriched (Table S2). In cluster 9, the DEGs were up-regulated from Lo 4 cm to Lo 6 cm and had high expression levels from Lo 6 cm to Lo 8 cm, with enriched pathways including peroxisome, fatty acid, ether lipid metabolism, citrate cycle and oxidative phosphorylation (Table S2). Based on these co-expression analyses, it can be suggested that there are many genes involved in tapetum degradation and independent microspore development; regulating protein biosynthesis, transport and degradation; lipid acid, plant hormones and secondary metabolism; and further controlling cell division, cell death and cell differentiation.

3.4. GO Enrichment and KEGG Pathway Analysis

We performed a GO (gene ontology) enrichment analysis and identified that DEGs were enriched within three major categories, namely, biological process (BP), cellular component (CC) and molecular function (MF) [31]. It was worth mentioning that we thought anther development was a continuous process, so we did not analyze the DEGs in the comparison between Lo 8 cm and Lo 4cm. Additionally, we used the volcano map to visualize the distribution of DEGs in the comparisons Lo 6 cm vs. Lo 4 cm and Lo 8 cm vs. Lo 6 cm (Figure 4C; Tables S3 and S4).
In the comparison Lo 6 cm vs. Lo 4 cm, BP category contained 14 annotations, of which the primary terms were: single-organism process and single-organism metabolic process (Figure 4A). The MF category contained 20 annotations, with the main ones being catalytic activity, transferase activity and binding functions. However, few DEGs were enriched in the CC category (Figure 4A). Compared with DEGs in Lo 4 cm vs. Lo 6 cm, fewer genes could be enriched in the comparison between Lo 6 cm and Lo 8 cm and DEGs were only enriched in the BP and MF categories. The BP category contained nine annotations, with the main ones being metabolic process, single-organism process and single-organism metabolic process, which is consistent with the results of Lo 6 cm vs. Lo 4 cm (Figure 4B). The MF category also contained nine annotations, including the primary terms catalytic activity, oxidoreductase activity and transporter activity (Figure 4B). The above GO enrichment analysis suggests that the occurrence of genes related to metabolism and enzyme catalysis could affect tapetum degradation and microspore development.
Congruently, we also analyzed pathways enrichment of DEGs in the KEGG (Kyoto Encyclopedia of Genes and Genomes) database, a resource for understanding high-level functions and utilities of the biological system [32]. Here, we listed the top 20 KEGG pathways enrichment results obtained from the comparisons of the Lo 6 cm vs. Lo 4 cm and Lo 8 cm vs. Lo 6 cm stages (Figure 5A,B). From our data, six pathways, ‘alpha-Linolenic acid metabolism’, ‘Pyruvate metabolism’, ‘Glycolysis/Gluconeogenesis’, ‘Phenylpropanoid biosynthesis’, ‘Fatty acid degradation’, ‘beta-Alanine metabolism’ and ‘Carbon fixation in photosynthetic organisms’ that were related to substance biosynthesis and metabolism, such as amino acids and lipids, were significantly enriched in both comparisons, which is consistent with the GO enrichment results (Figure 4 and Figure 5). In addition to common pathways, we also found that, during the period from the Lo 4 cm to Lo 6 cm stages, 47 genes were enriched in the TCA cycle, an essential pathway to provide energy for the organism (Figure 5A) [33]. In the comparison of the Lo 8 cm vs. Lo 6 cm stages, a total of 32 DEGs were enriched in the ‘Plant hormone signal transduction’ pathway (Figure 5B). Considering that plant endogenous hormones such as GA and IAA change rapidly during tapetum degradation [34], it may play conserved functions in regulating the anther development in lily.

3.5. Analysis of Transcription Factors Related to Tapetum Degradation and Microspore Development

Transcription factors (TFs) can regulate the expression of genes associated with anther development [35]. For example, an important MADS-box transcription factor, AGAMOUS (AG), is able to determine the emergence of stamen primordia. The ag mutant cannot form the tapetum and endothelium layer normally [36,37,38]. Several MYB proteins are also involved in tapetum development. In Arabidopsis, AtMYB33 is involved in tapetum differentiation and AtMYB103/80 can participate in the late stage of tapetum development [39,40,41]. In rice, two rice MYB proteins, OsCSA2 and OsCSA, control male sterility by affecting the sugar partitioning in rice anthers [42,43]. Therefore, we searched for differentially expressed TFs in our transcriptome data.
We first conducted expression profiling of all transcription factors. A total of 1192 transcription factors (TFs) were identified and then classified into 25 families (Figure 6A). Here, we put the TF families with relatively few members under the category of ‘others’ (Figure 6A). In addition, a heatmap representation of TF expression profiles and hierarchical clustering analysis is shown in Figure 6B. The heatmap shows that 1192 TFs were expressed differentially among the Lo 4 cm, Lo 6 cm and Lo 8 cm stages (Figure 6B).
Previous studies have shown that MADS-box transcription factors (such as AG) and MYB transcription factors (such as AtMYB33, AtMYB103/80, OsCSA) play an important role in the development of anthers [37,44]. In addition, bHLH transcription factors, DISFUNCTIONAL TAPETUM1(DYT1) and ABORTED MICROSPORES (AMS) from Arabidopsis are also important for tapetal and microspore development [9,45,46]. The tapetum differentiation and callose degradation processes are abnormal in the dyt1 mutant [47]. In rice, TAPETUM DEGENERATION RETARDATION (TDR), which is the ortholog of AMS, can interact with ETERNAL TAPETUM1 (EAT1) and regulates tapetum development [48,49]. It has been shown that the PHD-finger genes, such as the Arabidopsis MALE STERILITY 1 (MS1) and its rice ortholog PERSISTENT TAPETAL CELL1 (OsPTC1) gene, are critical for pollen formation in plants [38,50]. The homology of these genes (AtAMS and OsTDR; AtMS1 and OsPTC1) suggests a conserved switch in the regulation of anther development in both dicots and monocots. Therefore, we shed light on these four transcription factor families, MADS-box, MYB, bHLH and PHD, in the transcriptome data. Additionally, the NAC transcription factors, NAC SECONDARY WALL THICKENING PROMOTING FACTOR 1 and 2 (NST1 and NST2), are required for the anther dehiscence in Arabidopsis [51,52]. Similarly, in wheat, 11 TabZIP genes have been identified to be predominantly expressed in anthers, based on whole wheat genome sequencing [53]. These results suggest that members of the NAC family and bZIP family may also be involved in the process of tapetum and microspores development. Hence, we also took the NAC and bZIP families into account to find differentially expressed transcription factors.
From Figure 6A, these six families NAC, MYB, bZIP, bHLH, MADS and PHD, were found to have 23, 61, 52, 48, 39 and 36 members (259/1192), respectively, accounting for 1.93%, 5.12%, 4.36%, 4.03%, 3.27% and 3.02% of all TFs identified. Next, we used BLAST to obtain annotation information of these 259 transcription factors and combined information with the Swiss-Prot and Pfam databases for further screening, to obtain 45 candidate TFs to focus our further analyses on (Figure 7). The heatmap analysis showed the number of members of each family and their differential expression at different stages of anther development, respectively (Figure 7). Among the six TF families, the MADS-box family accounted for the largest percentage, containing nine TFs (Figure 7E). Two AGAMOUS-LIKE genes, AGL14 and AGL19, with higher expression levels may function at the late stages of pollen formation (Figure 7E). The NAC, MYB, bZIP and bHLH families had the same number of members, all of which had eight TFs (Figure 7A–D). In addition, the PHD family had only four TFs (Figure 7F). Among them, two MYB33 genes participated in the tapetum differentiation (Figure 7B). In rice, ABSCISIC ACID-INSENSITIVE 5 (OsABI5), which belong to the bZIP family, are involved in the regulation of pollen maturation. The rice abi5 mutant exhibits abnormal pollens causing low fertility rates [54]. We found three DEGs of ABI5 in the transcriptome and they may regulate pollen development in lily (Figure 7C). Considering their differential expression during anther development, it is possible that these 45 transcription factors are involved in the tapetum degradation and microspore development in the lily cultivar ‘Siberia’.

3.6. Identification of DEGs Related to Tapetum Degradation and Microspore Development

Based on the gene annotations, we checked a total of 6896 DEGs one by one, as well as the expression pattern of their homologous genes in anther using Arabidopsis ‘eFP browser’. We finally identified 59 DEGs related to tapetum degradation and microspore development and divided them into 5 groups—cell-cycle- and division-related genes, cell-wall-formation-related genes, transcription-factor-related genes, leucine-rich repeat receptor-like kinase (LRR-RLK)-related genes and plant-hormone-synthesis- and signal-transduction-related genes (Figure 8A–E). All the expression patterns and annotations of the DEGs are shown in the heatmap in Figure 8.
The tapetum plays a significant role in the development of anthers. When the microspores complete meiosis, the tapetum degrades the callose to release the microspores and, when the pollen grains enter the mitosis, the tapetum layer is further degraded to form a pollen coat [5,41]. Therefore, the degradation of the tapetum is accompanied by cell division and expansion; four DEGs were thus selected from the transcriptome database. Two cyclin-related DEGs were expressed more in the early stages, while the rest of the genes were expressed more in the Lo 6 cm and Lo 8 cm stages (Figure 8A). These genes may have an effect on mitosis and meiosis during microspore and tapetum development.
Except for the fact that the tapetum is generally involved in the cell division process, it has a vital effect on the formation of the pollen wall. The tapetum can synthesize the precursor material of the outer wall of pollen and release the required material for the inner wall of pollen after being degraded [55,56,57]. Given the relationship between the tapetum and the pollen wall, we identified seven DEGs that were related to cell-wall synthesis (Figure 8B). Three genes were highly expressed in the Lo 4 cm stage and may be involved in the process of precursors synthesis of the outer pollen wall. The other five genes were highly expressed in a later stage, which may be involved in the formation of the inner pollen wall (Figure 8B).
Previous research showed that transcription factors and leucine-rich repeat receptor-like kinases (LRR-RLKs) play important roles in anther development [35,58]. Thus, we identified 11 transcription factors (Figure 8C). The heat map shows that, except for AGL19, the rest of them were up-regulated in the period from Lo 4 cm to Lo 6 cm and could be used for subsequent analyses. It has been reported that these genes, LRR-RLKs, EMS1/EXS, BAM1/2, SERK1/2 and RPK2, can affect tapetum formation and degradation [59,60,61]. For example, the bam1/2 double mutant lacks endothelium and tapetum [62,63]. In the process of searching for DEGs, we found five BAM1, which were able to encode LRR-RLK proteins (Figure 8D). All these genes were highly expressed in the Lo 4 cm and Lo 6 cm stages, suggesting that they were possibly involved in tapetum development (Figure 8D). Meanwhile, we also considered that tapetum degradation and microspore development may also be regulated by plant hormones [34,64,65], so we found seven genes with high expression levels in the Lo 4 cm and Lo 6 cm stages (Figure 8E). They were involved in the synthesis and signal transduction of IAA, ABA, JA and GA. For example, one DEG encodes NCED (9-cis-epoxycarotenoid dioxygenase), a rate-limiting enzyme in the ABA biosynthesis, and GID1C, an important receptor in the GA signaling pathway. These DEGs were candidates for tapetum degradation in lily.

3.7. Differential Gene Expression Validation by qRT-PCR Analysis

To further verify the transcriptomic analysis, 16 DEGs were selected for qRT-PCR validation. These genes included transcription factors related to anther development and callose degradation, genes related to microspore maturation and genes related to plant hormone signaling. The results showed that the expression patterns of these DEGs were similar to the changes in FPKM values in the transcriptome data (Figure 9), thus supporting that our transcriptomic data and analyses are of high quality.

4. Discussion

Pollen is produced from the anthers in the stamens and reaches the pistils for fertilization [66]. The development of pollen goes through the stages of archesporial cells, sporogenous cells, dyads, tetrads and microspores; then, it gradually develops into mature pollen grains [67]. Mature pollen grains play a significant role in sexual reproduction. However, in some ornamental plants, such as lily flowers, the presence of a large number of pollen grains may cause allergies in some people and hinder the breeding process [2]. The tapetum is a special cell layer, the closest one to the anther chamber. It can provide a series of enzymes and nutrients for microspore development and secretes β-1,3-glucanase to dissolve the callose for microspore releasing [68,69]. The microspores further develop into mature pollen. Moreover, the degradation of the tapetum itself forms a pollen coating, which is necessary for pollen germination [70]. Therefore, the tapetum is crucial for pollen development and maturation.
In this study, we used a popular lily cultivar, ‘Siberia’, as plant material. We continuously observed the development of lily anthers and found that anthers from 4 cm, 6 cm and 8 cm (Lo 4 cm, Lo 6 cm and Lo 8 cm) buds were the key stages for tapetum degradation and microspore development. Then, anthers from Lo 4 cm, Lo 6 cm and Lo 8 cm were selected for subsequent transcriptome sequencing. We finally identified several genes possibly related to tapetum degradation and microspore development by a series of comparative analyses. By identifying the key genes regulating tapetum development and degradation, we can better understand the process of pollen development and contribute to the production of new lily varieties without pollen contamination by genetic engineering.
We first clarified the development process of the tapetum and microspores by paraffin-sectioning the anthers of lily buds at different lengths (Figure 1). It is reported that Sanders et al. divided the development of Arabidopsis stamens into 14 stages [67]. As Sanders et al. described, microspore mother cells differentiated and formed in the stages 1–5. The tetrads formed and were surrounded by callose in the stages 6–7, which corresponds to our observation of Lo 4 cm (Figure 1C). The microspores are released from the tetrad and develop further in the 8–11 period, which corresponds to Lo 6 cm (Figure 1D). Then, the pollen grains gradually mature with the tapetum degradation and are released from the anther in the 12–14 period, which is similar to what we observed in the Lo 8 cm stage (Figure 1E). Therefore, we thought that there was a certain similarity in tapetum and micropore development between Lily and Arabidopsis.
The degradation of tapetum and the development of microspores are regulated by strict cell fate determination and morphogenesis mechanisms [71,72]. Consequently, we found four genes, which were possibly involved in the cell cycle and division (Figure 8A). These genes may have an influence on the fate of tapetum cells and microspores. In Arabidopsis, TAPETUM DETERMINANT1 (TPD1) is a potential small protein-ligand of the EXCESS MICROSPOROCYTES1 (EMS1) LRR-RLK and the TPD1-EMS1 signal module is specifically required for anther cell differentiation during sexual reproduction [73,74]. Overexpression of TPD1 increases the number of carpel cells and delays the degradation of tapetum cells [74]. Its rice homolog, OsTPD1-like (OsTDL1A), binds to the receptor kinase MULTIPLE SPOROCYTES1 (MSP1) and is also required for cell proliferation in early pollen development [75,76]. In addition, studies have shown that changes in the expression level of TPD1 not only affected the endogenous auxin response but also enhanced the expression of cyclin genes CYCD3;3 and CYCA2;3, suggesting that TPD1 functions by affecting the auxin response pathway and cell-cycle process [73]. Similarly, we found one DEG, annotated as TPD1, and two cyclin genes, CYCA2;1 and CYCD4;1, with different expression levels at the Lo 4 cm, Lo 6 cm and Lo 8 cm stages (Figure 8A). This suggests that they may have a conserved function in lily tapetum and microspores development (Figure 8A). Apart from the TPD1 and cyclin genes, two DEGs, which were respectively aligned to EXO70A1 and DAR1 (DA1-RELATED 1), were identified to be expressed more in the Lo 6 cm and Lo 8 cm stages (Figure 8A). Previously, EXO70A1 has been shown to be required to control pollen hydration and pollen tube development and DAR1 can regulate endoreduplication in leaf epidermal tissue [77,78]. Through the eFP browser, we observed that DAR1 (AT4G36860) in Arabidopsis was also expressed differentially during anther development (Figure S2A). Therefore, these genes are possibly involved in cell fate determination in lily anthers.
The microspore mother cell undergoes meiosis to form microspores, accompanied by vacuolation of tapetum cells [79]. When meiosis is complete, tapetum cells can synthesize sporopollenin and other substances to provide the material basis for the formation of the microspore pollen wall. After the microspores are released, tapetum degradation can provide material and energy for the formation of the pollen wall [41,57]. Considering that the development of tapetum and microspores is closely related to the process of cell-wall formation, we identified seven enzymes possibly involved in these processes (Figure 8B). These genes encode cellulose synthase (CSLD2, CESA1 and CESA6), xylose isomerase (AT5G57655), pectin-esterase (PE11), sucrose synthase (SUS1) and beta-1,3-glucanase (BGL2), respectively, and have been shown to participate in the cell-wall formation process (Figure 8B) [80,81,82,83]. These seven genes were differentially expressed in the Lo 4 cm, Lo 6 cm and Lo 8 cm stages of lily bud development. We observed the expression levels of these genes in the eFP browser and found that, except for beta-1,3-glucanase-related genes with low expression levels, the other genes were differentially expressed in the anther development of Arabidopsis (Figure 8B and Figure S2B). Hence, we thought that these six enzyme-encoding genes had a potential impact on tapetum and microspores development.
LRR-RLKs (leucine-rich repeat receptor-like kinases) are key components in many plant signal-transduction processes and play an important role in plant growth and development, resistance to adversity stress, hormone signal transduction, etc. [84,85,86]. Notably, several LRR-RLKs work together to control anther-cell differentiation. EXCESS MICROSPOROCYTES1 (EMS1)/EXTRA SPOROGENOUS CELLS (EXS) and its small protein-ligand TAPETUM DETERMINANT1 (TPD1), SOMATIC EMBRYOGENESIS RECEPTOR-LIKE KINASES1/2 (SERK1/2) and BARELY ANY MERISTEM1/2 (BAM1/2) directly have a significant impact on the formation of the endothecium, middle layer, tapetum and other anther structures [87,88,89,90]. Besides, a novel factor, RECEPTOR-LIKE PROTEIN KINASE (RPK2), can control tapetum development in the late stage of the anther in Arabidopsis [91]. In rice, MULTIPLE SPOROCYTES1 (OsMSP1), which encodes an LRR-RLK, is closely similar to AtEXS/EMS1 in structure and function [92]. In all these mutants, the lack of these kinases results in the repression of tapetum formation and development. In our database, we found five DEGs that were annotated as BAM1 but no other LRR-RLKs were found (Figure 8D). They were expressed more in the early stage of the lily anther, which may affect the development of tapetum in lily (Figure 8D). This result also reflects that the regulation mechanism of anther development is somewhat conserved.
The physiological effects of plant hormones are complex and diverse. The dynamic changes in endogenous hormone content and the hormone signal transduction are involved in the various processes of plant growth and development [93,94,95]. Previous research has demonstrated that the development of anthers is regulated by plant hormones. Hirano K. et al. analyzed the endogenous content of plant hormones in different tissues of rice and found that the contents of IAA and GA in mature anthers were significantly higher than those in other tissues [64]. In the Petunia hybrida L., it was demonstrated that ABA and IAA were involved in the programmed cell death (PCD) of microsporocytes at the meiosis stage [96]. In addition, JA is also involved in the reproductive process of plants and the phenotype of the JA biosynthetic mutant is profound sporophytic male sterility [97]. Consequently, we paid attention to these four hormones, ABA, GA, IAA and JA, and found seven related DEGs. These genes were either related to hormone synthesis or as hormone receptors involved in hormone signal transduction and may participate in tapetum degradation and micropore development. As previous research has shown, AUX/IAA and ARF (AUXIN RESPONSE FACTOR) genes play a crucial role in pollen wall formation, tapetum development and anther dehiscence [64,98,99,100,101]. We found two DEGs, IAA8 and ARF23, whose expression levels were high in the Lo 4 cm period, indicating that they may play a regulatory role in the early stage of tapetum and microspore development (Figure 8E). NCED is a key rate-limiting enzyme in the ABA synthesis pathway and PYL8 is an important ABA receptor [102,103]. Two DEGs were identified as NCED and PYL8 in lily, respectively, and they were all highly expressed in the Lo 6 cm stage (Figure 8E). There are five NCED homologous genes in Arabidopsis, namely, NCED2, NCED3, NCED5, NCED6 and NCED9, and their expression levels in anthers were low in the eFP browser (Figure S2C). Conversely, PYL8, which corresponds to AT5G53160 in Arabidopsis thaliana, was expressed differently at different stages of anther development (Figure S2C), indicating that PYL8 may regulate the development of tapetum and microspores through the ABA signal transduction pathway. In addition, we also found two DEGs, which were respectively annotated as GID1C and JAR1 (Figure 8E); GID1C is an important receptor of the GA signaling pathway and JAR1 is involved in JA biosynthesis [104,105]. These DEGs were highly expressed in the Lo 6 cm period (Figure 8E). Through the eFP-browser tool, GID1C (AT5G27320) showed no differences in the expression of anthers development, but JAR1 (AT2G46370) was significantly less expressed in mature pollen grains than in other stages (Figure S2C), indicating that JAR1 may have an impact on the development of anthers by affecting the JA synthesis process. In fact, it is still unclear how plant hormones regulate tapetum degradation and microspore development and further research is required to clarify the underlying mechanisms.
In conclusion, we used ‘Siberia’, an important lily cut flower variety, as the experimental material and determined that the period from Lo 4 cm to Lo 8 cm is the key period for the degradation of anther tapetum and the development of microspores. The transcriptome analysis revealed that several DEGs, involved in cell division and expansion, cell-wall morphogenesis, transcription factors, LRR-RLKs, plant hormone biosynthesis and transduction, possibly regulated tapetum degradation and micropore development in lilies. The quantitative real-time PCR of the DEGs in three stages, Lo 4 cm, Lo 6 cm and Lo 8 cm confirmed the transcriptome data. This study helps us better understand the development process of lily anthers and provides valuable candidate genes to be further studied that may make it possible to cultivate new varieties without pollen contamination.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/genes13020366/s1, Figure S1: Biological replicates of three anther stages showed by the first two PCs which explained over 94% of variance among samples, Figure S2: The ‘eFP-seq’ browser of genes related to the tapetum degradation and microspore development. The color scale from yellow to red represents expression levels from low to high, Table S1: Primers used in this study, Table S2: KEGG enrichment analysis of clusted genes, Table S3: DEGs of LO6cm vs. LO4cm, Table S4: DEGs of LO8cm vs LO6cm, Supplementary Method S1.

Author Contributions

Conceptualization, J.S., W.J. and Y.Z.; data curation, J.S. and Y.Z.; formal analysis, J.S., Y.X. and Y.Z.; funding acquisition, J.S. and Y.Z.; methodology, J.S., W.J. and Y.Z.; supervision, Y.Z.; writing—original draft, J.S., W.J. and Y.Z.; writing—review & editing, J.S., Y.X. and Y.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by National Natural Science Foundation of China (grants 31902047 to J.S. and 32102536 to Y.Z.) and Natural Science Foundation of Anhui Province (grant 2008085MC79 to J.S.).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The transcriptome data are available at the National Genomics Data Center of China (project PRJCA008073; data access link: https://ngdc.cncb.ac.cn/gsa/browse/CRA005948, accessed on 23 June 2021). The codes used for data analysis methods are reported in Supplementary Method S1.

Acknowledgments

We thank Carina Carianopol (Biology Science at University of Toronto; e-mail, [email protected]) for help in revising the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Sui, J.; Cao, X.; Yi, M.; Wu, J.; He, J. Isolation and characterization of LoAMS gene in anther development of lily (Lilium oriental hybrids). N. Z. J. Crop. Hortic. Sci. 2020, 48, 257–269. [Google Scholar] [CrossRef]
  2. Wang, X.; Wu, Z.; Wang, L.; Wu, M.; Zhang, D.; Fang, W.; Chen, F.; Teng, N. Cytological and molecular characteristics of pollen abortion in lily with dysplastic tapetum. Hortic. Plant J. 2019, 5, 281–294. [Google Scholar] [CrossRef]
  3. Feng, J.; Wu, Z.; Wang, X.; Zhang, Y.; Teng, N. Analysis of Pollen Allergens in Lily by Transcriptome and Proteome Data. Int. J. Mol. Sci. 2019, 20, 5892. [Google Scholar] [CrossRef] [PubMed]
  4. Sui, J.; He, J.; Wu, J.; Gong, B.; Cao, X.; Seng, S.; Wu, Z.; Wu, C.; Liu, C.; Yi, M. Characterization and functional analysis of transcription factor lomyb80 related to anther development in lily (Lilium oriental hybrids). J. Plant Growth Regul. 2015, 34, 545–557. [Google Scholar] [CrossRef]
  5. Goldberg, R.B.; Beals, T.P.; Sanders, P.M. Anther development: Basic principles and practical applications. Plant Cell 1993, 5, 1217. [Google Scholar]
  6. Zhang, S.; Fang, Z.; Zhu, J.; Gao, J.; Yang, Z. OsMYB103 is required for rice anther development by regulating tapetum development and exine formation. Chin. Sci. Bull. 2010, 55, 3288–3297. [Google Scholar] [CrossRef]
  7. Mariani, C.; De Beuckeleer, M.; Truettner, J.; Leemans, J.; Goldberg, R.B. Induction of male sterility in plants by a chimaeric ribonuclease gene. Nature 1990, 347, 737–741. [Google Scholar] [CrossRef]
  8. Millar, A.A.; Gubler, F. The Arabidopsis GAMYB-like genes, MYB33 and MYB65, are microRNA-regulated genes that redundantly facilitate anther development. Plant Cell 2005, 17, 705–721. [Google Scholar] [CrossRef]
  9. Zhang, W.; Sun, Y.; Timofejeva, L.; Chen, C.; Grossniklaus, U.; Ma, H. Regulation of Arabidopsis tapetum development and function by DYSFUNCTIONAL TAPETUM1 (DYT1) encoding a putative bHLH transcription factor. Development 2006, 133, 3085–3095. [Google Scholar] [CrossRef]
  10. Zhu, J.; Chen, H.; Li, H.; Gao, J.F.; Jiang, H.; Wang, C.; Guan, Y.F.; Yang, Z.N. Defective in Tapetal development and function 1 is essential for anther development and tapetal function for microspore maturation in Arabidopsis. Plant J. 2008, 55, 266–277. [Google Scholar] [CrossRef]
  11. Sorensen, A.M.; Kröber, S.; Unte, U.S.; Huijser, P.; Dekker, K.; Saedler, H. The Arabidopsis ABORTED MICROSPORES (AMS) gene encodes a MYC class transcription factor. Plant J. 2003, 33, 413–423. [Google Scholar] [CrossRef] [PubMed]
  12. Xu, J.; Yang, C.; Yuan, Z.; Zhang, D.; Gondwe, M.Y.; Ding, Z.; Liang, W.; Zhang, D.; Wilson, Z.A. The ABORTED MICROSPORES regulatory network is required for postmeiotic male reproductive development in Arabidopsis thaliana. Plant Cell 2010, 22, 91–107. [Google Scholar] [CrossRef]
  13. Ferguson, A.C.; Pearce, S.; Band, L.R.; Yang, C.; Ferjentsikova, I.; King, J.; Yuan, Z.; Zhang, D.; Wilson, Z.A. Biphasic regulation of the transcription factor ABORTED MICROSPORES (AMS) is essential for tapetum and pollen development in Arabidopsis. New Phytol. 2017, 213, 778–790. [Google Scholar] [CrossRef] [PubMed]
  14. Zhang, Z.B.; Zhu, J.; Gao, J.F.; Wang, C.; Li, H.; Li, H.; Zhang, H.Q.; Zhang, S.; Wang, D.M.; Wang, Q.X. Transcription factor AtMYB103 is required for anther development by regulating tapetum development, callose dissolution and exine formation in Arabidopsis. Plant J. 2007, 52, 528–538. [Google Scholar] [CrossRef] [PubMed]
  15. Zhu, J.; Zhang, G.; Chang, Y.; Li, X.; Yang, J.; Huang, X.; Yu, Q.; Chen, H.; Wu, T.; Yang, Z. AtMYB103 is a crucial regulator of several pathways affecting Arabidopsis anther development. Sci. China Life Sci. 2010, 53, 1112–1122. [Google Scholar] [CrossRef] [PubMed]
  16. Wilson, Z.A.; Morroll, S.M.; Dawson, J.; Swarup, R.; Tighe, P.J. The Arabidopsis MALE STERILITY1 (MS1) gene is a transcriptional regulator of male gametogenesis, with homology to the PHD-finger family of transcription factors. Plant J. 2001, 28, 27–39. [Google Scholar] [CrossRef]
  17. Yang, C.; Vizcay-Barrena, G.; Conner, K.; Wilson, Z.A. MALE STERILITY1 is required for tapetal development and pollen wall biosynthesis. Plant Cell 2007, 19, 3530–3548. [Google Scholar] [CrossRef]
  18. Huang, G.; Yi, Q.; Tong, L.; Hui, Y.; Aide, W.; Dongmei, T. Comparative transcriptome analysis of actinidia arguta fruits reveals the involvement of various transcription factors in ripening. Hortic. Plant J. 2018, 4, 35–42. [Google Scholar] [CrossRef]
  19. Salmela, L.; Rivals, E. LoRDEC: Accurate and efficient long read error correction. Bioinformatics 2014, 30, 3506–3514. [Google Scholar] [CrossRef]
  20. Wu, T.D.; Watanabe, C.K. GMAP: A genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics 2005, 21, 1859–1875. [Google Scholar] [CrossRef]
  21. Fu, L.; Niu, B.; Zhu, Z.; Wu, S.; Li, W. CD-HIT: Accelerated for clustering the next-generation sequencing data. Bioinformatics 2012, 28, 3150–3152. [Google Scholar] [CrossRef] [PubMed]
  22. Zheng, Y.; Jiao, C.; Sun, H.; Rosli, H.; Pombo, M.A.; Zhang, P.; Banf, M.; Dai, X.; Martin, G.B.; Giovannoni, J.J. iTAK: A program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases. Mol. Plant 2016, 9, 1667–1670. [Google Scholar] [CrossRef] [PubMed]
  23. Xie, C.; Mao, X.; Huang, J.; Ding, Y.; Wu, J.; Dong, S.; Kong, L.; Gao, G.; Li, C.-Y.; Wei, L. KOBAS 2.0: A web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011, 39, W316–W322. [Google Scholar] [CrossRef] [PubMed]
  24. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef]
  25. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef]
  26. Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323. [Google Scholar] [CrossRef]
  27. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef]
  28. Sullivan, A.; Purohit, P.K.; Freese, N.H.; Pasha, A.; Esteban, E.; Waese, J.; Wu, A.; Chen, M.; Chin, C.Y.; Song, R. An ‘eFP-Seq Browser’for visualizing and exploring RNA sequencing data. Plant J. 2019, 100, 641–654. [Google Scholar] [CrossRef]
  29. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  30. Lutgen, D.; Ritter, R.; Olsen, R.A.; Schielzeth, H.; Gruselius, J.; Ewels, P.; García, J.T.; Shirihai, H.; Schweizer, M.; Suh, A. Linked-read sequencing enables haplotype-resolved resequencing at population scale. Mol. Ecol. Resour. 2020, 20, 1311–1322. [Google Scholar] [CrossRef]
  31. Consortium, G.O. The Gene Ontology (GO) database and informatics resource. Nucleic Acids Res. 2004, 32, D258–D261. [Google Scholar] [CrossRef] [PubMed]
  32. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [Google Scholar] [CrossRef] [PubMed]
  33. Martínez-Reyes, I.; Diebold, L.P.; Kong, H.; Schieber, M.; Huang, H.; Hensley, C.T.; Mehta, M.M.; Wang, T.; Santos, J.H.; Woychik, R.; et al. TCA Cycle and Mitochondrial Membrane Potential Are Necessary for Diverse Biological Functions. Mol. Cell 2016, 61, 199–209. [Google Scholar] [CrossRef] [PubMed]
  34. Weigt, D.; Niemann, J.; Siatkowski, I.; Zyprych-Walczak, J.; Olejnik, P.; Kurasiak-Popowska, D. Effect of zearalenone and hormone regulators on microspore embryogenesis in anther culture of wheat. Plants 2019, 8, 487. [Google Scholar] [CrossRef] [PubMed]
  35. Verma, N. Transcriptional regulation of anther development in Arabidopsis. Gene 2019, 689, 202–209. [Google Scholar] [CrossRef] [PubMed]
  36. Bowman, J.L.; Drews, G.N.; Meyerowitz, E.M. Expression of the Arabidopsis floral homeotic gene AGAMOUS is restricted to specific cell types late in flower development. Plant Cell 1991, 3, 749–758. [Google Scholar] [PubMed]
  37. Sieburth, L.E.; Meyerowitz, E.M. Molecular dissection of the AGAMOUS control region shows that cis elements for spatial regulation are located intragenically. Plant Cell 1997, 9, 355–365. [Google Scholar]
  38. Ito, T.; Ng, K.-H.; Lim, T.-S.; Yu, H.; Meyerowitz, E.M. The homeotic protein AGAMOUS controls late stamen development by regulating a jasmonate biosynthetic gene in Arabidopsis. Plant Cell 2007, 19, 3516–3529. [Google Scholar] [CrossRef]
  39. Higginson, T.; Li, S.F.; Parish, R.W. AtMYB103 regulates tapetum and trichome development in Arabidopsis thaliana. Plant J. 2003, 35, 177–192. [Google Scholar] [CrossRef]
  40. Alonso-Peral, M.M.; Li, J.; Li, Y.; Allen, R.S.; Schnippenkoetter, W.; Ohms, S.; White, R.G.; Millar, A.A. The microRNA159-regulated GAMYB-like genes inhibit growth and promote programmed cell death in Arabidopsis. Plant Physiol. 2010, 154, 757–771. [Google Scholar] [CrossRef]
  41. Liu, L.; Fan, X.-d. Tapetum: Regulation and role in sporopollenin biosynthesis in Arabidopsis. Plant Mol. Biol. 2013, 83, 165–175. [Google Scholar] [CrossRef] [PubMed]
  42. Wang, D.; Li, J.; Sun, L.; Hu, Y.; Yu, J.; Wang, C.; Zhang, F.; Hou, H.; Liang, W.; Zhang, D. Two rice MYB transcription factors maintain male fertility in response to photoperiod by modulating sugar partitioning. New Phytol. 2021, 231, 1612–1629. [Google Scholar] [CrossRef] [PubMed]
  43. Zhang, H.; Liang, W.; Yang, X.; Luo, X.; Jiang, N.; Ma, H.; Zhang, D. Carbon starved anther encodes a MYB domain protein that regulates sugar partitioning required for rice pollen development. Plant Cell 2010, 22, 672–689. [Google Scholar] [CrossRef] [PubMed]
  44. Parish, R.W.; Li, S.F. Death of a tapetum: A programme of developmental altruism. Plant Sci. 2010, 178, 73–89. [Google Scholar] [CrossRef]
  45. Xu, J.; Ding, Z.; Vizcay-Barrena, G.; Shi, J.; Liang, W.; Yuan, Z.; Werck-Reichhart, D.; Schreiber, L.; Wilson, Z.A.; Zhang, D. ABORTED MICROSPORES acts as a master regulator of pollen wall formation in Arabidopsis. Plant Cell 2014, 26, 1544–1556. [Google Scholar] [CrossRef]
  46. Guo, J.; Liu, C.; Wang, P.; Cheng, Q.; Sun, L.; Yang, W.; Shen, H. The aborted microspores (AMS)-like gene is required for anther and microspore development in pepper (Capsicum annuum L.). Int. J. Mol. Sci. 2018, 19, 1341. [Google Scholar] [CrossRef]
  47. Fiorio, M.; Gambarin, M.; Valente, E.M.; Liberini, P.; Loi, M.; Cossu, G.; Moretto, G.; Bhatia, K.P.; Defazio, G.; Aglioti, S.M. Defective temporal processing of sensory stimuli in DYT1 mutation carriers: A new endophenotype of dystonia? Brain 2007, 130, 134–142. [Google Scholar] [CrossRef]
  48. Li, N.; Zhang, D.S.; Liu, H.S.; Yin, C.S.; Li, X.X.; Liang, W.Q.; Yuan, Z.; Xu, B.; Chu, H.W.; Wang, J.; et al. The rice tapetum degeneration retardation gene is required for tapetum degradation and anther development. Plant Cell 2006, 18, 2999–3014. [Google Scholar] [CrossRef]
  49. Niu, N.; Liang, W.; Yang, X.; Jin, W.; Wilson, Z.A.; Hu, J.; Zhang, D. EAT1 promotes tapetal cell death by regulating aspartic proteases during male reproductive development in rice. Nat. Commun. 2013, 4, 1445. [Google Scholar] [CrossRef]
  50. Li, H.; Yuan, Z.; Vizcay-Barrena, G.; Yang, C.; Liang, W.; Zong, J.; Wilson, Z.A.; Zhang, D. PERSISTENT TAPETAL CELL1 encodes a PHD-finger protein that is required for tapetal cell death and pollen development in rice. Plant Physiol. 2011, 156, 615–630. [Google Scholar] [CrossRef]
  51. Mitsuda, N.; Seki, M.; Shinozaki, K.; Ohme-Takagi, M. The NAC transcription factors NST1 and NST2 of Arabidopsis regulate secondary wall thickenings and are required for anther dehiscence. Plant Cell 2005, 17, 2993–3006. [Google Scholar] [CrossRef] [PubMed]
  52. Zhang, Q.; Luo, F.; Zhong, Y.; He, J.; Li, L. Modulation of NAC transcription factor NST1 activity by XYLEM NAC DOMAIN1 regulates secondary cell wall formation in Arabidopsis. J. Exp. Bot. 2020, 71, 1449–1458. [Google Scholar] [CrossRef]
  53. Li, X.; Gao, S.; Tang, Y.; Li, L.; Zhang, F.; Feng, B.; Fang, Z.; Ma, L.; Zhao, C. Genome-wide identification and evolutionary analyses of bZIP transcription factors in wheat and its relatives and expression profiles of anther development related TabZIP genes. BMC Genom. 2015, 16, 976. [Google Scholar] [CrossRef] [PubMed]
  54. Zou, M.; Guan, Y.; Ren, H.; Zhang, F.; Chen, F. A bZIP transcription factor, OsABI5, is involved in rice fertility and stress tolerance. Plant Mol. Biol. 2008, 66, 675–683. [Google Scholar] [CrossRef] [PubMed]
  55. Vizcay-Barrena, G.; Wilson, Z.A. Altered tapetal PCD and pollen wall development in the Arabidopsis ms1 mutant. J. Exp. Bot. 2006, 57, 2709–2717. [Google Scholar] [CrossRef] [PubMed]
  56. Jiang, J.; Zhang, Z.; Cao, J. Pollen wall development: The associated enzymes and metabolic pathways. Plant Biol. 2013, 15, 249–263. [Google Scholar] [CrossRef] [PubMed]
  57. Li, D.-D.; Xue, J.-S.; Zhu, J.; Yang, Z.-N. Gene regulatory network for tapetum development in Arabidopsis thaliana. Front. Plant Sci. 2017, 8, 1559. [Google Scholar] [CrossRef] [PubMed]
  58. Wu, Y.; Xun, Q.; Guo, Y.; Zhang, J.; Cheng, K.; Shi, T.; He, K.; Hou, S.; Gou, X.; Li, J. Genome-wide expression pattern analyses of the Arabidopsis leucine-rich repeat receptor-like kinases. Mol. Plant 2016, 9, 289–300. [Google Scholar] [CrossRef]
  59. Ge, X.; Chang, F.; Ma, H. Signaling and transcriptional control of reproductive development in Arabidopsis. Curr. Biol. 2010, 20, R988–R997. [Google Scholar] [CrossRef]
  60. Yang, L.; Qian, X.; Chen, M.; Fei, Q.; Meyers, B.C.; Liang, W.; Zhang, D. Regulatory role of a receptor-like kinase in specifying anther cell identity. Plant Physiol. 2016, 171, 2085–2100. [Google Scholar] [CrossRef]
  61. Cai, W.; Zhang, D. The role of receptor-like kinases in regulating plant male reproduction. Plant Reprod. 2018, 31, 77–87. [Google Scholar] [CrossRef] [PubMed]
  62. Szakmary, A.; Cox, D.N.; Wang, Z.; Lin, H. Regulatory relationship among piwi, pumilio, and bag-of-marbles in Drosophila germline stem cell self-renewal and differentiation. Curr. Biol. 2005, 15, 171–178. [Google Scholar] [CrossRef] [PubMed]
  63. Hart, E.M.; Gupta, M.; Wühr, M.; Silhavy, T.J. The synthetic phenotype of Δ bamB Δ bamE double mutants results from a lethal jamming of the bam complex by the lipoprotein RcsF. MBio 2019, 10, e00662-19. [Google Scholar] [CrossRef] [PubMed]
  64. Hirano, K.; Aya, K.; Hobo, T.; Sakakibara, H.; Kojima, M.; Shim, R.A.; Hasegawa, Y.; Ueguchi-Tanaka, M.; Matsuoka, M. Comprehensive transcriptome analysis of phytohormone biosynthesis and signaling genes in microspore/pollen and tapetum of rice. Plant Cell Physiol. 2008, 49, 1429–1450. [Google Scholar] [CrossRef]
  65. Acosta, I.F.; Przybyl, M. Jasmonate signaling during Arabidopsis stamen maturation. Plant Cell Physiol. 2019, 60, 2648–2659. [Google Scholar] [CrossRef]
  66. Shivanna, K.R. Pollen Biology and Biotechnology; CRC Press: Boca Raton, FL, USA, 2019. [Google Scholar]
  67. Sanders, P.M.; Bui, A.Q.; Weterings, K.; McIntire, K.; Hsu, Y.-C.; Lee, P.Y.; Truong, M.T.; Beals, T.; Goldberg, R. Anther developmental defects in Arabidopsis thaliana male-sterile mutants. Sex. Plant Reprod. 1999, 11, 297–322. [Google Scholar] [CrossRef]
  68. Stieglitz, H.; Stern, H. Regulation of β-1, 3-glucanase activity in developing anthers of Lilium. Dev. Biol. 1973, 34, 169–173. [Google Scholar] [CrossRef]
  69. Wan, L.; Zha, W.; Cheng, X.; Liu, C.; Lv, L.; Liu, C.; Wang, Z.; Du, B.; Chen, R.; Zhu, L. A rice β-1, 3-glucanase gene Osg1 is required for callose degradation in pollen development. Planta 2011, 233, 309–323. [Google Scholar] [CrossRef]
  70. Pacini, E. Relationships between tapetum, loculus, and pollen during development. Int. J. Plant Sci. 2010, 171, 1–11. [Google Scholar] [CrossRef]
  71. Walbot, V.; Egger, R.L. Pre-meiotic anther development: Cell fate specification and differentiation. Annu. Rev. Plant Biol. 2016, 67, 365–395. [Google Scholar] [CrossRef]
  72. Zhang, D.; Yang, L. Specification of tapetum and microsporocyte cells within the anther. Curr. Opin. Plant Biol. 2014, 17, 49–55. [Google Scholar] [CrossRef] [PubMed]
  73. Huang, J.; Wijeratne, A.J.; Tang, C.; Zhang, T.; Fenelon, R.E.; Owen, H.A.; Zhao, D. Ectopic expression of TAPETUM DETERMINANT1 affects ovule development in Arabidopsis. J. Exp. Bot. 2016, 67, 1311–1326. [Google Scholar] [CrossRef] [PubMed]
  74. Yang, S.-L.; Jiang, L.; Puah, C.S.; Xie, L.-F.; Zhang, X.-Q.; Chen, L.-Q.; Yang, W.-C.; Ye, D. Overexpression of TAPETUM DETERMINANT1 alters the cell fates in the Arabidopsis carpel and tapetum via genetic interaction with excess microsporocytes1/extra sporogenous cells. Plant Physiol. 2005, 139, 186–191. [Google Scholar] [CrossRef] [PubMed]
  75. Hong, L.; Tang, D.; Shen, Y.; Hu, Q.; Wang, K.; Li, M.; Lu, T.; Cheng, Z. MIL2 (MICROSPORELESS2) regulates early cell differentiation in the rice anther. New Phytol. 2012, 196, 402–413. [Google Scholar] [CrossRef]
  76. Zhao, X.; de Palma, J.; Oane, R.; Gamuyao, R.; Luo, M.; Chaudhury, A.; Hervé, P.; Xue, Q.; Bennett, J. OsTDL1A binds to the LRR domain of rice receptor kinase MSP1, and is required to limit sporocyte numbers. Plant J. 2008, 54, 375–387. [Google Scholar] [CrossRef]
  77. Peng, Y.; Chen, L.; Lu, Y.; Wu, Y.; Dumenil, J.; Zhu, Z.; Bevan, M.W.; Li, Y. The ubiquitin receptors DA1, DAR1, and DAR2 redundantly regulate endoreduplication by modulating the stability of TCP14/15 in Arabidopsis. Plant Cell 2015, 27, 649–662. [Google Scholar] [CrossRef]
  78. Safavian, D.; Jamshed, M.; Sankaranarayanan, S.; Indriolo, E.; Samuel, M.A.; Goring, D.R. High humidity partially rescues the Arabidopsis thaliana exo70A1 stigmatic defect for accepting compatible pollen. Plant Reprod. 2014, 27, 121–127. [Google Scholar] [CrossRef]
  79. Wan, L.; Xia, X.; Hong, D.; Li, J.; Yang, G. Abnormal vacuolization of the tapetum during the tetrad stage is associated with male sterility in the recessive genic male sterile Brassica napus L. Line 9012A. J. Plant Biol. 2010, 53, 121–133. [Google Scholar] [CrossRef]
  80. Singh, P.; Rao, R.N.; Reddy, J.R.C.; Prasad, R.; Kotturu, S.K.; Ghosh, S.; Mukhopadhyay, S. PE11, a PE/PPE family protein of Mycobacterium tuberculosis is involved in cell wall remodeling and virulence. Sci. Rep. 2016, 6, 1–16. [Google Scholar]
  81. Jaquinod, M.; Villiers, F.; Kieffer-Jaquinod, S.; Hugouvieux, V.; Bruley, C.; Garin, J.; Bourguignon, J. A proteomics dissection of Arabidopsis thaliana vacuoles isolated from cell culture. Mol. Cell. Proteom. 2007, 6, 394–412. [Google Scholar] [CrossRef]
  82. Rodrıguez-Navarro, S.; Fischer, T.; Luo, M.-J.; Antúnez, O.; Brettschneider, S.; Lechner, J.; Pérez-Ortín, J.E.; Reed, R.; Hurt, E. Sus1, a functional component of the SAGA histone acetylase complex and the nuclear pore-associated mRNA export machinery. Cell 2004, 116, 75–86. [Google Scholar] [CrossRef]
  83. Sarthy, A.V.; McGonigal, T.; Coen, M.; Frost, D.J.; Meulbroek, J.A.; Goldman, R.C. Phenotype in Candida albicans of a disruption of the BGL2 gene encoding a 1, 3-β-glucosyltransferase. Microbiology 1997, 143, 367–376. [Google Scholar] [CrossRef] [PubMed]
  84. Furumizu, C.; Sawa, S. Insight into early diversification of leucine-rich repeat receptor-like kinases provided by the sequenced moss and hornwort genomes. Plant Mol. Biol. 2021, 107, 337–353. [Google Scholar] [CrossRef] [PubMed]
  85. Laohavisit, A.; Wakatake, T.; Ishihama, N.; Mulvey, H.; Takizawa, K.; Suzuki, T.; Shirasu, K. Quinone perception in plants via leucine-rich-repeat receptor-like kinases. Nature 2020, 587, 92–97. [Google Scholar] [CrossRef] [PubMed]
  86. Hosseini, S.; Schmidt, E.D.; Bakker, F.T. Leucine-rich repeat receptor-like kinase II phylogenetics reveals five main clades throughout the plant kingdom. Plant J. 2020, 103, 547–560. [Google Scholar] [CrossRef]
  87. Crook, A.D.; Willoughby, A.C.; Hazak, O.; Okuda, S.; Van Der Molen, K.R.; Soyars, C.L.; Cattaneo, P.; Clark, N.M.; Sozzani, R.; Hothorn, M. BAM1/2 receptor kinase signaling drives CLE peptide-mediated formative cell divisions in Arabidopsis roots. Proc. Natl. Acad. Sci. USA 2020, 117, 32750–32756. [Google Scholar] [CrossRef]
  88. Chen, W.; Lv, M.; Wang, Y.; Wang, P.-A.; Cui, Y.; Li, M.; Wang, R.; Gou, X.; Li, J. BES1 is activated by EMS1-TPD1-SERK1/2-mediated signaling to control tapetum development in Arabidopsis thaliana. Nat. Commun. 2019, 10, 1–12. [Google Scholar] [CrossRef]
  89. Zhao, D. Control of anther cell differentiation: A teamwork of receptor-like kinases. Sex. Plant Reprod. 2009, 22, 221–228. [Google Scholar] [CrossRef]
  90. Jia, G.; Liu, X.; Owen, H.A.; Zhao, D. Signaling of cell fate determination by the TPD1 small protein and EMS1 receptor kinase. Proc. Natl. Acad. Sci. USA 2008, 105, 2220–2225. [Google Scholar] [CrossRef]
  91. Mizuno, S.; Osakabe, Y.; Maruyama, K.; Ito, T.; Osakabe, K.; Sato, T.; Shinozaki, K.; Yamaguchi-Shinozaki, K. Receptor-like protein kinase 2 (RPK 2) is a novel factor controlling anther development in Arabidopsis thaliana. Plant J. 2007, 50, 751–766. [Google Scholar] [CrossRef]
  92. Nonomura, K.; Miyoshi, K.; Eiguchi, M.; Suzuki, T.; Miyao, A.; Hirochika, H.; Kurata, N. The MSP1 gene is necessary to restrict the number of cells entering into male and female sporogenesis and to initiate anther wall formation in rice. Plant Cell 2003, 15, 1728–1739. [Google Scholar] [CrossRef] [PubMed]
  93. Li, N.; Euring, D.; Cha, J.Y.; Lin, Z.; Lu, M.; Huang, L.-J.; Kim, W.Y. Plant hormone-mediated regulation of heat tolerance in response to global climate change. Front. Plant Sci. 2021, 11, 2318. [Google Scholar] [CrossRef] [PubMed]
  94. Wang, L.; Zou, Y.; Kaw, H.Y.; Wang, G.; Sun, H.; Cai, L.; Li, C.; Meng, L.-Y.; Li, D. Recent developments and emerging trends of mass spectrometric methods in plant hormone analysis: A review. Plant Methods 2020, 16, 1–17. [Google Scholar] [CrossRef] [PubMed]
  95. Blázquez, M.A.; Nelson, D.C.; Weijers, D. Evolution of plant hormone response pathways. Annu. Rev. Plant Biol. 2020, 71, 327–353. [Google Scholar] [CrossRef]
  96. Kovaleva, L.; Voronkov, A.; Zakharova, E.; Andreev, I. ABA and IAA control microsporogenesis in Petunia hybrida L. Protoplasma 2018, 255, 751–759. [Google Scholar] [CrossRef]
  97. Jewell, J.B.; Browse, J. Epidermal jasmonate perception is sufficient for all aspects of jasmonate-mediated male fertility in Arabidopsis. Plant J. 2016, 85, 634–647. [Google Scholar] [CrossRef]
  98. Bassa, C.; Mila, I.; Bouzayen, M.; Audran-Delalande, C. Phenotypes associated with down-regulation of Sl-IAA27 support functional diversity among Aux/IAA family members in tomato. Plant Cell Physiol. 2012, 53, 1583–1595. [Google Scholar] [CrossRef]
  99. Keller, M.; Schleiff, E.; Simm, S. miRNAs involved in transcriptome remodeling during pollen development and heat stress response in Solanum lycopersicum. Sci. Rep. 2020, 10, 1–13. [Google Scholar] [CrossRef]
  100. Richter, S.; Müller, L.M.; Stierhof, Y.-D.; Mayer, U.; Takada, N.; Kost, B.; Vieten, A.; Geldner, N.; Koncz, C.; Jürgens, G. Polarized cell growth in Arabidopsis requires endosomal recycling mediated by GBF1-related ARF exchange factors. Nat. Cell Biol. 2012, 14, 80–86. [Google Scholar] [CrossRef]
  101. Yao, X.; Chen, J.; Zhou, J.; Yu, H.; Ge, C.; Zhang, M.; Gao, X.; Dai, X.; Yang, Z.-N.; Zhao, Y. An essential role for miRNA167 in maternal control of embryonic and seed development. Plant Physiol. 2019, 180, 453–464. [Google Scholar] [CrossRef]
  102. Gavassi, M.A.; Silva, G.S.; da Silva, C.d.M.S.; Thompson, A.J.; Macleod, K.; Oliveira, P.M.R.; Cavalheiro, M.F.; Domingues, D.S.; Habermann, G. NCED expression is related to increased ABA biosynthesis and stomatal closure under aluminum stress. Environ. Exp. Bot. 2021, 185, 104404. [Google Scholar] [CrossRef]
  103. Qi, L.; Liu, S.; Li, C.; Fu, J.; Jing, Y.; Cheng, J.; Li, H.; Zhang, D.; Wang, X.; Dong, X. PHYTOCHROME-INTERACTING FACTORS interact with the ABA receptors PYL8 and PYL9 to orchestrate ABA signaling in darkness. Mol. Plant 2020, 13, 414–430. [Google Scholar] [CrossRef] [PubMed]
  104. Chen, H.; Wang, B.; Geng, S.; Arellano, C.; Chen, S.; Qu, R. Effects of overexpression of jasmonic acid biosynthesis genes on nicotine accumulation in tobacco. Plant Direct 2018, 2, e00036. [Google Scholar] [CrossRef] [PubMed]
  105. Vidhyasekaran, P. Gibberellin Signaling in Plant Innate Immunity. In Plant Hormone Signaling Systems in Plant Innate Immunity; Springer: Berlin/Heidelberg, Germany, 2015; pp. 383–401. [Google Scholar]
Genes 13 00366 g001 550
Figure 1. Morphological and histological analyses of lily flower buds. (A) Lily flower when pollen grain was released. (B) Lily flower buds and anthers at Lo 4 cm, Lo 6 cm and Lo 8 cm stages. The scale bar is 1 cm. (C) Anther sections at Lo 4 cm stage. The scale bar is 100 μm. (D) Anther sections at Lo 6 cm stage. The scale bar is 200 μm. (E) Anther sections at Lo 8 cm stage. The scale bar is 200 μm. T, tapetum; Tds, tetrads; N, nucleus; Msp, microspores; PG, pollen grain.
Figure 1. Morphological and histological analyses of lily flower buds. (A) Lily flower when pollen grain was released. (B) Lily flower buds and anthers at Lo 4 cm, Lo 6 cm and Lo 8 cm stages. The scale bar is 1 cm. (C) Anther sections at Lo 4 cm stage. The scale bar is 100 μm. (D) Anther sections at Lo 6 cm stage. The scale bar is 200 μm. (E) Anther sections at Lo 8 cm stage. The scale bar is 200 μm. T, tapetum; Tds, tetrads; N, nucleus; Msp, microspores; PG, pollen grain.
Genes 13 00366 g001
Genes 13 00366 g002 550
Figure 2. Summary of transcriptome sequencing and gene annotation. (A) Summary of transcripts before and after correction. (B) Length distribution map of unigenes. (C) Statistics of functional annotation in seven databases.
Figure 2. Summary of transcriptome sequencing and gene annotation. (A) Summary of transcripts before and after correction. (B) Length distribution map of unigenes. (C) Statistics of functional annotation in seven databases.
Genes 13 00366 g002
Genes 13 00366 g003 550
Figure 3. DEGs among Lo 4 cm, Lo 6 cm and Lo 8 cm. (A) Venn diagram showing the number of DEGs in different comparisons. (B) Numbers of DEGs in pairwise comparisons. (C) Heatmap diagrams showing the relative expression levels of DEGs. The color scale on the right from blue to red represents FPKM values from low to high. (D) The trend analysis of DEG expression fold-change divided into nine clusters based on expression patterns.
Figure 3. DEGs among Lo 4 cm, Lo 6 cm and Lo 8 cm. (A) Venn diagram showing the number of DEGs in different comparisons. (B) Numbers of DEGs in pairwise comparisons. (C) Heatmap diagrams showing the relative expression levels of DEGs. The color scale on the right from blue to red represents FPKM values from low to high. (D) The trend analysis of DEG expression fold-change divided into nine clusters based on expression patterns.
Genes 13 00366 g003
Genes 13 00366 g004 550
Figure 4. GO classification of DEGs in Lo 6 cm vs. Lo 4 cm and Lo 8 cm vs. Lo 6 cm stages. (A) Enhanced GO terms in Lo 6 cm vs. Lo 4 cm. (B) Enhanced GO terms in Lo 8 cm vs. Lo 6 cm. (C) Volcano chart showing the distribution of DEGs for each comparison. The x-axis represents the log2FC (FoldChange) value, the y-axis the −log10 p value and the dotted line represents the threshold line of the DEGs screening criteria. The red dots and blue dots represent up- and down-regulated genes, respectively, and genes with no significant differential expression are represented by black dots. BP, biological process; MF, molecular function; CC, cellular component.
Figure 4. GO classification of DEGs in Lo 6 cm vs. Lo 4 cm and Lo 8 cm vs. Lo 6 cm stages. (A) Enhanced GO terms in Lo 6 cm vs. Lo 4 cm. (B) Enhanced GO terms in Lo 8 cm vs. Lo 6 cm. (C) Volcano chart showing the distribution of DEGs for each comparison. The x-axis represents the log2FC (FoldChange) value, the y-axis the −log10 p value and the dotted line represents the threshold line of the DEGs screening criteria. The red dots and blue dots represent up- and down-regulated genes, respectively, and genes with no significant differential expression are represented by black dots. BP, biological process; MF, molecular function; CC, cellular component.
Genes 13 00366 g004
Genes 13 00366 g005 550
Figure 5. The top 20 enriched KEGG pathways among the DEGs. (A) Enriched KEGG pathways in comparison Lo 6 cm vs. Lo 4 cm. (B) Enriched KEGG pathways in Lo 8 cm vs. Lo 6 cm stages. The y-axis presents the different KEGG pathways and the x-axis presents the enrichment factor. Dot colors represent the q-value and dot sizes represent gene numbers enriched in the corresponding pathway.
Figure 5. The top 20 enriched KEGG pathways among the DEGs. (A) Enriched KEGG pathways in comparison Lo 6 cm vs. Lo 4 cm. (B) Enriched KEGG pathways in Lo 8 cm vs. Lo 6 cm stages. The y-axis presents the different KEGG pathways and the x-axis presents the enrichment factor. Dot colors represent the q-value and dot sizes represent gene numbers enriched in the corresponding pathway.
Genes 13 00366 g005
Genes 13 00366 g006 550
Figure 6. The overall expression of TFs in DEGs. (A) Numbers of total differentially expressed transcription factors. (B) Heatmap diagrams showing the relative expression levels of differentially expressed transcription factors. The color scale on the right from blue to red represents FPKM values from low to high.
Figure 6. The overall expression of TFs in DEGs. (A) Numbers of total differentially expressed transcription factors. (B) Heatmap diagrams showing the relative expression levels of differentially expressed transcription factors. The color scale on the right from blue to red represents FPKM values from low to high.
Genes 13 00366 g006
Genes 13 00366 g007 550
Figure 7. Heat map diagrams showing the relative expression levels of TFs possibly associated with the development of tapetum and microspore. (A) NAC (NAM/ATAF/CUC) family member. (B) MYB family members. (C) bZIP (basic region/leucine zipper) family members. (D) bHLH (basic helix-loop-helix) family members. (E) MADS-box family members. (F) PHD (plant homeodomain)-finger protein family members. The color scale on the top from blue to red represents FPKM values from low to high.
Figure 7. Heat map diagrams showing the relative expression levels of TFs possibly associated with the development of tapetum and microspore. (A) NAC (NAM/ATAF/CUC) family member. (B) MYB family members. (C) bZIP (basic region/leucine zipper) family members. (D) bHLH (basic helix-loop-helix) family members. (E) MADS-box family members. (F) PHD (plant homeodomain)-finger protein family members. The color scale on the top from blue to red represents FPKM values from low to high.
Genes 13 00366 g007
Genes 13 00366 g008 550
Figure 8. The heatmap expression representation of candidate DEGs involved in the tapetum degradation and the development of microspores. (A) Cell-cycle- and division-related genes. (B) Cell-wall-formation-related genes. (C) Transcription-factor-related genes. (D) Leucine-rich repeat receptor-like kinase (LRR-RLK)-related genes. (E) Plant-hormone-synthesis- and signal-transduction-related genes. The color scale on the top from blue to red represents FPKM values from low to high.
Figure 8. The heatmap expression representation of candidate DEGs involved in the tapetum degradation and the development of microspores. (A) Cell-cycle- and division-related genes. (B) Cell-wall-formation-related genes. (C) Transcription-factor-related genes. (D) Leucine-rich repeat receptor-like kinase (LRR-RLK)-related genes. (E) Plant-hormone-synthesis- and signal-transduction-related genes. The color scale on the top from blue to red represents FPKM values from low to high.
Genes 13 00366 g008
Genes 13 00366 g009 550
Figure 9. Relative expression analysis of the candidate DEGs. The values of qRT-PCR were determined using the 2−ΔΔCT method, the lowercase letters “a”, “b” and “c” represent statistically significant differences calculated by an ANOVA and post hoc Tukey’s HSD (p < 0.05). FPKM values, listed below the x-axis, are the counts from RNA-seq data.
Figure 9. Relative expression analysis of the candidate DEGs. The values of qRT-PCR were determined using the 2−ΔΔCT method, the lowercase letters “a”, “b” and “c” represent statistically significant differences calculated by an ANOVA and post hoc Tukey’s HSD (p < 0.05). FPKM values, listed below the x-axis, are the counts from RNA-seq data.
Genes 13 00366 g009
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Sui, J.; Jia, W.; Xin, Y.; Zhang, Y. Transcriptomics-Based Identification of Genes Related to Tapetum Degradation and Microspore Development in Lily. Genes 2022, 13, 366. https://doi.org/10.3390/genes13020366

AMA Style

Sui J, Jia W, Xin Y, Zhang Y. Transcriptomics-Based Identification of Genes Related to Tapetum Degradation and Microspore Development in Lily. Genes. 2022; 13(2):366. https://doi.org/10.3390/genes13020366

Chicago/Turabian Style

Sui, Juanjuan, Wenjie Jia, Yin Xin, and Yuanyuan Zhang. 2022. "Transcriptomics-Based Identification of Genes Related to Tapetum Degradation and Microspore Development in Lily" Genes 13, no. 2: 366. https://doi.org/10.3390/genes13020366

APA Style

Sui, J., Jia, W., Xin, Y., & Zhang, Y. (2022). Transcriptomics-Based Identification of Genes Related to Tapetum Degradation and Microspore Development in Lily. Genes, 13(2), 366. https://doi.org/10.3390/genes13020366

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop