Skip to main content

Small RNA sequencing provides candidate miRNA-target pairs for revealing the mechanism of apomixis in Zanthoxylum bungeanum



Apomixis is a form of asexual reproduction that produces offspring without the need for combining male and female gametes, and the offspring have the same genetic makeup as the mother. Therefore, apomixis technology has great application potential in plant breeding. To identify the apomixis types and critical period, embryonic development at different flower development stages of Zanthoxylum bungeanum was observed by cytology.


The results show that the S3 stage is the critical period of apomixis, during which the nucellar cells develop into an adventitious primordial embryo. Cytological observations showed that the type of apomixis in Z. bungeanum is sporophytic apomixis. Furthermore, miRNA sequencing, miRNA-target gene interaction, dual luciferase reporter assay, and RT-qPCR verification were used to reveal the dynamic regulation of miRNA-target pairs in Z. bungeanum apomixis. The miRNA sequencing identified 96 mature miRNAs, of which 40 were known and 56 were novel. Additionally, 29 differentially expressed miRNAs were screened according to the miRNAs expression levels at the different developmental stages. Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses showed that the target genes of the differentially expressed miRNAs were mainly enriched in plant hormone signal transduction, RNA biosynthetic process, and response to hormone pathways.


During the critical period of apomictic embryonic development, miR172c significantly reduces the expression levels of TOE3 and APETALA 2 (AP2) genes, thereby upregulating the expression of the AGAMOUS gene. A molecular regulation model of miRNA-target pairs was constructed based on their interactions and expression patterns to further understand the role of miRNA-target pairs in apomixis. Our data suggest that miR172c may regulates AGAMOUS expression by inhibiting TOE3 in the critical period of apomixis.

Peer Review reports


The most common form of reproduction in flowering plants is double fertilization. Here, the male and female gametophytes combine to form a new plant [1]. However, another mode of plant reproduction can produce offspring without sexual reproduction, known as apomixis. This is a reproductive mode that directly produces offspring, without combining the male and female gametes. The genetic makeup of the offspring in apomictic species is completely consistent with that of the female parent, and so retains any superior traits the female parent might possess. Hence, apomixis has considerable application potential in crop genetic breeding.

There are three types of apomixis, depending on the source of the embryo: diplospory, apospory, and sporophyte apomixis. Diplospory is when an egg cell in the unreduced diploid embryo sac develops into an embryo [2]. This type can be found in plants such as fleabanes [3]. Apospory is the reproductive process in which the nucellus cells form an unreduced embryo sac through mitosis, and the initial cells develop into the embryo. This type occurs in plants such as bahiagrass [4]. Both these forms of reproduction belong to gametophytic apomixis. Sporophyte apomixis is a process in which an adventitious polyembryony is formed by the mitotic division of nucleated cells. This type can be found in plants such as Cyrtogonellum caducum [5,6,7]. In such cases, the endosperm may not form autonomously by fertilization during the process of apomixis, and pollen and polar nuclei might be required to form a triploid endosperm through pseudogamy [8, 9].

Recent research on the mechanisms of apomixis has made unprecedented progress, and various regulatory factors have been shown to be involved in the apomictic process. Gibberellin-insensitive dwarf protein 1 (BbGID1) is specifically expressed in the nucellus of Brachiaria brizantha. Overexpression of AtGID1 results in the differentiation of Arabidopsis thaliana megaspore mother cell-like cells [10]. Baby boom 1 (BBM1) is a member of the AP2 family of transcription factors expressed in rice sperm cells. Its ectopic expression in egg cells can cause the rice to reproduce apomictically, without fertilization [11]. In addition, a long non-coding RNA (lncRNA), PN_LNC_N13, has differential expression in sexually reproducing and apomictic Paspalum notatum, indicating that lncRNAs might be involved in apomixis regulation [12]. Furthermore, small RNA (sRNA) is an indispensable regulator in apomixis. Argonaute 9-dependent sRNA silencing plays a key role in the development of Arabidopsis ovules, and the sRNA-dependent silencing of plant gametes is required for epigenetic reprogramming of the accompanying cells [13]. In the ovules of apomictic citrus, the expression level of miRN23-5P in monoembryonic varieties was reported to be significantly higher than in polyembryonic varieties, and its expression is negatively correlated with the expression of its regulated target genes [14]. In Eragrostis curvula, miRNA-mRNA interaction analysis showed that miRNA might participate in apomixis by regulating a MADS-box transcription factor gene and a transposon [15]. Therefore, expression analysis of miRNA-target genes would likely be a valuable tool for elucidating the apomictic regulatory networks. Although some progress has been made in the study of apomixis, its full regulatory network remains obscure and incomplete.

