Open access peer-reviewed chapter

Plant Comparative Transcriptomics Reveals Functional Mechanisms and Gene Regulatory Networks Involved in Anther Development and Male Sterility

Written By

Xiangyuan Wan and Ziwen Li

Submitted: March 23rd, 2019 Reviewed: July 1st, 2019 Published: July 30th, 2019

DOI: 10.5772/intechopen.88318

Chapter metrics overview

1,195 Chapter Downloads

View Full Metrics


Gene transcription and transcriptional regulation are crucial biological processes in all cellular life. Through the next-generation sequencing (NGS) technology, transcriptome data from different tissues and developmental stages can be easily obtained, which provides us a powerful tool to reveal the transcriptional landscape of investigated tissue(s) at special developmental stage(s). Anther development is an important process not only for sexual plant reproduction but also for genic male sterility (GMS) used in agriculture production. Plant comparative transcriptomics has been widely used to uncover molecular mechanism of GMS. Here, we focused on researches of anther developmental process and plant GMS genes by using comparative transcriptomics method. In detail, the contents include the following: (1) we described the commonly used flowchart in comparative transcriptomics; (2) we summarized the comparative strategies used to analyze transcriptome data; (3) we presented a case study on a maize GMS gene, ZmMs33; (4) we described the methods and results previously reported on gene co-expression and gene regulatory networks; (5) we presented the workflow of a case study on gene regulatory network reconstruction. The further development of comparative transcriptomics will provide us more powerful theoretical and application tools to investigate molecular mechanism underlying anther development and plant male sterility.


  • plant comparative transcriptomics
  • gene regulatory network
  • anther development
  • genic male sterility
  • molecular mechanism

1. Introduction

Gene transcription is an important biological process by which genetic information stored within DNA molecules is transmitted to RNA molecules according to the “genetic central dogma” in molecular biology [1]. After completion of the human genome project, the researchers began to reveal the transcriptional landscape of all genes in a genome to further investigate the functional mechanisms underlying phenotypic variations at a genome-wide transcriptional level. Therefore, biological studies on high-throughput omics data run from the genomic level into the transcriptomic level. Transcriptome data includes biological information of gene transcriptional activities in a certain cell, a tissue, or an individual (a population of cells) and even in a pool of samples under a certain developmental stage, an environmental condition, or an experimental treatment. Compared with other omics data (e.g., data of genome, epigenome, proteome, metabolome, or phenome), the primary characteristic of transcriptome data is that it includes temporal–spatial bioinformation affected by diverse developmental stages, tissue types, and internal/external environment events. Therefore, transcriptome data is more complex than genome data.

Transcriptomic studies usually focus on the transcriptional content and gene regulations in a genome. Gene expression microarray (GEM) is an early developed but still-utilized biotechnology by which the genome-wide transcription information can be obtained for genome-sequenced or transcriptional loci available species. In 1995, Schena et al. monitored expression levels of 48 genes by GEM in Arabidopsis thaliana [2], and then GEM was gradually and widely used for the estimation of gene expression levels. Until 2013, the amount of transcripts monitored by one microarray had been reached to more than 285,000 in human transcriptomics studies (the human transcriptome array). GEM is a hybrid-based method, while the sequencing-based method has been developed much faster and became one of the most commonly used biotechnologies in scientific studies and applications related to disease diagnosis [3]. Serial analysis of gene expression (SAGE) proposed by Velculescu et al. [4] and massively parallel signature sequencing (MPSS) reported by Brenner et al. [5] are two earlier developed sequencing-based methods to estimate the transcription information at a genome level. Nowadays, the majority of transcriptome data are generated by the NGS-based RNA sequencing (RNA-seq). RNA-seq technology combining with the following developed comparative transcriptomics analysis flowchart that is mainly based on digital gene expression profile (DGEP) is a commonly used research strategy in biological studies at molecular and genomic levels.

Anther is an important organ in sexual plant reproduction. Anther development is a dynamic process from the identity of the stamen to the production of mature pollen grains. During this period, two-thirds of protein-coding genes are transcribed, and more than 6% of them are anther specific (a reanalyzed result based on [6]). Thus, the anther transcriptome is specific and complex compared with transcriptomes of other plant organs. Plant comparative transcriptomics is an effective strategy used to investigate the molecular mechanism underlying anther developmental process. The comparative method based on anther transcriptomes can be performed between different genotypes, different developmental stages, different types of anther cells, and different biotic or abiotic treatments and even between different plant species. Consequently, differentially expressed genes (DEGs) are identified from above comparisons. Based on the comparison results, functionally important coding genes and noncoding transcripts including long noncoding RNAs (lncRNAs), microRNAs (miRNAs), and other small RNAs could be uncovered. However, the goal of plant comparative transcriptomics is not only to identify DEGs but also to reconstruct gene regulatory relationships of the upstream regulators and the downstream regulated targets of the investigated genes. In this review, based on anther transcriptomes, we first summarized the research workflow commonly used in the experimental design and data analyses in plant transcriptomics studies, and then we described several types of comparison strategies in comparative transcriptomics using anther transcriptome data as the analyzed example. In the following section, we generally discussed gene regulatory and co-expression networks used to investigate the molecular foundation of developing anther in a network-based perspective. Additionally, we described two case studies in our laboratory to explain the detailed analysis processes and applications of comparative transcriptomics in plant GMS gene studies.


2. Comparative analysis using transcriptome data

In comparative transcriptomics, the commonly used pipeline to identify potential functional genes and to reveal the gene functions, as well as to investigate the regulatory relationships between these genes, includes five aspects. They are data preparation, DGEP analysis, DEG analysis, gene set enrichment (GSE) analysis, and gene regulatory network (GRN) analysis, respectively (Figure 1). These five aspects are closely connected in the whole pipeline, and the corresponding analyses mainly depend on data management skills in bioinformatics.

Figure 1.

A flowchart of comparative transcriptomics analysis.

The basic application of comparative transcriptomics is to obtain a transcriptional landscape of the investigated biological sample. It is composed of not only the estimated transcription levels of annotated transcribed loci along the genomes (the known genomic loci with reported or predicted transcription abilities) but also the identification of novel transcribed loci (the stably transcribed loci not annotated or identified in previous studies). More importantly, in current biological studies, transcribed loci identified by researchers include not only the protein-coding genes but also lncRNAs and other noncoding RNAs. Both GEM and RNA-seq technologies can be used to uncover the genome-wide profiles of transcription levels of annotated genes. However, the identification of novel transcribed loci can be only effectively performed by RNA-seq method and the following DGEP analysis. This is one reason why RNA-seq is more commonly used in transcriptomics studies. Moreover, GEM method depends on hybridization probes that are designed based on known whole genome sequence or an appreciable set of sequenced transcripts (e.g., expressed sequence tags) of the investigated species, which restrict its application on some species without whole genome information or sequence resource. On the contrary, the sequencing-based method of RNA-seq can be applicable for species without sequenced genomes. This is another reason for the popularity of RNA-seq. In genome available species, RNA-seq data should be firstly mapped to the reference genome (Figure 1).

