MicroRNA and Putative Target Discoveries in Chrysanthemum Polyploidy Breeding

MicroRNAs (miRNAs), around 22 nucleotides (nt) in length, are a class of endogenous and noncoding RNA molecule that play an essential role in plant development, either by suppressing the transcription of target genes at a transcriptional level or inhibiting translation at a posttranscriptional level. To understand the roles of miRNAs and their target genes in chrysanthemum polyploidy breeding, three sRNA libraries of normal and abnormal embryos after hybridization were performed by RNA-Seq. As a result, a total of 170 miRNAs were identified and there are 41 special miRNAs in cross of paternal chromosome doubling, such as miR169b, miR440, and miR528-5p. miR164c and miR159a were highly expressed in a normal embryo at 18 days after pollination, suggesting the regulatory role at the late stage of embryonic development. miR172c was only detected in the normal embryo at 18 days after pollination, which means that miR172c mainly mediates gene expression in postembryonic development and these genes may promote embryo maturation. Other miRNAs, including miR414, miR2661, and miR5021, may regulate the genes participated in pathways of auxin response and energy metabolism; then they regulate the complex embryonic development together.


Introduction
In plants, microRNAs (miRNAs) are a major class of small noncoding RNAs (sRNAs) with 20-22 nt. sRNAs have been identified to control the developmental processes in plants by regulating gene expression [1]. They have the potential to regulate a gene by two ways: (1) posttranscriptional gene silencing (PTGS) by binding to 3 ′ untranslated region (UTR) of messenger RNAs (mRNAs) and repressing the translation of target mRNAs; (2) transcriptional gene silencing (TGS) by epigenetic modifications [2,3]. PTGS is the main strategy used by miRNAs to regulate gene expression in plant development. Plant mature miRNAs are generated from miRNA precursors that are processed by a ribonuclease DICER-LIKE1, which negatively regulates specific target mRNAs [2].
A number of studies have investigated the key regulatory role of miRNAs in a wide range of growth and development processes in plants, including the regulation of embryonic development [4,5]. Plant embryo development includes two stages, embryo morphogenesis and seed maturation. In the model plant Arabidopsis thaliana, embryonic pattern formation has been mainly concerned and studied [5]. Based on the embryo shape, the cell division goes through several stages from preglobular to mature embryo, tightly regulated by multiple genes [6]. As a class of small regulatory RNA, the function of miRNAs during plant embryo development has been considerably reported recently [7]. 28,645 mature miRNAs have been discovered and deposited in the public miRNA database miRBase (Release 21, http://www. mirbase.org/) and hundreds of miRNA target in plant embryo. A large number of Arabidopsis mutants of miRNA biogenesis genes have revealed the crucial roles of miRNA during seed morphogenesis and maturation. Willmann et al. reported the earlier timing of embryo maturation in Arabidopsis mutant for strong alleles of DCL1 (DICER-LIKE1) that are required for miRNA biogenesis and demonstrated the negative regulatory role of specific miRNAs during early embryogenesis and later in embryonic development and that miRNAs are key regulators during seed maturation program [8].
Embryonic miRNAs mediate plant embryo development by regulating transcription factors located in special spatial organization and other key developmental regulators [5]. In plants, miR165/166 is one of the best characterized miRNA families, which regulates five class III homeodomain leucine zipper genes (HD-ZIP IIIs) (PHB, PHV, REV, CNA, and ATHB8) [9,10]. The HD-Zip III gene family regulates apical embryo patterning and organ polarity as well as controls shoot and root apical meristem (SAM and RAM) formation in embryogenesis [11][12][13]. The family of miR160/miR167 regulates the target mRNAs of auxin response factors (ARFs) associated with auxin homeostasis [14,15]. Auxin response is an important signaling pathway during embryonic pattern formation, embryo development, and seed maturation [16].
Chrysanthemum (Chrysanthemum morifolium) is an economically important flower around the world, with the increase of chrysanthemum consumption, breeders are driven to improve the cultivars' traits, such as color, size, shape, and tolerance [17]. Artificial distant hybridization is one of the most effective methods to improve and create new cultivars. However, embryo abortion commonly happens during hybridization [18,19]. Although previous studies analyzed many different reasons leading to plant embryo abortion, including maternal genotypes [20], parent ploidy [21,22], and gene regulation [6,23], the molecular mechanisms regulating embryo development are poorly understood. Zhang et al. revealed the gene, protein, and miRNA change in the stage of chrysanthemum embryo development [19,24] and provided numerous information of embryo abortion in chrysanthemum hybridization breeding. Evidence supporting a major role for chromosome doubling in overcoming chrysanthemum embryo abortion has been obtained from some studies [22,25], but the functions of miRNA remain largely uncharacterized.
Next generation sequencing (NGS) approaches have been employed to identify individual miRNAs in various samples, and bioinformatics analyses have offered the technical support to predict the miRNA targets [26]. In the present study, we sequenced three small RNA libraries from chrysanthemum embryo in cross C. morifolium "Yuhualuoying" × tetraploid C. nankingense and identified the key miRNAs and targets that may facilitate embryo development in cross of paternal chromosome doubling.
2.2. RNA Extraction, Small RNA Library Preparation, and Illumina Sequencing. Total RNA was extracted with TRIzol reagent (Takara Bio Inc., Otsu, Japan) according to the manufacturer's protocol. Small RNAs with 18-30 nt fragments were enriched by 15% denaturing polyacrylamide gel electrophoresis. After purification, they were ligated to 5 ′ and 3′ adaptors and reversed transcribed into cDNA by reverse transcription-PCR (RT-PCR). Three small RNA libraries were constructed and sequenced using the Illumina HiSeq™ 2000 by the Beijing Genomics Institute (BGI) (Shenzhen, Guangdong Province, China).

Conserved and Novel miRNA Prediction.
Firstly, the data cleaning analysis was performed by getting rid of low-quality reads, reads with 5 ′ primer contaminants and poly A, reads without 3 ′ primer and the insert tag, and reads shorter than 18 nt. The small RNA tags with miRNA, rRNA, snRNA, snoRNA, and tRNA were annotated by aligning to GenBank (http://www.ncbi.nlm.nih.gov/genbank/) and Rfam database (http://rfam.xfam.org/) using all clean reads of 18~30 nt. The number and proportion of each type of sRNAs were calculated in three libraries. Then, to identify the miRNAs in chrysanthemum embryo, miRBase 19.0 (http://www. mirbase.org/) was used to search the conserved miRNAs by BlantN. Only 90% matched sequences were considered to be conserved miRNAs. To allow the unambiguous mapping of small RNAs to annotations, the priority rule was followed: rRNA (in which GenBank > Rfam) > known miRNA > repeat > exon > intron. Finally, the novel miRNAs were predicted using the Mireap software (https://sourceforge.net/ projects/mireap/); here, the chrysanthemum embryo transcriptome library obtained from the same sample with the present study (the NCBI accession number PRJNA315793) was used as the reference database. In the chrysanthemum embryo, a total of 99,119 unigenes were assembled with a mean length of 550-580 nt, which were used for prediction in the present study.

Target Prediction and Functional
Annotation for miRNAs. The potential target genes of the known miRNAs were predicted by the web tool psRNATarget (http:// plantgrn.noble.org/psRNATarget/) with parameters suggested by Allen et al. [27]. The target genes were identified in chrysanthemum embryo transcriptome dataset [22], and the function of these potential target genes was annotated using the two protein databases, Gene Ontology (GO) (http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/).

Differential Expression of Known miRNA and Their
Targets. To find out the differentially expressed miRNAs in three samples, we normalized the expression of miRNA in three samples and obtained the expression of transcript per million (TPM). Normalized expression = actual miRNA count/total count of clean reads * 1000000. If the miRNA count was zero, it was revised to 0.01 for analysis of differential expression. The fold-change of log 2 (sample 1/sample 2) and p value from the normalized counts were calculated to determine significant expression changes. Finally, those miR-NAs with fold-change > 1 and p value <0.05 were considered to be differentially expressed in the two samples. Heatmap is a visual tool reflecting the expression of differential miRNAs. In the present study, the OmicShare tool, a free online platform for data analysis (http://www.omicshare.com/tools), was plotted. Based on the transcriptome library, we searched the expression pattern of these target genes regulated by the differentially expressed miRNAs.
2.6. Real-Time Quantitative PCR (qRT-PCR) Validation of miRNAs and Target Genes. We randomly selected 12 miRNAs and 9 target genes to validate the reliability of sRNA sequencing by qRT-PCR. RNA samples used for sequencing were reverse transcribed using PrimeScript miRNA qPCR starter kit ver 2.0 (Takara Bio, Dalian, China). qRT-PCR was performed using the SYBR Premix EX Taq Kit (Takara, Dalian, China), and the PCR amplification was done as described by Zhang et al. [24]. Three biological replicates were performed for each sample, and relative expression levels were calculated by the 2 −△△CT method. The chrysanthemum gene EF1α (elongation factor 1a) (GenBank accession number KF305681) was used as a reference, which is stably expressed in chrysanthemum [19]. Special primers were designed using PRIMER3 RELEASE 2.3.4. All of the primers were shown in Table S1.

Sequencing Analysis of sRNA.
To explore the regulation of miRNAs in chrysanthemum embryo development when paternal chromosomes were doubled, three cDNA libraries of small RNAs were constructed, which were named NE12, NE18, and AE18. All clean reads were obtained by filtering the low-quality sequences, adapter sequences, and poly-A sequences shorter than 18 nt, which altogether resulted in more than 99.79% of raw reads. When these small RNA tags were mapped to genome, about 20% of reads in each sample were matched. The unique sRNAs matched to genome in the three libraries were 3,864,037, 3,780,667, and 3,734,164, accounting for 9.62%, 10.04%, and 10.29% of all unique sRNAs (Table 1).
When these unique sRNAs were aligned to the GenBank and miRBase, many types of sRNAs were identified, including miRNA, rRNA, snRNA, snoRNA, and tRNA, but the vast majority of sequences were unannotated. The length of different types of sRNAs was discrepant; since in chrysanthemum embryo, the most abundant classes of sRNA showed the length of 24 nt (dominant siRNA), then 21 nt (mainly miRNAs), and 22 nt ( Figure 1). In the present study, the unique miRNAs in each library were taken into consideration for subsequent analysis, and the proportion of miRNAs was 0.26%, 0.24%, and 0.26% in libraries of NE12, NE18, and AE18 (Table 1).

miRNAs and Target Genes Identified in Three Libraries.
In all samples, a total of 170 conversed miRNAs were identified in miRBase, and 130, 131, and 132 miRNAs were expressed in NE12, NE18, and AE18, respectively. 100 miRNAs (accounting for 58.8%) were detected in three samples; however, some were in two different samples, such as 10 miRNAs in NE12 and NE18, 7 miRNAs in NE18 and AE18, and 6 miRNAs in NE12 and AE18. Venn diagram ( Figure 2) presented the quantity distribution of conserved miRNAs in chrysanthemum embryo.
Expression level is an important feature to explain the regulation function of miRNAs, and they commonly varied greatly in different samples. Here, the expression of 170 miRNAs with sequences was normalized and analyzed. The most abundant miRNAs in three samples were miR156a, miR157a, and miR166a, with the expression more than a thousand (Table S2). However, some miRNAs highly expressed in a particular sample, such as the expression of miR398b-5p, are 3524 in AE18 but neither in NE12 nor NE18, miR5721 only expressed in NE18, and miR5662 only in NE12.
To study the biological function of miRNA in chrysanthemum embryo development, the sequences of target mRNA were paid attention. A total of 770 target genes were identified and regulated by 88 miRNAs (51.76% of all miR-NAs), and miRNA414 had the most targets (347 unigenes), followed were miR5293 (61 targets) and miR5021 (44 targets) (Table S3). 34 target genes were regulated by two miRNAs, such as CL1999. Contig2 was targeted by miR165a and miR166a; unigene10412 was the target of miR156a and miR157a (Table S3).

Functional Annotation of Target
Genes. Gene Ontology (GO) is a standardized classification system of gene used to describe the characteristics of genes and gene products in organisms. The result showed that 517 target genes were classified into three gene ontology categories: biological processes, cellular components, and molecular functions ( Figure 3). Of the 40 functional categories, the number of genes in each sample is different, but in most of categories, NE12 had the largest number but least in AE18. In the third ontology of molecular function, some target genes were expressed in a specific sample, such as the genes related to "electron carrier activity" only expressed in NE12, two genes in the category of "enzyme regulator activity" expressed in NE18, but there was no gene regulated "molecular transducer activity" in NE18 (Figure 3).
To understand the biological function of target genes, KEGG analysis, as with GO, was also used to analyze candidate targets. In NE12, a total of 337 target genes had the biological function on 207 KEGG pathways, but there were less target genes and pathways in NE18 and AE18. 211 target genes with 118 pathway annotation in NE18 and 179 annotated target genes with 113 pathways in AE18 (Table S4). Some target genes were annotated in dozens of pathways that are only in library of NE12, including some energy metabolism pathways, such as ko00280 (valine, leucine, and isoleucine degradation), ko00310 (lysine degradation), and ko00410 (beta-alanine metabolism).

Differentially Expressed miRNAs in Chrysanthemum
Embryo. The aim of this study is to analyze the different miR-NAs possibly involved in chrysanthemum embryo development, so we identified 112 differentially expressed miRNAs with at least 1.5-fold after standardized expression. The heatmap showed the expression pattern of these miRNAs ( Figure 4). They were assembled in three groups depending on the expression trend. In group NE18/AE18, the number of upregulated miRNAs was almost the same as downregulated, and similar quantity distribution occurred in NE12/ NE18. However, the difference is in NE12/AE18, in which only 1/3 of the miRNAs were upregulated in NE12 and means more negative genes in AE18.

Characteristics of Target Genes in Chrysanthemum
Embryo Development. miRNAs regulate the plant development by mediating the expression of target genes. Apart from those redundancy regulated by several miRNAs, a total of 770 target genes were identified and annotated by chrysanthemum transcriptome dataset. The expression level and annotation were presented in Table S5. Most of them were regulated by one miRNA, and some were negatively regulated by miRNAs, such as the transcription factor MYB11 (Unigene2183) highly expressed in NE12, but the targeted miR858b was the lowest expression level in it. Another transcription factor WRKY48 (Unigene25838) was expressed to be the highest in AE18 and lowest in NE12, showing the negative regulation by miR414 ( Figure 5, Tables S2 and S5). Several targets were regulated by two miRNAs, and there was no difference in the level of expression among samples. For example, unigene21756 was the target of miR5215 and miR5373, and its FPKM is near the three samples, without negative regulation by these two miRNAs.
In order to study the target genes that regulate the embryo development, some differentially expressed target genes (the fold of FPKM > 1.5) were selected according to their Nr/Nt annotation, which may be involve in chrysanthemum embryo development (Table 2). They contained some transcription factor, genes related to energy metabolism and protein synthesis, and some uncharacterized protein.
The expression of transcription factor was various; WRKY and NAC were highest in AE18, but MYB was lowest. Some of the genes associated with auxin and ATP synthesis were downregulated in AE18 (Table 2 and Figure 5).
3.6. Validation of qRT-PCR. qRT-PCR is an efficient and accurate way to examine the result of RNA-Seq. In the present study, a total of 12 miRNAs were double tested by qRT-PCR, suggesting the high compatibility of the expression pattern between the two methods ( Figure 6). Seven of them had the highest expression in AE18, such as miR169b, miR159a, and miR858b. miR167c-3p was verified with the highest expression in NE18 by two methods. miR414 that regulated the most target genes with various expression levels was expressed dominantly in NE12 ( Figure 6).

Identified miRNAs in Chrysanthemum
Embryo. miRNAs can regulate the plant embryonic development, which important and diverse roles have been studied in various species; however, the specific function of individual miRNA is still uncharacterized for stage-specific embryo [5]. Next-  Figure 2: Venn diagram presented the quantity distribution of conserved miRNAs in chrysanthemum embryos. There were 130, 131, and 132 miRNAs expressed in NE12, NE18, and AE18, respectively. In the middle, 100 means that the number of miRNAs was detected in three samples. NE12, NE18, and AE18 mean the normal embryo at 12 DAP, normal embryo at 18 DAP, and abnormal embryo at 18 DAP, respectively. generation sequencing makes it easier to identify individual miRNA families from different organs or treated plants. In Pinus taeda [15] and Brassica napus [28], miRNAs were identified from zygotic embryos at late developmental stages. More and more researchers have provided abundant evidences that miRNAs are required for the majority of embryonic cell differentiation and development in Arabidopsis [29,30]. In chrysanthemum cross breeding, embryonic development is a crucial stage for seed formation, but embryo abortion is prevalent in chrysanthemum distant hybridization [31,32]. Previous studies explored this barrier at the level of cell structure, gene expression, and miRNA regulation [19,24], in which 227 miRNAs were identified in hybrid embryos from cross C. morifolium and diploid C. nankingense. Because of the importance of chromosome doubling in cross breeding, in the present study, we performed the hybridization using C. morifolium and tetraploid C. nankingense; 179 miRNAs were identified in three embryo samples (Table S2). Compared with two crosses, less miRNAs were identified after paternal chromosome doubling; as a result, 135 miRNAs were identified simultaneously in two crosses, 44 new miRNAs in the present cross, such as miR169b, miR440, and miR528-5p, but 92 miRNAs expressed in the previous cross were not detected here, such as miR172a, miR172b, and miR391 [24]. These similarities and differences suggest that most of the identified miRNAs in two crosses are the main factor regulating embryonic development in chrysanthemum and those expressed in specific cross may regulate the development depending on whether the male chromosome doubled. To some extent, it implicated that chromosome doubling has an effect on embryo development in distant hybridization by regulation of miRNAs.

The Role of miRNA Chrysanthemum Embryo.
In plant embryo, miRNAs can mediate their downstream targets and regulate the expression of some transcription factors or other key developmental regulators [5]. The overexpression of miRNAs and their target genes have allowed assignment of developmental roles in embryonic, vegetative, and floral organ boundary formation [33,34], such as miR159, miR164, miR165/166, miR172, and miR319 families. It has demonstrated that miR164 directly regulates NAC domain genes for making the function of normal plant morphogenesis and normal embryonic [34]. In the present study, there was no miR164 detected, but miR164c was expressed in three samples and higher expression level in abnormal embryo (Table S2) Figure 4: Heatmap showed the expression pattern of differentially expressed miRNAs in chrysanthemum embryos. The bar represents the ratio of expression levels of miRNAs between two samples. Three ratios of NE18/AE18, NE12/NE18, and NE12/AE18 indicated the expression pattern of each miRNA in three samples. Red rectangle means upregulation, and green means downregulation. All information for each miRNA list can be found in Table S2. NE12, NE18, and AE18 mean the normal embryo at 12 DAP, normal embryo at 18 DAP, and abnormal embryo at 18 DAP, respectively. diploid C. nankingense [24], suggesting that miR164c may have the important regulatory role at the late stage of embryonic development, and the lower expression level of miR164c in normal embryo at 18 DAP may be beneficial to chrysanthemum embryo development normally. Similar situation happened in miR159a; there was no miR159 detected in chrysanthemum embryo, but the miR159a was differentially expressed in normal and abnormal embryo at 18 DAP, suggesting the connected function with miR164a during chrysanthemum embryonic development. miR172 showed the regulatory function for early flowering and floral organ identity defects by regulating AP2 and TOE [35]. miR172c was detected in this study, and the expression was zero in samples of NE12 and AE18, which means that miR172c mainly mediates gene expression in postembryonic development and these genes may promote embryo maturation. DCL1 was a key gene for embryo development by embryo lethality, which was regulated by miR163 [36,37]. However, there was no miR163 identified in chrysanthemum embryo whether the cross has paternal chromosome doubling or not. This result exhibits that DCL1 was not regulated by miR163 and participated in embryonic lethality in chrysanthemum.

miRNA-Mediated Target Genes in Chrysanthemum
Embryo. Research in the last decade has demonstrated that miRNAs make the crucial roles during plant embryogenesis by regulating various genes and pathways [38]. It contains the process of spatial control of differentiation, regulation of auxin responses, and temporal control of differentiation [5]. A study on miR160-resistant ARF17 transgenes showed the defected cotyledons, suggesting that miR160 negatively regulated genes involved in auxin signaling that is critical for proper development of the embryo and cotyledons [39]. Auxin response is a critical biological pathway for embryonic development, and it has also been reported in chrysanthemum [22]. In the present study, miR160 was identified in NE12 and NE18, highly expressed in NE18 (Table S2). The result of target gene prediction showed that the targets of miR160 were auxin response factor with differential  [19], and the value of FPKM was shown on the right of y-axis. NE12, NE18, and AE18 mean the normal embryo at 12 DAP, normal embryo at 18 DAP, and abnormal embryo at 18 DAP, respectively.   Table 2). Transcriptome provided the evidence of importance of energy synthesis for normal embryo development [22]; here, these identified target genes were regulated by miRNAs, such as miR414, miR2661, and miR5021, also support the significance of energy metabolism for chrysanthemum embryo development.

Conclusion
Polyploid breeding will pay more attention from genomic research in the future as rapid advances in the next generation sequencing technology, which makes unprecedented opportunities to explore and understand the regulatory of genomic or transcriptomic changes. As a critical regulatory factor, the function of miRNAs has attracted a lot of attention during plant growth and development. In chrysanthemum distant hybridization, breeders always faced the barriers existed in embryo development. The present study provided some explanation about the embryo abortion and the miRNAs related to embryo development even their target genes. We propose that late embryonic miR-NAs, especially miR164a, regulate NAC transcription factor and thereby affect the embryonic development. miR160 mediated the auxin response, and miR414, miR2661, and miR5021 regulate the genes involved in energy metabolism; together, they regulate the embryo development in chrysanthemum hybridization.