Zanthoxylum bungeanum (ZB) is a species of the genus Zanthoxylum, which is widely distributed in Shaanxi, Gansu, Sichuan, Shandong, and Hebei Provinces, and many other places around China [16]. The skin of the ZB fruit is a well-known traditional Chinese seasoning herb. It is widely used in cooking because of its unique numbing properties. It is also used in traditional Chinese medicine to treat various diseases [17]. ZB has apomictic characteristics, and so it is a suitable model in which to study apomixis. In this study, cytological observations, transcriptomics, dual-luciferase reporter assay, RT-qPCR, and miRNA-target gene interaction verification were used to study the molecular mechanism of ZB apomixis. Our results provide references for the study of the apomixis developmental process and mechanism.


Identification of Apomixis in Zanthoxylum bungeanum

The critical period of apomixis embryo development needs to be identified when studying the mechanism of apomixis. To do that, we performed cytological observations on the collected flowers at different developmental stages (Fig. 1a). The results showed that the ovule appears before flowering (S1) and contains a megaspore mother cell (Figs. 1b i and ii). In S2, the megaspore mother cell in the ovule divides to form a binuclear embryo sac (Fig. 1b iii). In S3, the binuclear embryo sac develops into a mature embryo sac, but the embryo sac decomposes and cannot function normally (Fig. 1b iv). Near the embryo sac, the nucellus cells continue to differentiate to form the adventitious primordial embryo (Fig. 1b v). Finally, the adventitious primordial embryo develops into a nucellus embryo in S4 (Fig. 1b vi). A model diagram of adventitious embryogenesis was established to show apomixis in ZB. Based on this diagram, ZB was found to have the characteristics of sporophytic apomixis. The differentiation of nucellar cells into the adventitious primordial embryo at the S3 stage is the most critical period for apomixis.

Fig. 1
figure 1

Cytological observations of apomixis in Zanthoxylum bungeanum. a Different developmental stages of the Zanthoxylum bungeanum fruit. b Cytological observation of apomixis in Zanthoxylum bungeanum. (i-ii) The megaspore mother cell is formed in the S1. (iii) The binuclear embryo sac is formed in the S2. (iv-v) The embryo sac decomposes, and the nucellus cells develop into the adventitious primordial embryo in the S3. (vi) The adventitious primordial embryo develops into a nucellus embryo in the S4. Scale bars = 20 μm. c A model diagram of adventitious embryogenesis in Zanthoxylum bungeanum. ae: adventitious embryo, ap: antipodals, ape: adventitious primordial embryo, bes: binuclear embryo sac, eg: egg, es: embryo sac, mmc: megaspore mother cell, nc: nucellar cell, ne: nucellus embryo, neic: nucellar embryo initial cells, ov: ovule, pn: polar nuclei, sy: synergids

Small RNA sequencing profile

Small RNA sequencing was performed by the Illumina HiSeq 2500 sequencing platform (Illumina, San Diego, CA, USA) on 12 ZB samples at the four apomixis developmental stages. Sequencing obtained 62,810,130 reads and 60,347,610 clean reads, among which high-quality reads accounted for more than 99.95% (Table 1). The Bowtie software was used to analyze the location and distribution of miRNAs on the reference sequence. The analysis results showed that more than 62% of the small RNAs could be located on the reference sequence (Zanthoxylum bungeanum non-reference genome transcriptome spliced transcript as reference sequence). The number of reads mapped to the same strand in the reference sequence exceeded 50%, while the number of reads mapped to the opposite strand of the reference sequence comprised less than 10%.

Table 1 Summary of small RNA sequencing results
Table 2 miRNAs and their target genes

A total of 96 mature miRNAs were identified, of which 40 were known, and 56 were novel. A total of 106 miRNA hairpins were also identified, 49 known and 57 novel (Table S1). Transcripts per million reads (TPM) can represent the miRNAs expression level. According to the log10 (TPM + 1) value of miRNAs, the differentially expressed miRNAs in the different groups were clustered to determine their expression patterns at the four ZB fruit developmental stages (Fig. 2a). The miRNAs expression pattern was directly related to their functions. We primarily analyzed the differentially expressed miRNAs during the critical apomictic embryogenesis period (S3) to study the relationship between miRNAs and apomixis. Additionally, large differences in the TPM distribution of miRNAs were noted between the four ZB fruit developmental stages (Fig. 2b). The length range of animal sRNA is 18–35 nt, while the length range of plant sRNA is 18 ~ 30 nt. The length distribution peak can determine the type of sRNA. For example, the length of miRNA is concentrated in the 21–22 nt range, while the siRNA length is mostly of 24 nt. The miRNAs in the four ZB fruit development stages were 18–30 nt long, with miRNA of 24 nt being the most abundant in all stages (Fig. 2c). The differences in the length distribution of ZB miRNAs offer indirect evidence that there are different regulatory mechanisms at different apomixis stages.

Fig. 2
figure 2