A gene with its transcription levels significantly different between two groups of samples is defined as a DEG under a certain comparison condition (Figure 1). It is notable that the concept of DEG specially represents the expression changes of protein-coding genes at the earlier stages of expression data analysis. However, along with the rapid development of molecular biology and the deeper understanding on the functional element on the genome, the concept of DEG has been expanded to noncoding transcripts, for example, the differentially expressed (DE) miRNA and the DE lncRNA. Furthermore, if both coding and noncoding transcripts are considered in the comparative analysis of transcriptome data, transcriptional alterations between control and treated samples should be defined as DE transcribed loci or DE loci. Thus, DE loci is a broad concept used to describe transcriptional alterations of genetic element. There are several strategies for comparing transcriptomes from different research subjects to identify DE loci (described in Section 3, “Plant comparative transcriptomics in anther”).

Identified DEG set or DE loci should be appropriately annotated with functional descriptions to determine which biological process or pathway these DEGs are involved in. In comparative transcriptomics, this step is a critical bridge linking transcriptional changes to gene functions and even gene regulation networks. Two commonly utilized methods to annotate DEGs consist of the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene enrichment analyses. Both of them belong to GSE analysis (Figure 1). The GO database includes tens of thousands of GO terms, and each GO term contains several genes with the same biological function. Each gene has three functional aspects, including molecular function (the molecular activities of gene products), cellular component (the cellular locations of functional gene products), and biological process (the gene products’ molecular functions with biological process). GSE analysis based on GO database provides some basic functional descriptions for the investigated DEG set. KEGG analysis is a pathway-based enrichment method. The KEGG database has accumulated hundreds of metabolism pathways in plants, animals, and other species. Thus, KEGG analysis can reveal significant pathways the DEGs participated in. GO-based methods can annotate more genes than KEGG-based method, as the GO terms are more flexible and include a larger number of genes. On the other hand, because most metabolic pathways are conserved across species and more significant in biological processes, annotated results obtained from KEGG-based method may be more conserved and stable. In comparative transcriptomics, GO- and KEGG-based analyses are together utilized in gene function studies.

The locations of transcribed loci on the genome, their transcription levels, and the changed expression can be identified through comparative transcriptomics analysis. The detected DEG set represents a functional gene set related to the function of investigated gene, the phenotype variation, the stress resistance ability, or the development process. Furthermore, gene regulation relationships are the underlining molecular mechanism of altered transcriptomes, and novel gene regulatory networks could be uncovered by comparative transcriptomics analysis (Figure 1). Several types of gene regulatory relationships and the reconstructions of gene regulatory networks based on plant comparative transcriptomics are described and discussed in Section 5 (“Gene co-expression and regulatory networks reconstructed by comparative transcriptomics method”).


3. Plant comparative transcriptomics in anther

One of the major subjects of modern molecular biology is to uncover the functions of genes in the genome and reveal the molecular mechanism of phenotypic variation. Gene transcription levels and their changes in different conditions are important information that can reflect the functions and transcriptional regulation relationships of investigated genes. How to estimate the transcription levels of genes and how to obtain the transcriptional landscape of a genome are two major subjects in biological studies on gene expression. DGEP and DEG analyses are powerful tools to solve these questions. In DEG analysis, according to the scientific or application questions, the comparison strategies between investigated biological samples are classified into six types including (1) different genotypes, (2) different developmental stages, (3) different tissues, (4) different cell types, (5) different treatments, and (6) different species (Figure 2). Here, as we mainly focus on comparative transcriptomics analysis on the developmental anther tissues and the interspecies analysis on anther transcriptome data being rare, the third and sixth types will not be discussed.

Figure 2.

Comparative transcriptomics strategies.

3.1 Different genotypes

There are two types of genotype-based transcriptome data between wild type (WT) and mutant lines in GMS studies, which are based on whether the causal mutation is known or not (Table 1). For transcriptomes of male sterility (MS) lines with known causal mutations, the comparison of transcriptomes between WT and MS lines will identify many DEGs associated with the function loss or expression change of the investigated mutation locus. If the causal mutation has not been identified from the MS line, comparative transcriptomics analyses will provide the researchers important results related to the unsettled genetic difference, such as how many genes are changed in expression levels in the MS lines and what the functions of these genes are, even though the causal mutation candidates can be inferred from these genes if the researchers have primary mapping results.

Plants MS gene or MS line Method Tissue Data sourcea Reference
Arabidopsis thaliana ROXY1 and ROXY2 Microarray Young inflorescences SD [7]
AMS Microarray Anthers GSE18225 [8]
AMS Microarray Anthers SD [9]
AMS Microarray Floral buds SD [10]
EMS1 Microarray Anthers SD [9]
EMS1 Microarray Anthers SD [11]
Ms1 Microarray Young closed buds SD [12]
Ms1 Microarray Floral buds GSE8864 [13]
ICE1 RNA-seq Anthers GSE107260 [14]
DYT1 Microarray Anthers GSE18225 [8]
CDM1 Microarray Young floral buds GSE55799 [15]
TEK Microarray Closed floral buds GSE56497 [16]
bHLH010, bHLH08, bHLH091 RNA-seq Anthers SRS838170, SRS838173 [17]
Oryza sativa PTC1 Microarray Anthers SD [18]
UDT1 Microarray Anthers GSE2619 [19]
OsGAMYB Microarray Anthers SD [20]
TDR Microarray Spikelets SD [21]
MADS3 Microarray Anthers SD [22]
Zea mays Ms23 RNA-seq Anthers GSE90849 [23]
Ms32 Microarray Anthers GSE90968 [23]
MAC1 Microarray Anthers SD [24]
Triticum aestivum TaMs1 RNA-seq Anthers SRP113349 [25]
TaMs2 RNA-seq Anthers SRP092366 [26]
Solanum lycopersicum ms1035 RNA-seq Floral buds SD [27]
MS line 7B-1 RNA-seq Anthers GSE85859 [28]
Brassica napus MS line WSLA RNA-seq Young flower buds SRR2192464, SRR2192489 [29]
MS line SP2S RNA-seq Young flower buds GSE69638 [30]
MS line TE5A RNA-seq Young flower buds SRP068170 [31]
Citrullus lanatus MS line DAH3615-MS RNA-seq Floral buds and flowers GSE69073 [32]

Table 1.

Published studies on anther transcriptome data between WT and MS lines.

“SD” indicates the raw data is unavailable, while the up- and downregulated genes are listed in the supplemental data (SD) in references cited.

3.2 Different developmental stages

The phenotypic differences among tissues and organs (e.g., root, leaf, and flower in plant) due to their differences of transcriptome landscape are well known. Furthermore, it is a developmental process for most types of plant organs from the organ identity (e.g., meristematic cells) to the final mature organ. Thus, how to reveal the dynamic changes of gene transcription levels and how to explain the morphological alterations regulated by gene expression changes are important tasks in plant comparative transcriptomics studies.

Meiosis is an important step in gametophyte generation process and sexual plant reproduction. Morphologic changes during cell meiosis process have been well described by cellular level investigations, while the molecular level alterations and their corresponding gene regulatory networks are not well understood. Plant transcriptomes are a powerful dataset to estimate the gene expression changes and infer the regulatory roles of key genes. Based on GEM technology, Ma et al. investigated maize anther transcriptomes during seven developmental stages and found that transcriptomes during meiosis stages exhibited the lowest complexity [33]. Hollender et al. surveyed the gene transcription profiles of anther of woodland strawberry (Fragaria vesca) from developmental stages 7–12 and identified numerous F-Box genes induced in transcription levels at meiosis stage [34]. Besides, tapetum is the inner cell layer of anther with important functions in anther development and gametocyte maturation. The generation, development, and degradation of tapetum are fine regulated during the anther development, while the regulatory framework and the details are far from complete. Yue et al. identified 243 DEG and 108 stage-specific genes during four anther developmental stages in Hamelia patens [35]. Chen et al. investigated the expression of genes involving in tapetum development of male floral bud during eight developmental stages in Populus tomentosa [36]. Thus, anther transcriptome data during different developmental stages provide valuable data sources for anther development studies. By the combination of comparative transcriptomics and bioinformatics analyses, more key functional genes and the underlying regulatory mechanisms for anther development will be further revealed.

3.3 Different types of anther cells

The cytological structure of anther consists of four cell layers, including the epidermis, endothecium, middle layer, and tapetum, and the archesporial cells are directly surrounded by the tapetum. Thus, the transcriptome data of a whole anther tissue is a mixed gene expression data from diverse cell types with different functions in the anther development process. It is necessary to obtain transcriptional dynamics from different cell layers separately to investigate anther development and the underlying molecular mechanism at a cell type-specific level. Several studies have identified cell layer-specifically expressed genes (e.g., tapetum cells or microgametes). Ma et al. identified 104 MS-related and non-pollen expressed genes most specifically expressed in tapetum by comparative transcriptomics analysis on four diverse MS lines in Brassica oleracea [37]. The other way to obtain cell layer-specific transcriptome in anther is firstly separating the investigated cell layer by laser capture microdissection (LCM) technology and then performing RNA-seq or GEM experiment on the separated samples. This strategy has been successfully used in rice, maize, and woodland strawberry to identify the tapetum- or microgamete-specifically expressed genes and their expression dynamics [34, 38, 39]. A recent published research has investigated maize male meiosis using single-cell RNA sequencing (scRNA-seq) technology on pre-meiotic and meiotic cells from maize anthers, which greatly promoted studies on plant anther scRNA-seq [40]. The comparative studies on transcriptomic dynamics between different types of cells facilitate the deeper understanding of functions of specific cell layers on anther development.

3.4 Different treatments

At the reproductive stage, plant is more sensitive to external environment conditions. The abiotic stresses, such as high temperature, drought, and cold and freezing stresses, will critically affect the developmental process of anther and pollen in flowing plants. Though there have been numerous studies on stress resistance and response in plant, the regulatory pathways of stress response and their cross talk at molecular level should be further investigated for anther development. Additionally, more effective stress-resistant genes should be identified for the purpose of crop improvement. Plant comparative transcriptomics between normal and stress-treated plants provide a wide insight into the stress response mechanisms of plant during sexual reproductive stage. Zhang et al. investigated the genome-wide transcriptional changes of rice panicle under heat treatment (40°C) and found thousands of DEGs participating in transcriptional regulation, transport, cellular homeostasis, and stress response [41]. Studies on photosensitive or thermosensitive GMS lines can also reveal a lot of genes responding to environmental changes.


4. A case study: revealing the molecular functions of a MS gene, ZmMs33, by comparative transcriptomics

The discoveries of genes that play key roles in the development of maize anther provide important genetic resources for the utilization of heterosis in maize. Analysis of functional mechanism of GMS genes can effectively promote researches on anther development biology and deepen our understanding of molecular mechanism controlling sexual plant reproduction [42]. There are several published case studies containing comparative transcriptomics analysis on maize GMS genes in our laboratory, including ZmMs7 [43], ZmMs20 [44], ZmMs30 [45], and ZmMs33 [4647]. We used comparative transcriptomics analysis based on developmental anthers of ZmMs33 wild type and ms33–6038 mutant to analyze the transcription changes corresponding to male sterility phenotype and to further investigate the underlying molecular mechanisms of GMS regulated by ZmMs33 gene.

This ms33–6038 mutant is complete male sterility and displays small and pale-yellow anthers (Figure 3A). Transmission electron microscope (TEM) observation and dynamic scanning electron microscopy (SEM) analysis were performed to analyze the phenotypic alteration of anther wall layers, microspores, Ubisch bodies, and exine between wild type and ms33–6038 mutant during anther developmental stages (Figure 3A–C).

Figure 3.

Reveal ZmMs33 gene functions for anther development by comparative transcriptomics analysis. (A) Phenotype of whole plants (A1), anthers (A2), pollen grains (A3), and outer surface of anther wall (A4) of WT and ms33–6038 mutant. (B) TEM analysis of anther wall layers, microspores, Ubisch bodies, and exine in WT and ms33–6038 mutant. (C) SEM analysis of microspores and pollen grains in WT and ms33–6038 mutant. (D) Map-based cloning of ZmMs33 gene. (E) Phenotypes of tassels, anthers, and pollen grains in three ms33 knockout lines generated by a CRISPR/Cas9 system. (F) PCA analysis of RNA-seq data from WT and ms33–6038 mutant. (G) Venn plot of DEGs at each developmental stage. Figure 3AC was cited from [46]. Figure 3D and E was cited from [47].

Maize Zm00001d007714 was identified as ZmMs33 via a map-based cloning approach (Figure 3D). ZmMs33 encodes an esterase that belongs to gene family of glycerol-3-phosphate acyltransferase (GPAT) in maize. To further confirm gene function of Zm00001d007714, a CRISPR/Cas9 system was used to generate targeted knockout lines. Three types of T0-generation maize plants homozygous for null alleles of Zm00001d007714 were observed to be complete male sterility (Figure 3E), suggesting that function loss of Zm00001d007714 is the causal mutation for male sterile phenotype of the ms33 mutant.