The expression pattern of miRNAs in Zanthoxylum bungeanum apomictic fruit developmental stages. a Cluster analysis of differentially expressed sRNA. b Length distribution of transcripts per million reads (TPM). c Length distribution of the miRNAs

Differentially expressed miRNA analysis in Apomixis

Principal component analysis (PCA) was performed on 12 ZB samples from the four fruit developmental stages According to the miRNAs expression level or TPM. The results showed that PC1 and PC2 explained together 81.0% of the sample differences, and sample clustering into groups showed that they had good repeatability and met the data analysis requirements (Fig. 3a). Analysis of the differentially expressed miRNAs in the experimental groups can screen for the miRNAs that play key roles in the process. Therefore, we used the DESeq software to screen and analyze the miRNAs at the four ZB apomixis stages. We found 29 differentially expressed miRNAs (Fig. 3b). Of these, 15 were known miRNAs, and 14 were novel. An Upset plot was drawn based on the number distribution of 96 mature miRNAs at the different developmental stages (Fig. 3c). The results showed that the number of mature miRNAs detected in S3 was the largest (82), while that in S4 was the smallest (64). Furthermore, the numbers of miRNAs unique to developmental stages S1, S2, S3, and S4 were 6, 2, 5, and 2, respectively. The number of miRNAs shared by the four developmental stages was 50.

Fig. 3
figure 3

Functional enrichment analysis of differentially expressed miRNAs and their target genes at the four Zanthoxylum bungeanum fruit developmental stages. a Principal component analysis of fruits at the four developmental stages. b Venn diagram of differentially expressed miRNAs at the different developmental stages. c Quantitative distribution of mature miRNAs at the different developmental stages. The bottom panel shows the number distribution of miRNAs at different developmental stages, and the connected black dots represent the number of miRNAs shared by different samples. d Gene Ontology (GO) enrichment analysis of the differentially expressed miRNAs target genes

miRNA participates in various life processes by regulating the expression of target genes or inhibiting their translation. Therefore, it is necessary to determine the miRNA functions through their interactions with target genes and their expression levels. We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) to analyze pathways enrichment of differentially expressed genes during apomixis in ZB. The results showed that many differentially expressed genes were associated with RNA transport, protein synthesis on the plant hormone signal transduction (ko04075), RNA transport (ko03013), and steroid biosynthesis (ko00100; Figure S1). Many substances were also synthesized, including fatty acids, starch, sucrose, and folate. In addition, Gene Ontology (GO) enrichment analysis was performed on the target genes of differentially expressed miRNAs (Fig. 3d). The results showed that many differentially expressed genes were enriched in several pathways, including RNA biosynthetic process (GO:0032774), transcription factor activity (GO:0003700), response to hormone (GO:0009725), and more. Among them, 457, 167, and 81 differentially expressed genes were enriched in the RNA biosynthetic process, transcription factor activity, and response to hormone pathways, respectively, indicating that scores of RNA and transcription factors were synthesized, activating the hormone response pathways during ZB apomixis. KEGG and GO enrichment analyses of the differentially expressed miRNAs target genes showed that many genes were enriched in pathways related to hormone synthesis and regulation, indicating hormones that might be involved in ZB apomixis regulation.

Interaction analysis of the miRNAs and their target genes

A comprehensive analysis of miRNA and mRNA expression profiles helped identify functional miRNA-mRNA interaction pairs involved in the apomictic processes (Fig. 4a and Table 2). Eight miRNAs in the four developmental stages were significantly differentially expressed. These were detected by quantitative reverse transcription PCR (RT-qPCR) to verify the validity of the sequencing results (Fig. 5). The results showed that the miRNAs sequencing and RT-qPCR results were essentially consistent, confirming their reliability. Additionally, we predicted the target genes of differentially expressed miRNAs, and the results showed that they were involved in genes and regulatory factors of multiple life processes, including hormone-related ABI5, IAA26, and TOE3, and several transcription factors such as WRKY75, MYB3, and TCP2. Among them, miR172c had two target genes, TOE3 and AP2. The mature miR172c can form complementary pairs with both target genes (Fig. 4b). It had a high expression level in the key ZB apomixis stage (S3), implying that it and its target genes played an important role in apomixis.

Fig. 4
figure 4

miRNAs and their target genes. a Sankey diagram of miRNAs and their target genes. b miR172c precursor, mature body, and binding site with target genes

Fig. 5
figure 5

Relative expression levels of miRNAs and their target genes involved in apomixis by RT-qPCR

The relative expression levels of miRNAs and their target genes were analyzed by RT-qPCR to further determine the interactions between them (Fig. 5). The results showed that during ZB fruit development, the expression levels of the eight differentially expressed miRNAs and their target genes were negatively correlated, preliminarily verifying the mode of interaction between them.

Verification of the interaction between miRNA and its target genes