Subsequently, RNA-seq was performed using anther tissues during developmental stages 5–9 to obtain a comprehensive transcriptional profile of WT and ms33-6038. Three biological samples were collected at each developmental stage for sequencing. After data preparation and transcription level estimation, we compared similarities of transcriptional profiles of protein-coding genes by principal component analysis (PCA) (Figure 3F) and found good repeatability among three biological repeats.

Finally, we identified DEGs between WT and mutant and between adjacent developmental stages, separately. We found that the amount of DEGs between WT and mutant at stages 5–7 was significantly smaller in magnitude than that at stages 8a–9 (Figure 3G), indicating that ms33 mutant transcriptomes are significantly divergent from WT starting from stage 8a. The transcriptome landscapes of WT were similar to those of ms33 mutant at stages 5–7. Besides, DEG amounts were various between adjacent developmental stages. It is worth noting that the DEG amount between WT and mutant exceeded that between adjacent stages from stage 8a–9. This result implied that the transcriptomes were significantly changed at the later three stages. Therefore, we compared the transcriptomes between genotypes at the former three and the later three stages, separately. In contrast to a limited number of DEGs (only two genes) shared by the former three stages, there were thousands of shared DEGs at the later three stages. GSE analysis based on KEGG database suggested that the upregulated gene set was firstly enriched in the function of biosynthesis of secondary metabolites, while the downregulated genes were significantly related to the photosynthesis process. This pathway enrichment analysis partly represents the alterations in metabolisms and physiological activities closely associated with the transcriptional changes caused by function defect of ms33.


5. Gene co-expression and regulatory networks reconstructed by comparative transcriptomics method

Though DEGs are mainly identified by pairwise comparisons between transcriptomes of tissues, stages, or treatment conditions and can reflect most of the transcriptional changes between two sets of samples, these transcriptional alterations are not sufficient to explain the detailed molecular mechanism underlying tissue-specific development processes and stress-resistant pathways. Moreover, the molecular functions of genes act under GRNs. All the biological processes of growth, development, stress response, and reproduction are regulated by GRNs. The prediction of gene regulatory relationships and the reconstruction of the GRNs by using the transcriptome data are also the major aims in transcriptomics studies, except for the DGEP and DEG analyses.

5.1 Gene co-expression analysis

Function-related genes tend to co-express in a cell, either to form a complex or to involve in the same biological pathway. Thus, the similar pattern of gene expressions can be used as an indicator to predict gene functions. Gene co-expression (GCE) analysis is a powerful tool to discover important functional genes in biological processes including anther development. A relatively early study identified two functional GMS genes, POLYKETIDE SYNTHASE A (PKSA) and PKSB, through detecting co-expressed genes with ACOS5, a GMS gene belonging to fatty acyl-CoA synthetase gene family, based on microarray data in A. thaliana [48]. Similarly, ABORTED MICROSPORES (AMS) gene was reported participating in the pollen wall formation in rice by the analyses of 98 co-expressed genes with AMS in flower development [49]. GCE analysis can be also used to investigate the biological functions and the regulatory targets of a gene. This genome-wide analysis on GCE networks has been performed based on microarray data from A. thaliana anther tissues, and 254 complete GCE groups containing 10,513 anther-transcribed genes were revealed [50]. Another microarray-based GCE network was reconstructed in A. thaliana anther by using 10,797 genes expressed in anther/flora, and transcriptional landscape of GMS mutant was included in the stable examination of this newly constructed network [51]. In rice, microarrays from WT anther tissue across stages 2–14 and nine GMS lines were integrated to reconstruct a big GCE network containing more than 9000 genes and 0.4 million pairs of co-expression relationships [52].

RNA-seq data-based GCE network analysis was performed in anther when high-throughput sequencing technology was developed. In woodland strawberry, stages 1–12 floral samples dissected by LCM or hand, including stages 6–12 anther tissues, were sequenced by RNA-seq. Gene co-expression network analysis was used to reconstruct GCE networks in strawberry’s flower development, and 23 modules were discovered from the GCE networks including 4584 pollen-specific genes [34]. These genome-wide GCE networks are useful for characterization of genes associated with anther development and floral reproduction.

5.2 TF-encoding gene regulatory network

Genes with their products forming one protein complex, genes encoding transcription factor (TF) and TF target genes, and genes functioning in the same metabolic pathway or stress-resistant process often tend to be co-expressed in a cell. Therefore, the expression-associated genes in GCE network may be not directly functionally linked. A more accurate and robust gene regulatory network is needed for both the biological function and network researches at molecular and genome levels. One way to improve the gene regulatory network is to introduce gene regulatory types into the network. Several TF gene regulatory networks (TF-GRN), also called as transcriptional regulatory network (TRN), were reconstructed based on expression patterns of TF-encoding genes and TF target genes from transcriptome data. One TF-GRN comprised 19 TFs and their 101 target genes involving in A. thaliana pollen development [53]. Another GRN of early anther development was constructed by interactively analyzing transcriptome data from three GMS lines of TF-encoding gene knockout mutants [9]. In the maize genome, there are 2298 TF-encoding genes identified which belonged to 56 diverse families [54]. A total of 3078 TF-encoding genes belonging to 59 families are predicted in silico analysis in rice genome [55]. These TF databases, combining with increased amount of transcriptome data from mutants of TF-encoding genes and other omics data (e.g., Chip-seq, DAP-seq), provide abundant data for the reconstruction of TF-GRN with increased credibility, applicability, and completeness.

5.3 miRNA target gene regulatory network

Both transcriptional and posttranscriptional regulations are crucial in controlling the normal development and stress-resistant process in cellular life. The miRNA-mediated regulation model on target genes is a well-studied posttranscriptional gene regulation pathway that plays important roles in floral identification and the following development of flower organs [56, 57, 58] as well as male fertility [5960]. Beyond numerous case studies on functional miRNAs in anther development and GMS genes [61, 62, 63, 64], the expression profile of miRNAs and the regulatory networks were investigated to elevate our understanding on the transcriptional regulatory mechanism of miRNAs. GRNs between miRNA and their target genes have been constructed via flower/anther transcriptomics in the model plant species, A. thaliana, and some other plants [65, 66, 67, 68]. Furthermore, comparative transcriptomics analysis on small miRNAs has been commonly used as a research method to reveal the transcriptional alterations between fertility and sterility lines in economic and food plant species, such as maize [45], tomato [69], cotton [70, 71], wheat [72, 73], pine [74], lycium [75], watermelon [32], and Brassica campestris [76].