We constructed miRNA and target gene sequences on pCA1301 and pBI121 vectors, respectively, to demonstrate the interactions between them. These were transferred to tobacco leaves, and the expression of the GUS protein was evaluated to determine whether there was an interaction between the miRNA and its target gene [18]. The 35S::AP2-GUS, 35S::TOE3-GUS, and 35S::miR172c-GUS produced a large amount of GUS protein after injection into tobacco leaves, while 35S::AP2-GUS + 35S::miR172c-GUS and 35S::TOE3-GUS + 35S::miR172c-GUS significantly reduced GUS protein activity (Fig. 6a and b). These results suggest that miR172c could cleave TOE3 and APETALA 2 (AP2) genes and thus block the synthesis of the GUS protein, confirming that TOE3 and AP2 were its target genes.

Fig. 6
figure 6

Verification of the interactions between miR172c and its target genes. a GUS immunohistochemistry. b GUS activity. c The miRNA-target interactions, and their expression patterns during apomixis

A molecular regulation model was constructed based on miR172c-AP2 and miR172c-TOE3 interactions to understand the molecular process in ZB apomixis (Fig. 6c). AGAMOUS belongs to the MADS-Box transcription factor class C genes, and AP2 belongs to the MADS-Box transcription factor class A genes. Class A and C genes have antagonistic effects on each other. AGAMOUS can control the development of stamens and carpels and regulate the activity of floral meristems, playing an important role in apomixis [19, 20]. AGAMOUS expression level gradually increased during fruit development, especially during the critical period of ZB apomictic embryogenesis (S3). The expression level of AP2, on the other hand, gradually decreased during development. The miRNA-target genes interaction experiment showed that miR172c could cleave TOE3, and TOE3 could inhibit AGAMOUS expression (Fig. 6c). Therefore, miR172c can increase AGAMOUS expression abundance by inhibiting TOE3 transcription level during ZB apomixis. Besides, during fruit development, the expression level of the auxin response factor (ARF) increased gradually, suggesting that hormones might play an important role in apomixis.


Apomixis results in the offspring inheriting the complete genotype of the female parent, thus having broad application potential in agricultural heterosis fixation. However, the mechanism of apomixis is still unclear. Clarifying the apomixis embryogenesis process is a prerequisite for studying its mechanism. We analyzed the apomixis embryogenesis process in ZB by cytological observations, concluding that ZB was the sporophytic apomixis type and that S3 was a critical period of apomixis. In this study, miRNA sequencing, miRNA-target gene interactions, dual luciferase reporter assay, and RT-qPCR verification were used to reveal the dynamic regulation of miRNA-target pairs in ZB apomixis. The miRNAs are involved in various biological processes and are among the most important factor groups regulating organisms [21,22,23]. Many differentially expressed miRNAs were detected in ZB apomixis, indicating that they probably play an important role in the process. Studying their interactions with their target genes during apomixis could identify candidate genes and new ideas to help elucidate the apomixis mechanism. Using miRNA sequencing, miRNA-target gene interaction experiments, and RT-qPCR verification, we showed that miR172 inhibits the activity of TOE3 and AP2 during S3, thereby increasing the expression level of AGAMOUS. Because of that, miR172 might be a key regulator of apomixis (Fig. 7). Moreover, miR172 is evolutionarily conserved among angiosperms and is widely involved in plant floral organ development. It controls the flowering time of floral organs by regulating multiple MADS-box transcription factors [24, 25].

Fig. 7
figure 7

miRNA-targets phytohormone regulation model during apomixis in Zanthoxylum bungeanum

The functional analysis of differentially expressed miRNAs and their target genes in apomixis showed that many target genes play an active role in hormone synthesis and regulation. To further understand the possible role of miRNAs in apomixis, based on the analysis results and previous studies, a miRNA interaction regulation model was established. Moreover, numerous studies have shown that hormones can induce plant somatic cells to form embryos and activate plant cell totipotency [26, 27]. Studies have also detected many hormonal changes during seed embryo development [28, 29]. In Brassica napus, jasmonic acid is a tissue-specific hormone involved in gene and protein expression during embryonic development. Abscisic acid (ABA) has been shown to induce cell differentiation during somatic embryogenesis in Cunninghamia lanceolate, thus promoting embryonal mass growth [30]. The auxin-responsive gene, ABA-related gene, and gibberellin-related gene are all regulated by miR393. This miRNA plays an important role in barley hormone synthesis and embryonic development [31]. Based on our analysis, it is inferred that miRNAs regulate the hormone synthesis by interacting with their target genes, thereby activating the totipotency of plant cells that participate in apomixis.

WRKY75 is a target gene of ath-miR166a-5p, which binds and degrades it to block its transcription and inhibit its function. RT-qPCR showed that WRKY75 and ath-miR166a-5p expression trends during ZB apomixis were negatively correlated, and that WRKY75 relative expression level was significantly increased. WRKY75 is a transcription factor that interacts with many hormones [32]. AtWRKY75 can activate the pathways of jasmonic acid and ethylene [33], and promote the synthesis of PHR1 and MYB62, members of the MYB transcription factor family. They inhibit gibberellin synthesis and its signaling pathway [34, 35]. ABA-insensitive 5 protein (ABI5) binds to the embryonic specification element (ESE) of the Dc3 gene promoter and ABA response element (ABRE). By controlling the expression of ABA regulatory genes during seed development, it directly affects ABA synthesis [36]. Moreover, ERF7 can directly inhibit ABA synthesis by interacting with SNL1 (SIN3-LIKE1) and SNL2 in the presence of ethylene [37, 38]. By analyzing the interaction between mRNA and miRNA, we found that ath-miR167d binds to and inhibits ABI, limiting its synthesis. However, our relative expression levels analysis demonstrated that ath-miR167d expression level is low during apomixis, so it could not restrict ABA synthesis at this stage. Furthermore, miRNAs play a role in the auxin regulatory pathway. AVP1 encodes a proton pump that regulates the pH in plant vacuoles. AVP1 and AUX1/IAA participate in gradient regulation and transport of auxin [39, 40]. The expression of AVP1 is restricted by ath-miR858a, and the relative expression level of the two is negatively correlated, as was shown by RT-qPCR analysis. An increase in ath-miR858a expression level in the ZB apomixis process inhibits the function of AVP1. In summary, we suggest that miRNA and its target genes participate in apomixis by regulating hormone synthesis.


We have demonstrated through cytological observation that ZB employs the sporophyte apomixis type and that the S3 stage is a critical period in the process. We performed miRNA sequencing and RT-qPCR to screen candidate miRNAs involved in apomixis. Our results showed that miR172 had a higher expression level during the apomixis S3 period. Furthermore, the dual-luciferase reporter assay verified that miR172 regulates TOE3 and AP2 expression levels. These interact with AGAMOUS to regulate the apomixis process. This study showed evidence for the apomixis process in ZB and provided insights into its mechanism. It is necessary to clarify the roles of the regulatory factors at multiple levels and by multiple means to understand the apomixis mechanism more clearly, and, ultimately, describe its regulatory interaction network.


Plant materials

Plant materials were collected from the ZB Experimental Station, Northwest Agriculture and Forestry University, Fengxian, Shaanxi, China. Healthy, five-year-old ZB (Zanthoxylum bungeanum cv. Hanchengdahongpao) shrubs with uniform growth were selected. ZB is a dioecious plant. Unopened ZB female flowers were bagged and collected at four stages of development: S1, pre-flowering, flowers not to open; S2, mid-flowering, 7 days after flowering starts; S3, young fruit, 15 days after flowering starts; S4, fruit expansion, 30 days after flowering starts. Six-hundred flowers were collected from three ZB trees at each stage. The samples were immediately stored in liquid nitrogen for genetic analysis. Three biological replicates (three ZB trees) were collected at each developmental stage. The flowers or fruits used for cytological observations were stored in formalin-acetic acid-alcohol (FAA) fixative.

miRNAs library construction and sequencing

Plant samples were sent to the Beijing Novogene Genomics Institute for miRNAs sequencing. The HiSeq and MiSeq techniques were used for miRNAs sequencing. The Small RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) was used to construct sequencing libraries. The 12 miRNAs libraries established included three biological replicates at each of the four fruit development stages. The sRNA terminals were directly added with adaptors and reverse-transcribed to synthesize cDNA. PAGE gel electrophoresis was used to separate the target DNA fragments after PCR amplification. The cDNA libraries were obtained by cutting the gel and recovering the DNA. The adaptors and low-quality reads were removed from the raw reads to obtain the clean reads. Clean reads with a length of 18–30 nucleotides were kept for downstream analysis. The clean reads were compared with the sequence in miRBase, and the secondary structures of known miRNA, miRNA sequence, length, abundance, and other information were obtained. miREvo [41] and miRDeep2 [42] were used to predict and analyze novel miRNA. We then predicted the target genes of differentially expressed miRNAs and performed GO and KEGG enrichment analyses on these target genes.

Differentially expressed miRNA analysis

Transcripts per million reads [TPM = (Read count*1,000,000)/ Mapped Reads count] represent the expression level of miRNAs. The DESeq software was used to analyze differentially expressed miRNAs between developmental stages. The screening criteria used were False Discovery Rate (FDR) ≤ 0.05 and │log2 FC│ ≥ 1, where FC is fold change.

Total RNA extraction and cDNA synthesis

The TaKaRa MiniBEST Plant RNA Extraction Kit (TaKaRa, Beijing, China) was used to extract total RNA from the flowers or fruits at the different developmental stages. The first-strand cDNA synthesis for the miRNA was carried out using the Mir-X miRNA First-Strand Synthesis Kit (TaKaRa) according to the manufacturer’s instructions. These were used as templates for the RT-qPCR reactions.