5.4 ceRNA-miRNA regulatory network

It is well known that miRNAs are crucial regulators on gene expressions that control key biological functions including anther development, since miRNA was firstly found in nematodes in 1993 [77]. It is noteworthy that a novel type of gene regulatory model, the competing endogenous RNA (ceRNA) hypothesis, was recently proposed [78]. According to the ceRNA hypothesis, some endogenous transcripts have abilities to adsorb miRNA molecules; subsequently, the expression levels of miRNA target genes can be derepressed [78, 79]. A typical ceRNA in plant, a long noncoding RNA, IPSI, was found in A. thaliana. It could completely sponge miRNA ath-miR399 and indirectly increase the transcription levels of an important gene involved in phosphate homeostasis [80]. The following studies revealed that transcripts of protein-coding genes, pseudogene, transposable elements, simple sequence repeat, and circular RNAs have molecular functions as ceRNAs [79, 81, 82], indicating that the ceRNA-miRNA relationship is an essential gene regulatory mechanism in the growth and development of plants and animals. Consequently, it is necessary to introduce ceRNA regulators into GRN construction. Here, we present our recent study on reconstructing ceRNA regulatory network mainly based on RNA-seq and small RNA-seq transcriptomes from developmental maize anther.


6. A case study: reconstructing ceRNA-miRNA target gene regulatory networks using transcriptome data of maize anther

Here we summarized the research progress of one recently completed research related to the ceRNA-mediated GRN in our laboratory. Generally speaking, this is the first study introducing ceRNA regulation into miRNA target gene regulatory pathway for deeply dissecting the mechanism of anther development and sexual plant reproduction at a network level. This provides a fresh example for GRN research by plant comparative transcriptomics and has dual significance in both theoretical and practical senses. It may also provide new thoughts and strategies for further transcriptome-based GRN studies.

It is well known that gene expressions are controlled by the GRN in cellular life. Newly found regulatory patterns (e.g., miRNA pathway and epigenetic modification) have enhanced our understanding on the GRN. Recently, “ceRNA hypothesis” was proposed as a novel type of gene regulatory relationship and was found to participate in different development and stress response processes of organisms by a number of case studies. However, the network level study on ceRNA regulatory functions is still rare, which limited our deep understanding on the GRN. In addition, studies on the GRN of sexual plant reproduction and male sterility are crucial for both fundamental biological significance and applications in plant hybrid breeding and seed production. We investigated ceRNA-miRNA target gene regulatory network in maize anther developmental process by plant comparative transcriptomics method. Six steps were performed from raw sequencing data preparation to the finally constructed GRN (Figure 4). Firstly, we performed RNA- and small RNA-seq using anther tissues at three developmental stages from two maize lines to obtain a relative broad transcriptional landscape in anther development and transcribed loci that are stably expressed in maize species. Secondly, we identified stably transcribed loci based on the maize reference genome and estimated their transcription levels. In this step, we only used shared transcription loci identified from RNA-seq data between two maize lines (Figure 4A). Notably, these transcribed loci were divided into five groups such as protein-coding genes, lncRNAs, transposable elements, and unassigned loci. Thirdly, we identified known miRNAs and predicted potential novel miRNAs that may be involved in maize anther development. Sequenced small RNA data were obtained from the same samples that were used in RNA-seq. A matched dataset (e.g., matched RNA and small RNA sequenced dataset here) is important in experimental design and more powerful to reveal the investigated biological questions. Though the analysis workflow of small RNA-seq data is similar to that of RNA-seq data in general (Figure 1), there are some differences between them. In our analysis, we reanalyzed two sets of published small RNA data to compare with their results from our own sequenced data for credible known and potential novel miRNAs involved in maize anther development [23, 83] (Figure 4B). This is an important check method to confirm the stability of research results and conclusions. Fourthly, we predicted ceRNA-miRNA interaction pairs and miRNA target gene regulatory pairs by computational approach (Figure 4C). Bioinformatics analysis in this step is mainly based on genome sequence but not the transcriptomes. Fifthly, we reconstructed ceRNA-miRNA target gene regulatory networks by predicted interaction pairs and transcription correlation patterns from transcriptomics data (Figure 4D). It is well known that miRNAs could repress the transcription levels of their target genes. Additionally, ceRNA was demonstrated to negatively regulate the transcription levels of matched miRNAs. The negatively associated gene pairs in transcription levels may be more credible in mutual interactions. By integrating ceRNA-miRNA and miRNA target gene interactions, we reconstructed ceRNA-miRNA target gene regulatory networks in maize anther. Finally, we generally investigated the functional significance of genes in the regulatory network by GO enrichment analysis. In these networks, we found a number of well-studied genes and miRNA target gene pairs involved in maize anther development and male sterility, suggesting that the ceRNA-miRNA target gene regulatory networks contribute to anther development in maize. Besides, GO analysis of target genes in the network revealed that they are functionally enriched in flower development process (Figure 4E) [84].

Figure 4.

A flowchart of reconstructing the ceRNA-miRNA target gene regulatory network in developmental maize anther. (A) Identification and classification of stably transcribed loci in maize anther. (B) Identification of known miRNA in maize anther. (C) Prediction of ceRNA-miRNA and miRNA target gene interaction pairs. (D) Reconstruction of ceRNA-miRNA target gene regulatory networks. (E) GSE analysis of target genes in the networks.


7. Conclusions

Here, we summarized major points in comparative transcriptomics analysis from the commonly utilized workflow to the closely related research cases and from the single gene-based function analysis to GRN-based gene function investigation. In GMS gene studies, the research experiments using comparative transcriptomics method to investigate key functional genes and the genome-wide GRNs in developmental anther will facilitate our systematical understanding on the biological processes and molecular regulatory networks for anther development and sexual plant reproduction. More importantly, case studies illustrated here have a general meaning on technologies and methodologies for functional researches of other biological pathways and processes. With the fast advancement of sequencing technology, plant comparative transcriptomics has achieved considerable development. However, our understanding on the transcriptional dynamics and gene regulatory relationships of biological processes are far from being completed. Consequently, more efforts are needed for the further improvement of comparative transcriptomics in plant biological studies.



The research in our lab was supported by the National Key Research and Development Program of China (2018YFD0100806, 2017YFD0102001, 2017YFD0101201), the National Transgenic Major Program of China (2018ZX0800922B, 2018ZX0801006B), the National Natural Science Foundation of China (31,771,875, 31,871,702), the Fundamental Research Funds for the Central Universities of China (06500060), and the “Ten Thousand Plan” of National High-level Talents Special Support Plan (For Xiangyuan Wan).


Conflict of interest

The authors declare that they have no conflict of interest. Advertisement



competing endogenous RNA


differentially expressed


digital gene expression profile


gene co-expression


differentially expressed gene


gene expression microarray


genic male sterility


gene ontology


gene regulatory network


gene set enrichment


Kyoto encyclopedia of genes and genomes


laser capture microdissection


long noncoding RNA




massively parallel signature sequencing


male sterility


the next-generation sequencing


RNA sequencing


serial analysis of gene expression


single-cell RNA sequencing


scanning electron microscopy


transmission electron microscope


transcription factor


TF gene regulatory network


transcriptional regulatory network


wild type


  1. 1. Crick F. Central dogma of molecular biology. Nature. 1970;227:561-563
  2. 2. Schena M, Shalon D, Davis RW, Brown PO. Quantitative monitoring of gene expression patterns with a complementary DNA microarray. Science. 1995;270:467-470
  3. 3. Byron SA, Van Keuren-Jensen KR, Engelthaler DM, Carpten JD, Craig DW. Translating RNA sequencing into clinical diagnostics: Opportunities and challenges. Nature Reviews. Genetics. 2016;17:257-271
  4. 4. Velculescu VE, Zhang L, Vogelstein B, Kinzler KW. Serial analysis of gene expression. Science. 1995;270:484-487
  5. 5. Brenner S, Johnson M, Bridgham J, Golda G, Lloyd DH, Johnson D, et al. Gene expression analysis by massively parallel signature sequencing (MPSS) on microbead arrays. Nature Biotechnology. 2000;18:630-634
  6. 6. Stelpflug SC, Sekhon RS, Vaillancourt B, Hirsch CN, Buell CR, de Leon N, et al. An expanded maize gene expression atlas based on RNA sequencing and its use to explore root development. Plant Genome. 2016;9:1-15
  7. 7. Xing S, Zachgo S. ROXY1 and ROXY2, two Arabidopsis glutaredoxin genes, are required for anther development. The Plant Journal. 2008;53:790-801
  8. 8. Feng B, Lu D, Ma X, Peng Y, Sun Y, Ning G, et al. Regulation of the Arabidopsis anther transcriptome by DYT1 for pollen development. The Plant Journal. 2012;72:612-624
  9. 9. Ma X, Feng B, Ma H. AMS-dependent and independent regulation of anther transcriptome and comparison with those affected by other Arabidopsis anther genes. BMC Plant Biology. 2012;12:23
  10. 10. Xu J, Yang C, Yuan Z, Zhang D, Gondwe MY, Ding Z, et al. The aborted microspores regulatory network is required for postmeiotic male reproductive development in Arabidopsis thaliana. Plant Cell. 2010;22:91-107
  11. 11. Wijeratne AJ, Zhang W, Sun Y, Liu W, Albert R, Zheng Z, et al. Differential gene expression in Arabidopsis wild-type and mutant anthers: Insights into anther cell differentiation and regulatory networks. The Plant Journal. 2007;52:14-29
  12. 12. Yang C, Vizcay-Barrena G, Conner K, Wilson ZA. Male sterility1 is required for tapetal development and pollen wall biosynthesis. Plant Cell. 2007;19:3530-3548
  13. 13. Alves-Ferreira M, Wellmer F, Banhara A, Kumar V, Riechmann JL, et al. Global expression profiling applied to the analysis of Arabidopsis stamen development. Plant Physiology. 2007;145:747-762
  14. 14. Wei D, Liu M, Chen H, Zheng Y, Liu Y, Wang X, et al. Inducer of CBF expression 1 is a male fertility regulator impacting anther dehydration in Arabidopsis. PLoS Genetics. 2018;14:e1007695
  15. 15. Lu P, Chai M, Yang J, Ning G, Wang G, Ma H. The Arabidopsis callose defective microspore1 gene is required for male fertility through regulating callose metabolism during microsporogenesis. Plant Physiology. 2014;164:1893-1904
  16. 16. Lou Y, Xu XF, Zhu J, Gu JN, Blackmore S, Yang ZN. The tapetal AHL family protein TEK determines nexine formation in the pollen wall. Nature Communications. 2014;5:3855
  17. 17. Zhu E, You C, Wang S, Cui J, Niu B, Wang Y, et al. The DYT1-interacting proteins bHLH010, bHLH089 and bHLH091 are redundantly required for Arabidopsis anther development and transcriptome. The Plant Journal. 2015;83:976-990
  18. 18. Li H, Yuan Z, Vizcay-Barrena G, Yang C, Liang W, Zong J, et al. Persistent tapetal cell1 encodes a PHD-finger protein that is required for tapetal cell death and pollen development in rice. Plant Physiology. 2011;156:615-630
  19. 19. Jung KH, Han MJ, Lee YS, Kim YW, Hwang I, Kim MJ, et al. Rice undeveloped tapetum1 is a major regulator of early tapetum development. Plant Cell. 2005;17:2705-2722
  20. 20. Aya K, Ueguchi-Tanaka M, Kondo M, Hamada K, Yano K, Nishimura M, et al. Gibberellin modulates anther development in rice via the transcriptional regulation of GAMYB. Plant Cell. 2009;21:1453-1472
  21. 21. Zhang DS, Liang WQ , Yuan Z, Li N, Shi J, Wang J, et al. Tapetum degeneration retardation is critical for aliphatic metabolism and gene regulation during rice pollen development. Molecular Plant. 2008;1:599-610
  22. 22. Hu L, Liang W, Yin C, Cui X, Zong J, Wang X, et al. Rice MADS3 regulates ROS homeostasis during late anther development. Plant Cell. 2011;23:515-533
  23. 23. Nan GL, Zhai J, Arikit S, Morrow D, Fernandes J, Mai L, et al. MS23, a master basic helix-loop-helix factor, regulates the specification and development of the tapetum in maize. Development. 2017;144:163-172
  24. 24. Zhang H, Egger RL, Kelliher T, Morrow D, Fernandes J, Nan GL, et al. Transcriptomes and proteomes define gene expression progression in pre-meiotic maize anthers. G3 (Bethesda). 2014;4:993-1010
  25. 25. Wang Z, Li J, Chen S, Heng Y, Chen Z, Yang J, et al. Poaceae-specific MS1 encodes a phospholipid-binding protein for male fertility in bread wheat. Proceedings of the National Academy of Sciences of the United States of America. 2017;114:12614-12,619
  26. 26. Ni F, Qi J, Hao Q , Lyu B, Luo MC, Wang Y, et al. Wheat Ms2 encodes for an orphan protein that confers male sterility in grass species. Nature Communications. 2017;8:15121
  27. 27. Jeong HJ, Kang JH, Zhao M, Kwon JK, Choi HS, Bae JH, et al. Tomato male sterile 1035 is essential for pollen development and meiosis in anthers. Journal of Experimental Botany. 2014;65:6693-6709
  28. 28. Omidvar V, Mohorianu I, Dalmay T, Zheng Y, Fei Z, Pucci A, et al. Transcriptional regulation of male-sterility in 7B-1 male-sterile tomato mutant. PLoS One. 2017;12:e0170715
  29. 29. Qu C, Fu F, Liu M, Zhao H, Liu C, Li J, et al. Comparative transcriptome analysis of recessive male sterility (RGMS) in sterile and fertile Brassica napus lines. PLoS One. 2015;10:e0144118
  30. 30. Liu XQ , Liu ZQ , Yu CY, Dong JG, Hu SW, Xu AX. TGMS in rapeseed (Brassica napus) resulted in aberrant transcriptional regulation, asynchronous microsporocyte meiosis, defective tapetum, and fused sexine. Frontiers in Plant Science. 2017;8:1268
  31. 31. Yan X, Zeng X, Wang S, Li K, Yuan R, Gao H, et al. Aberrant meiotic prophase I leads to genic male sterility in the novel TE5A mutant of Brassica napus. Scientific Reports. 2016;6:33955
  32. 32. Rhee SJ, Seo M, Jang YJ, Cho S, Lee GP. Transcriptome profiling of differentially expressed genes in floral buds and flowers of male sterile and fertile lines in watermelon. BMC Genomics. 2015;16:914
  33. 33. Ma J, Skibbe DS, Fernandes J, Walbot V. Male reproductive development: Gene expression profiling of maize anther and pollen ontogeny. Genome Biology. 2008;9:R181
  34. 34. Hollender CA, Kang C, Darwish O, Geretz A, Matthews BF, Slovin J, et al. Floral transcriptomes in woodland strawberry uncover developing receptacle and anther gene networks. Plant Physiology. 2014;165:1062-1075
  35. 35. Yue L, Twell D, Kuang Y, Liao J, Zhou X. Transcriptome analysis of Hamelia patens (Rubiaceae) anthers reveals candidate genes for tapetum and pollen wall development. Frontiers in Plant Science. 2016;7:1991
  36. 36. Chen Z, Rao P, Yang X, Su X, Zhao T, Gao K, et al. A global view of transcriptome dynamics during male floral bud development in Populus tomentosa. Scientific Reports. 2018;8:722
  37. 37. Ma Y, Kang J, Wu J, Zhu Y, Wang X. Identification of tapetum-specific genes by comparing global gene expression of four different male sterile lines in Brassica oleracea. Plant Molecular Biology. 2015;87:541-554
  38. 38. Hirano K, Aya K, Hobo T, Sakakibara H, Kojima M, Shim RA, et al. Comprehensive transcriptome analysis of phytohormone biosynthesis and signaling genes in microspore/pollen and tapetum of rice. Plant and Cell Physiology. 2008;49:1429-1450
  39. 39. Yuan TL, Huang WJ, He J, Zhang D, Tang WH. Stage-specific gene profiling of germinal cells helps delineate the mitosis/meiosis transition. Plant Physiology. 2018;176:1610-1626
  40. 40. Nelms B, Walbot V. Defining the developmental program leading to meiosis in maize. Science. 2019;364:52-56
  41. 41. Zhang X, Li J, Liu A, Zou J, Zhou X, Xiang J, et al. Expression profile in rice panicle: Insights into heat response mechanism at reproductive stage. PLoS One. 2012;7:e49652
  42. 42. Wan X, Wu S, Li Z, Dong Z, An X, Ma B, et al. Maize genic male-sterility genes and their applications in hybrid breeding: Progress and perspectives. Molecular Plant. 2019;12:321-342
  43. 43. Zhang D, Wu S, An X, Xie K, Dong Z, Zhou Y, et al. Construction of a multicontrol sterility system for a maize male-sterile line and hybrid seed production based on the ZmMs7 gene encoding a PHD-finger transcription factor. Plant Biotechnology Journal. 2018;16:459-471
  44. 44. Wang Y, Liu D, Tian Y, Wu S, An X, Dong Z, et al. Map-based cloning, phylogenetic, and microsynteny analyses of ZmMs20 gene regulating male fertility in maize. International Journal of Molecular Sciences. 2019;20:1411
  45. 45. An X, Dong Z, Tian Y, Xie K, Wu S, Zhu T, et al. ZmMs30 encoding a novel GDSL lipase is essential for male fertility and valuable for hybrid breeding in maize. Molecular Plant. 2019;12:343-359
  46. 46. Zhu T, Wu S, Zhang D, Li Z, Xie K, An X, et al. Genome-wide analysis of maize GPAT gene family and cytological characterization and breeding application of ZmMs33/ZmGPAT6 gene. Theoretical and Applied Genetics. 2019;132:2137-2154
  47. 47. Xie K, Wu S, Li Z, Zhou Y, Zhang D, Dong Z, et al. Map-based cloning and characterization of Zea mays male sterility33 (ZmMs33) gene, encoding a glycerol-3-phosphate acyltransferase. Theoretical and Applied Genetics. 2018;131:1363-1378
  48. 48. Kim SS, Grienenberger E, Lallemand B, Colpitts CC, Kim SY, Souza Cde A, et al. LAP6/polyketide synthase A and LAP5/polyketide synthase B encode hydroxyalkyl alpha-pyrone synthases required for pollen development and sporopollenin biosynthesis in Arabidopsis thaliana. Plant Cell. 2010;22:4045-4066
  49. 49. Xu J, Ding Z, Vizcay-Barrena G, Shi J, Liang W, Yuan Z, et al. Aborted microspores acts as a master regulator of pollen wall formation in Arabidopsis. Plant Cell. 2014;26:1544-1556
  50. 50. Jiao QJ, Huang Y, Shen HB. Large-scale mining co-expressed genes in Arabidopsis anther: From pair to group. Computational Biology and Chemistry. 2011;35:62-68
  51. 51. Pearce S, Ferguson A, King J, Wilson ZA. FlowerNet: A gene expression correlation network for anther and pollen development. Plant Physiology. 2015;167:1717-1730
  52. 52. Lin H, Yu J, Pearce SP, Zhang D, Wilson ZA. RiceAntherNet: A gene co-expression network for identifying anther and pollen development genes. The Plant Journal. 2017;92:1076-1091
  53. 53. Wang J, Qiu X, Li Y, Deng Y, Shi T. A transcriptional dynamic network during Arabidopsis thaliana pollen development. BMC Systems Biology. 2011;5(Suppl 3):S8
  54. 54. Jiang Y, Zeng B, Zhao H, Zhang M, Xie S, Lai J. Genome-wide transcription factor gene prediction and their expressional tissue-specificities in maize. Journal of Integrative Plant Biology. 2012;54:616-630
  55. 55. Chen W, Chen Z, Luo F, Liao M, Wei S, Yang Z, et al. RicetissueTFDB: A genome-wide identification of tissue-specific transcription factors in rice. Plant Genome. 2019;12:1-11
  56. 56. Luo Y, Guo Z, Li L. Evolutionary conservation of microRNA regulatory programs in plant flower development. Developmental Biology. 2013;380:133-144
  57. 57. Li X. Next-generation sequencing sheds new light on small RNAs in plant reproductive development. Current Issues in Molecular Biology. 2018;27:143-170
  58. 58. Li ZF, Zhang YC, Chen YQ. miRNAs and lncRNAs in reproductive development. Plant Science. 2015;238:46-52
  59. 59. Ru P, Xu L, Ma H, Huang H. Plant fertility defects induced by the enhanced expression of microRNA167. Cell Research. 2006;16:457-465
  60. 60. Chuck G, Meeley R, Irish E, Sakai H, Hake S. The maize tasselseed4 microRNA controls sex determination and meristem cell fate by targeting Tasselseed6/indeterminate spikelet1. Nature Genetics. 2007;39:1517-1521
  61. 61. Millar AA, 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
  62. 62. Ding Y, Ma Y, Liu N, Xu J, Hu Q , Li Y, et al. MicroRNAs involved in auxin signaling modulate male sterility under high-temperature stress in cotton (Gossypium hirsutum). The Plant Journal. 2017;91:977-994
  63. 63. Field S, Thompson B. Analysis of the maize dicer-like1 mutant, fuzzy tassel, implicates microRNAs in anther maturation and dehiscence. PLoS One. 2016;11:e0146534
  64. 64. Xing S, Salinas M, Hohmann S, Berndtgen R, Huijser P. miR156-targeted and nontargeted SBP-box transcription factors act in concert to secure male fertility in Arabidopsis. Plant Cell. 2010;22:3935-3950
  65. 65. Feng N, Song G, Guan J, Chen K, Jia M, Huang D, et al. Transcriptome profiling of wheat inflorescence development from spikelet initiation to floral patterning identified stage-specific regulatory genes. Plant Physiology. 2017;174:1779-1794
  66. 66. Chen J, Su P, Chen P, Li Q , Yuan X, Liu Z. Insights into the cotton anther development through association analysis of transcriptomic and small RNA sequencing. BMC Plant Biology. 2018;18:154
  67. 67. Srivastava S, Zheng Y, Kudapa H, Jagadeeswaran G, Hivrale V, Varshney RK, et al. High throughput sequencing of small RNA component of leaves and inflorescence revealed conserved and novel miRNAs as well as phasiRNA loci in chickpea. Plant Science. 2015;235:46-57
  68. 68. Wei LQ , Yan LF, Wang T. Deep sequencing on genome-wide scale reveals the unique composition and expression patterns of microRNAs in developing pollen of Oryza sativa. Genome Biology. 2011;12:R53
  69. 69. Omidvar V, Mohorianu I, Dalmay T, Fellner M. Identification of miRNAs with potential roles in regulation of anther development and male-sterility in 7B-1 male-sterile tomato mutant. BMC Genomics. 2015;16:878
  70. 70. Yang X, Zhao Y, Xie D, Sun Y, Zhu X, Esmaeili N, et al. Identification and functional analysis of microRNAs involved in the anther development in cotton genic male sterile line Yu98-8A. International Journal of Molecular Sciences. 2016;17:1677
  71. 71. Wei M, Wei H, Wu M, Song M, Zhang J, Yu J, et al. Comparative expression profiling of miRNA during anther development in genetic male sterile and wild type cotton. BMC Plant Biology. 2013;13:66
  72. 72. Sun L, Sun G, Shi C, Sun D. Transcriptome analysis reveals new microRNAs-mediated pathway involved in anther development in male sterile wheat. BMC Genomics. 2018;19:333
  73. 73. Tang Z, Zhang L, Xu C, Yuan S, Zhang F, Zheng Y, et al. Uncovering small RNA-mediated responses to cold stress in a wheat thermosensitive genic male-sterile line by deep sequencing. Plant Physiology. 2012;159:721-738
  74. 74. Niu SH, Liu C, Yuan HW, Li P, Li Y, Li W. Identification and expression profiles of sRNAs and their biogenesis and action-related genes in male and female cones of Pinus tabuliformis. BMC Genomics. 2015;16:693
  75. 75. Shi J, Chen L, Zheng R, Guan C, Wang Y, Liang W, et al. Comparative phenotype and microRNAome in developing anthers of wild-type and male-sterile Lycium barbarum L. Plant Science. 2018;274:349-359
  76. 76. Jiang J, Lv M, Liang Y, Ma Z, Cao J. Identification of novel and conserved miRNAs involved in pollen development in Brassica campestris ssp. chinensis by high-throughput sequencing and degradome analysis. BMC Genomics. 2014;15:146
  77. 77. Lee RC, Feinbaum RL, Ambros V. The C. elegans heterochronic gene lin-4 encodes small RNAs with antisense complementarity to lin-14. Cell. 1993;75:843-854
  78. 78. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: The rosetta stone of a hidden RNA language? Cell. 2011;146:353-358
  79. 79. Thomson DW, Dinger ME. Endogenous microRNA sponges: Evidence and controversy. Nature Reviews. Genetics. 2016;17:272-283
  80. 80. Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, et al. Target mimicry provides a new mechanism for regulation of microRNA activity. Nature Genetics. 2007;39:1033-1037
  81. 81. Tay Y, Rinn J, Pandolfi PP. The multilayered complexity of ceRNA crosstalk and competition. Nature. 2014;505:344-352
  82. 82. Paschoal AR, Lozada-Chavez I, Domingues DS, Stadler PF. ceRNAs in plants: Computational approaches and associated challenges for target mimic research. Briefings in Bioinformatics. 2018;19:1273-1289
  83. 83. Zhai J, Zhang H, Arikit S, Huang K, Nan GL, Walbot V, et al. Spatiotemporally dynamic, cell-type-dependent premeiotic and meiotic phasiRNAs in maize anthers. Proceedings of the National Academy of Sciences of the United States of America. 2015;112:3146-3151
  84. 84. Li Z, An X, Zhu T, Yan T, Wu S, Tian Y, et al. Discovering and constructing ceRNA-miRNA-target gene regulatory networks during anther development in maize. International Journal of Molecular Sciences. 2019;20:3480

Written By

Xiangyuan Wan and Ziwen Li

Submitted: March 23rd, 2019 Reviewed: July 1st, 2019 Published: July 30th, 2019