Quantitative real-time PCR (RT-qPCR)

The CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) was used to confirm the relative expression of miRNAs during apomixis. A reaction system of 10 μL contained 5 μL of 2× SYBR Premix Ex Taq II (TaKaRa), 1 μL of cDNA, 1 μL of each of the forward and universal reverse primers, and 2 μL of ddH2O. U6 was used as an endogenous reference gene. Similarly, the relative expression levels of the target genes were detected by RT-qPCR, with ZBUBQ and ZBTIF as reference genes [43]. Primers were designed using Primer Premier 5.0 (Premier, Palo Alto, CA, USA). The primer sequences are shown in Table S3. The RT-qPCR reaction protocol was: 95 °C for 30 s followed by 40 cycles of 94 °C for 5 s, 54 °C for 30 s, and 72 °C for 45 s. Three technical replicates were made for each sample.

Prediction of miRNAs target genes

The psRNATarget online tool ( uses the full-length transcriptome as the target template to predict the target genes of known and novel miRNAs. The specific parameters were set as follows: expectation: 5, penalty for G:U pair: 0.5, penalty for other mismatches: 1, extra weight in seed region: 1.5, seed region: 2–13 nucleotides, mismatches allowed in seed region: 2, HSP size: 19, penalty for opening gap: 2, penalty for extending gap: 0.5, and translation inhibition range: 10–11 nucleotides.

Dual-luciferase reporter assay system

The dual-luciferase reporter assay system was used to verify the interaction between the miRNA and its target genes. The miRNA precursors and their target genes were constructed on pCA1301 and pBI121 vectors, respectively. A GUS (β-glucuronidase) reporter gene was inserted downstream of the pBI121 vector. If the miRNA cleaves the mRNA, the downstream GUS gene cannot be expressed, and thus the GUS protein cannot be detected. The interaction between miRNAs and target genes can thus be verified by detecting GUS protein expression [18].

The GUS staining experiment was carried out on four-week-old tobacco plant leaves (Nicotiana benthamiana). We detected the GUS activity in leaves injected with LB, 35S::AP2-GUS, 35S::miR172c + 35S::GUS, 35S::AP2-GUS + 35S::miR172c, 35S::TOE3-GUS, and 35S::TOE3-GUS + 35S::miR172c. We thoroughly grounded 0.5 g of the liquid nitrogen-frozen tobacco leaves, which were used as experimental material. We then used the GUS histochemical staining Kit (MaoKang, Shanghai, China), following the manufacturer’s instructions, to detect the GUS protein expression in them.

Availability of data and materials

All the data and materials that are required to reproduce these findings can be shared by contacting the corresponding author, Prof. Anzhi Wei (



Abscisic acid

ABI5 :

ABA-Insensitive 5 protein


ABA response element

AP2 :



Auxin response factor


Gibberellin-insensitive dwarf protein 1

BBM1 :

Baby boom 1


Embryonic specification element


Gene Ontology


Kyoto Encyclopedia of Genes and Genomes


Long non-coding RNA


Principal component analysis


small RNA


Transcripts per million reads


Zanthoxylum bungeanum






  1. Raghavan V. Some reflections on double fertilization, from its discovery to the present. New Phytol. 2003;159(3):565–83.

    Article  CAS  Google Scholar 

  2. Hojsgaard D. Apomixis technology: separating the wheat from the chaff. Genes. 2020;11(4):411.

    Article  CAS  PubMed Central  Google Scholar 

  3. Noyes RD. Inheritance of apomeiosis (diplospory) in fleabanes (Erigeron, Asteraceae). Heredity. 2005;94(2):193–8.

    Article  CAS  PubMed  Google Scholar 

  4. Martínez EJ, Urbani MH, Quarin CL, Ortiz JPA. Inheritance of Apospory in Bahiagrass, Paspalum notatum. Hereditas. 2010;135(1):19–25.

    Article  Google Scholar 

  5. Koltunow AM. Apomixis: embryo sacs and embryos formed without meiosis or fertilization in ovules. Plant Cell. 1993;5(10):1425–37.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Gupta P, Shivanna KR, Hym R. Apomixis and Polyembryony in the Guggul plant, Commiphora wightii. Ann Botany. 1996;78(1):67–72.

    Article  Google Scholar 

  7. Liu HM, Dyer RJ, Guo ZY, Zhen M, Li JH, Harald S. The evolutionary dynamics of apomixis in ferns: a case study from polystichoid ferns. J Botany. 2012;2012:1–11 (2012-11-5).

    Article  Google Scholar 

  8. Yao JL, Hu ZCG. Apomixis in Eulaliopsis binata: characterization of reproductive mode and endosperm development. Sex Plant Reprod. 2007;20(3):151–8.

    Article  Google Scholar 

  9. Ogawa D, Johnson SD, Henderson ST, Koltunow AM. Genetic separation of autonomous endosperm formation (AutE) from the two;other components of apomixis in Hieracium. Plant Reprod. 2013;26(2):113–23.

    Article  PubMed  Google Scholar 

  10. Ferreira LG, de Alencar Dusi DM, Irsigler AST, Gomes A, Mendes MA, Colombo L, et al. GID1 expression is associated with ovule development of sexual and apomictic plants. Plant Cell Rep. 2018;37(2):293–306.

  11. Khanday I, Skinner D, Yang B, Mercier R, Sundaresan V. A male-expressed rice embryogenic trigger redirected for asexual propagation through seeds. Nature. 2019;565(7737):91–5.

    Article  CAS  PubMed  Google Scholar 

  12. Ochogavia A, Galla G, Seijo JG, Gonzalez AM, Bellucci M, Pupilli F, et al. Structure, target-specificity and expression of PN_LNC_N13, a long non-coding RNA differentially expressed in apomictic and sexual Paspalum notatum. Plant Mol Biol. 2018;96(1–2):53–67.

  13. Vianey OM, Noé DF, Mario AV, Edgar DA, Daphné A, Daniel G. R Keith S, Martienssen RA, Jean-Philippe VC: control of female gamete formation by a small RNA pathway in Arabidopsis. Nature. 2010;464(7288):628–32.

    Article  Google Scholar 

  14. Long J-M, Liu Z, Wu X-M, Fang Y-N, Jia H-H, Xie Z-Z, et al. Genome-scale mRNA and small RNA transcriptomic insights into initiation of citrus apomixis. J Exp Bot. 2016;67(19):5743–56.

  15. Garbus I, Selva JP, Pasten MC, Bellido AM, Carballo J, Albertini E, et al. Characterization and discovery of miRNA and miRNA targets from apomictic and sexual genotypes of Eragrostis curvula. BMC Genomics. 2019;20(1):839.

  16. Li-Chen Y, Rong L, Jin T, Zi-Tao J. Polyphenolics composition of the leaves of Zanthoxylum bungeanum maxim. Grown in Hebei, China, and their radical scavenging activities. J Agric Food Chem. 2013;61(8):1772–8.

    Article  Google Scholar 

  17. He F, Li D, Wang D, Deng M. Extraction and purification of Quercitrin, Hyperoside, Rutin, and Afzelin from Zanthoxylum Bungeanum maxim leaves using an aqueous two-phase system. J Food Sci. 2016;81(7):C1593–602.

    Article  CAS  PubMed  Google Scholar 

  18. Heng C, Chunxia Y, Sian L, Haoran Q, Ling W, Li-An X, et al. MiRNA-target pairs regulate adventitious rooting in Populus: a functional role for miR167a and its target Auxin response factor 8. Tree Physiol. 2019;39(11):1922–36.

  19. Liu DD, Dong QL, Sun C, Wang QL, You CX, Yao YX, et al. Functional characterization of an apple apomixis-related MhFIE gene in reproduction development. Plant Sci. 2012;185–186(none):111.

  20. Uemura A, Yamaguchi N, Xu Y, Wee WY, Ichihashi Y, Suzuki T, et al. Regulation of floral meristem activity through the interaction of AGAMOUS, SUPERMAN, and CLAVATA3 in Arabidopsis. Plant Reprod. 2018;31(1):89–105.

  21. Lionel N, Patrice D, Florence J, Benedict A, Nihal D, Mark E, et al. A plant miRNA contributes to antibacterial resistance by repressing auxin signaling. Science. 2006;312(5772):436.

  22. Zhao S, Wang X, Yan X, Guo L, Mi X, Xu Q, et al. Revealing of MicroRNA involved regulatory gene networks on Terpenoid biosynthesis in Camellia sinensis in different growing time points. J Agric Food Chem. 2018;66(47):12604–16.

  23. Zhang X, Li K, Xing R, Liu S, Chen X, Yang H, et al. miRNA and mRNA expression profiles reveal insight into chitosan-mediated regulation of plant growth. J Agric Food Chem. 2018;66(15):3810–22.

  24. Huijser P, Schmid M. The control of developmental phase transitions in plants. Development. 2011;138(19):4117–29.

    Article  CAS  PubMed  Google Scholar 

  25. Qian-Hao, Zhu, Chris, a., Helliwell: regulation of flowering time and floral patterning by miR172. J Exp Bot 2010, 62(2):487–495, DOI:

  26. Vasil IK, Vasil V. Totipotency and embryogenesis in plant cell and tissue cultures. Vitro. 1972;8(3):117–27.

    Article  CAS  Google Scholar 

  27. Oinam GS, Kothari SL. Totipotency of coleoptile tissue in indica rice ( Oryza sativa L. cv. ch 1039). Plant Cell Rep. 1995;14(4):245.

    Article  CAS  Google Scholar 

  28. Miransari M, Smith DL. Plant hormones and seed germination. Environ Exp Botany. 2014;99(3):110–21.

    Article  CAS  Google Scholar 

  29. Luckwill LC. Studies of fruit development in relation to plant hormones: I. hormone production by the developing apple seed in relation to fruit drop. J Pomol Horticult Sci. 2015;28(1):14–24.

    Article  Google Scholar 

  30. Zhou X, Zheng R, Liu G, Xu Y, Zhou Y, Laux T, et al. Desiccation Treatment and Endogenous IAA Levels Are Key Factors Influencing High Frequency Somatic Embryogenesis in Cunninghamia lanceolata (Lamb.) Hook. Front Plant Sci. 2017;8:2054.

  31. Bai B, Shi B, Hou N, Cao Y, Meng Y, Bian H, et al. microRNAs participate in gene expression regulation and phytohormone cross-talk in barley embryo during seed development and germination. BMC Plant Biol. 2017;17(1):150.

  32. Baek D, Chun HJ, Yun DJ, Kim MC. Cross-talk between phosphate starvation and other environmental stress signaling pathways in plants. Mol Cells. 2017;40(10):697–705.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Chen X, Liu J, Lin G, Wang A, Wang Z, Lu G. Overexpression of AtWRKY28 and AtWRKY75 in Arabidopsis enhances resistance to oxalic acid and Sclerotinia sclerotiorum. Plant Cell Rep. 2013;32(10):1589–99.

    Article  CAS  PubMed  Google Scholar 

  34. Doerks T, Copley RR, Schultz J, Ponting CP, Bork P. Systematic identification of novel protein domain families associated with nuclear functions. Genome Res. 2002;12(1):47–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Devaiah BN, Madhuvanthi R, Karthikeyan AS, Raghothama KG. Phosphate starvation responses and gibberellic acid biosynthesis are regulated by the MYB62 transcription factor in Arabidopsis. Mol Plant. 2009;2(1):43–58.

    Article  CAS  PubMed  Google Scholar 

  36. Chak RKF, Thomas TL, Quatrano RS, Rock CD. The genes ABI1 and ABI2 are involved in abscisic acid- and drought-inducible expression of the Daucus carota L. Dc3 promoter in guard cells of transgenic Arabidopsis thaliana (L.) Heynh. Planta. 2000;210(6):875.

    Article  CAS  Google Scholar 

  37. Liu W, Wang Y, Gao C. The ethylene response factor (ERF) genes from Tamarix hispida respond to salt, drought and ABA treatment. Trees. 2014;28(2):317–27.

    Article  CAS  Google Scholar 

  38. Zhi W, Hong C, Yongzhen S, Xiaoying L, Fengying C, Annaick C, et al. Arabidopsis paired amphipathic helix proteins SNL1 and SNL2 redundantly regulate primary seed dormancy via abscisic acid-ethylene antagonism mediated by histone deacetylation. Plant Cell. 2013;25(1):149–66.

  39. Liu K, Feng S, Pan Y, Zhong J, Chen Y, Yuan C, et al. Transcriptome analysis and identification of genes associated with floral transition and flower development in sugar apple (Annona squamosa L.). Front Plant Sci. 2016;7:1695.

  40. Wei L, Xiao Z, Chao F, Jiang N, Meng X, Xu X. Cloning and characterization of a Flavonol synthase gene from Litchi chinensis and its variation among Litchi cultivars with different fruit maturation periods. Front Plant Sci. 2018;9:567.

    Article  Google Scholar 

  41. Wen M, Shen Y, Shi S. Tang T: miREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. Bmc Bioinform. 2012;13(1):140.

    Article  CAS  Google Scholar 

  42. Mackowiak SD. Identification of novel and known miRNAs in deep-sequencing data with miRDeep2. Curr Protoc Bioinform. 2011;12(12):12.

    Google Scholar 

  43. Fei X, Shi Q, Yang T, Fei Z, Wei A. Expression stabilities of ten candidate reference genes for RT-qPCR in Zanthoxylum bungeanum Maxim. Molecules. 2018;23(4):802.

    Article  Google Scholar 

Download references


The authors thank the experimental technical support provided by Dr. Huang Dong.


This study was financially supported by the National Key Research and Development Program Project Funding (2018YFD1000605). The funders had no role in the experimental design, data analysis, decision to publish, or preparation of the manuscript.

Author information

Authors and Affiliations



AW conceived the project; X.F. designed and performed the experiment; XF wrote the paper. All authors (XF, YL, YQ, SW, HH, YH, YM, and AW) discussed the results and commented on the manuscript. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Anzhi Wei.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fei, X., Lei, Y., Qi, Y. et al. Small RNA sequencing provides candidate miRNA-target pairs for revealing the mechanism of apomixis in Zanthoxylum bungeanum. BMC Plant Biol 21, 178 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: