Development of Incompletely Fused Carpels in Maize Ovary Revealed by miRNA, Target Gene and Phytohormone Analysis

Although the molecular basis of carpel fusion in maize ovary development remains largely unknown, increasing evidence suggests a critical role of microRNAs (miRNAs). In this study, a combination of miRNA sequencing, degradome and physiological analyses was used to characterize carpel fusion development in maize ovaries showing incompletely (IFC) and completely fused carpels (CFC). A total of 162 known miRNAs distributed across 33 families were identified, of which 20 were differentially expressed. In addition, 53 miRNA candidates were identified, of which 10 were differentially expressed in the IFC and CFC ovaries. In degradome analysis, a total of 113 and 11 target genes were predicted for the known and novel miRNAs, respectively. Moreover, 24 (60%) target genes of the differentially expressed known miRNAs were found to code transcription factors, including auxin response factor (ARF), TB1-CYC-PCFs (TCP), APETALA2 (AP2), growth regulating factor (GRF), MYB, NAC, and NF-YA, all of which have been shown to play a role in carpel fusion development. Correlation analysis of these differentially expressed known miRNAs and their targets with phytohormone signals revealed significant correlations with at least one phytohormone signal, the main regulator of carpel fusion development. These results suggest that incomplete carpel fusion is partly the result of differential expression of certain miRNAs and their targets. Overall, these findings improve our knowledge of the effect of miRNA regulation on target expression, providing a useful resource for further analysis of the interactions between miRNAs, target genes and phytohormones during carpel fusion development in maize.


INTRODUCTION
Maize (Zea mays L.) is one of the most productive cereals in the world, and is widely used as a model organism for genetic research in cultivated crop plants (Wallace et al., 2014). Maize is monocotyledon, with two different types of inflorescence in a single plant: the tassel (male inflorescence) and the ear (female inflorescence) (McSteen et al., 2000). Development of the pistil in female inflorescences is crucial to maize kernel development and grain yield formation. The inner whorl of the floret, the maize pistil is derived from fusion of three carpels, and is composed of a silk and an ovary (McSteen et al., 2000). Incomplete carpel fusion results in a hole in the pericarp of the kernel. Without the protection of the intact pericarp, the kernel is highly susceptible to pathogen infection, threatening food safety (Duncan and Howard, 2010). In addition, these kernels rot more easily after sowing, and thus, complete carpel fusion is essential for food security as well as seed vigor quality. However, despite this, few studies have examined carpel fusion development in maize. Although several genes associated with carpel organ identity have been identified in maize using reverse genetics and molecular studies (Mena et al., 1996), knowledge of carpel fusion remains limited and regulation of the genes involved is unknown.
Incomplete carpel fusion is thought to be regulated by transcription factors (TFs) and hormonal balance (Reyes-Olalde et al., 2013;Marsch-Martínez and de Folter, 2016). Mutations in TFs such as APETALA2 (AP2), CUP-SHAPED COTYLEDON2 (CUC2) and SPATULA (SPT) can lead to partially or almost completely unfused carpels in the gynoecium of Arabidopsis (Ripoll et al., 2011;Nahar et al., 2012). The Arabidopsis mutant HECATE has a low auxin (IAA) content and presents a carpel fusion-deficient phenotype in the gynoecium (Schuster et al., 2015). In line with this, a high concentration of IAA was found to be important for apical fusion of the two carpels in the stigma of the Arabidopsis gynoecium (Larsson et al., 2013;Sehra and Franks, 2015). Moreover, plant hormones are known to affect transcriptional regulation through hormone sensing, synthesis and transport (Marsch-Martínez and de Folter, 2016). To understand carpel fusion in maize, the mechanisms regulating carpel fusion and limiting intact ovary wall formation need to be determined, based not only on genetic research and conventional physiological studies, but also at the level of post-transcriptional regulation.
MicroRNAs (miRNAs) are small non-coding RNA molecules that negatively regulate gene expression mainly through mRNA cleavage or translational inhibition, or DNA methylation of miRNA genes (Voinnet, 2009). MiRNAs, which are generated from single-strand RNA precursors able to form hairpin structures, have been widely studied as essential regulators of diverse aspects of plant development (Larue et al., 2009), including flower development. For example, maize ts4, which encodes miR172 and targets the AP2 homolog indeterminate spikelet 1, has two tandem AP2 domains and plays an important role in regulating maize inflorescence development (Chuck et al., 2007). In ts4 mutants, the floret of the male inflorescence fails to form stamens and develops unfused carpels (Chuck et al., 2007), suggesting a role of miR172 in regulation of carpel fusion development in maize female inflorescences. Furthermore, in Arabidopsis, miR164 has been shown to target a subset of NAC TFs that includes CUC1 and CUC2, which contribute to organ boundary formation including carpel marginal tissue development (Nahar et al., 2012). Moreover, miR160, which is believed to target mRNA coding ARF DNA-binding protein, is thought to be involved in female and male flower development in poplar through regulation of auxin signaling (Song et al., 2013). Thus, miRNAs have great potential as a tool for elucidating floral organ development in maize and other plant species.
In this study, we used high-throughput sequencing to detect miRNA activity at the moment of incomplete carpel fusion could being morphological distinguished, and identified their targets through degradome analysis, which was previously used to identify miRNA-mRNA target pairs in tomato, maize and grapevine (Pantaleo et al., 2010;Karlova et al., 2013;Liu et al., 2014). Subsequently, the results were validated by quantitative RT-PCR (qRT-PCR) analysis of ovary formation and growth development in IFC and CFC ovaries. In addition, contents of six phytohormones during ovary development, prior to pollination, were determined as well. Data showing the interactions between differentially expressed miRNAs, target genes and phytohormone signals will contribute to our understanding of the molecular foundation of carpel fusion during maize ovary development.

Phenotypic Differences between IFC and CFC ovaries and kernels
Compared with CFC ovaries, IFC ovaries showed an incomplete carpel wall with a hole in the top, near the silk. Moreover, with growth and expansion, the inner tissue (nucellus) extruded from the carpel wall in the IFC ovaries (Figures 1A-D). After pollination, IFC kernels lacking complete pericarp wrapping became easily deformed on contact with other kernels, resulting in an irregular shape compared to CFC kernels at maturity (Figures 1E,F).

Global Analysis of Small RNAs in IFC and CFC Ovaries
To determine the involvement of regulatory miRNAs in carpel fusion during maize ovary formation, we profiled miRNA variation in IFC and CFC ovaries. An average of 11,762,937 and 11,545,117 raw reads were obtained from IFC and CFC ovaries, respectively, which after filtering, represented an average of 6,049,566 valid reads representing 2,658,536 unique sequences and 7,085,972 valid reads representing 3,565,628 unique sequences, respectively (Table 1). These unique sequences were subsequently used to identify known and novel miRNAs by alignment against miRBase (Version 21.0). The length of the sRNAs ranged from 18 to 25 nt, and in both the IFC and CFC libraries, the 24 nt category was most abundant (average of 51.07%; Figure 2A). This is consistent with the typical lengths of plant sRNAs reported in other studies Hackenberg et al., 2015). Of the known miRNAs, 21-nt miRNAs were most abundant (59.7%) (Figure 2B), representing the dominant size of mature miRNAs in plants.
The 5 ′ terminal nucleotides of sRNA sequences influence classification of their AGO complexes and is an important feature affecting function. Most miRNAs are incorporated into the AGO1 effector complex, resulting in sequence specificity that either cleaves or translationally represses their targets (Rogers and Chen, 2013). AGO1 harbors miRNAs that favor the 5 ′ terminal uridine (Mi et al., 2008); therefore, we examined the 5 ′ nucleotide distribution of known miRNAs obtained from IFC and CFC ovaries. Of 172 mature miRNA sequences, 58.15% FIGURE 1 | Phenotypes of IFC and CFC ovaries and kernels at different stages of maize kernel development. (A,B), longitudinal section through a normal pistil and an IFC pistil, respectively. (C,D), CFC and IFC ovaries at the silking stage, respectively. The extruding nucellus is indicated by an arrow in the IFC ovary. (E,F), mature kernels of CFC and IFC ovaries with whole and incomplete pericarps, respectively. Bars = 1 mm. IFC, incompletely fused carpels; CFC, completely fused carpels. started with uridine at the 5 ′ -end, and 25.56% started with Guanine ( Figure 2C). We also examined the sequence frequency and distribution of known and novel miRNA candidates across precursor sequences in the two libraries (Supplementary  Table S1). Similar to other deep sequencing studies (Peng et al., 2011;Tian et al., 2011), unique miRNAs outnumbered unique miRNA precursors, suggesting that two or more mature miRNAs are distributed in different parts of the same precursor body, their expression levels determining the dominant miRNA. Differentially Expressed miRNAs between IFC and CFC Ovaries miRNAs have important functions in plant development and stress responses, and are a well-studied class of regulatory sRNAs.
To determine differential expression between IFC and CFC ovaries, miRNA expression was normalized to transcripts per 1,000,000 and simplified as normalized expression (NE). In a given sample, a miRNA was considered if the NE value was greater than one for a known miRNA. Based on this criterion, a total of 162 (50.47% of the total) mature known miRNAs belonging to 33 families were observed in at least one of the four samples (Supplementary Table S2). In correlation analysis, NE values of the IFC and CFC ovaries were found to be highly correlated between repeats (r = 0.95 and 0.90, respectively; Supplementary Figure S1), indicating good reproducibility of the RNA sequencing results. A one-tailed t-test was used to identify differentially expressed miRNAs with a fold change >1.5 between IFC and CFC as well as a p value < 0.05 and an NE value > 5 in at least one of the samples. As a result, a total of 20 miRNAs were found to be differentially expressed between the two phenotypes ( Table 2). Compared with the CFC ovary, 15 miRNAs were found to be down-regulated and 5 up-regulated in the IFC ovary ( Table 2). Since the IFC and CFC ovaries were taken from the same area on the same ear, most of the detected miRNAs showed similar expression patterns; however, nevertheless, differences did exist.
To determine detailed expression patterns of these miRNAs in the IFC and CFC ovaries, real-time qRT-PCR was performed. Eight differentially expressed known miRNAs were selected for validation. Overall, the results corresponded to the deep sequencing results (Figure 3A), indicating reliability of the miRNA expression levels determined by high-throughput sequencing. In addition, we surveyed expression of these miRNAs in the IFC and CFC ovaries after the initial observation of carpel fusion deficiency, prior to pollination. Changes in expression levels of the selected miRNAs exhibited a consistent tendency with differential expression during IFC and CFC development. It should be noted that miR396 expression increased on the day of silking and thereafter decreased ( Figure 3A).

Target Prediction and Identification of miRNAs by Degradome Sequencing
To identify miRNA targets, two cleaved miRNA target libraries (degradomes) were generated for the IFC and CFC ovaries, respectively. After high-throughput sequencing, 16.9 and 16.7 million raw reads representing the 5 ′ ends of uncapped, polyadenylated RNAs were obtained from the IFC and CFC libraries, respectively. After initial processing, 75.81% (75.88 and 75.73% for the IFC and CFC ovaries, respectively) of the short sequencing reads were mapped to the maize transcriptome, suggesting that some of the filtered reads mapped to unannotated genes (Supplementary Table S3).
In plants, miRNAs have been shown to bind with almost complete complementarity to their mRNA targets, with miRNAmediated cleavage occurring precisely between the 10th and 11th nucleotide from the 5 ′ end of the miRNA in the complementary region of the target mRNA. In the CleaveLand pipeline program (Addo-Quaye et al., 2009), sequenced tags were first mapped to cDNA then the number of tags with a 5 ′ nucleotide corresponding to each position in the mRNA sequence counted, as depicted by a t (target)-plot. That is, cleaved transcripts have distinct t-plot peaks at the predicted cleavage site in the degradome sequence tags. As previously reported (Karlova et al., 2013), the cleavage products can be categorized into three classes: Class I t-plots, category 0 and 1 peaks with a P-value ≤ 0.05; Class II, category 0 and 1 peaks with a P-value > 0.05, and all category 2 peaks regardless of their P-value; and Class III, all category 3 and 4 peaks regardless of their P-value. Thus, Class I includes the most credible miRNA target genes among the three classes. Representative examples of t-plots of cleavage transcripts belonging to the different classes are shown in Figure 4.
As the results, a total of 113 (471 transcripts) target genes were predicted for 73 known miRNAs, 46 (181 transcripts) of which were categorized as Class I (Supplementary Table S4). Of 20 differentially expressed known miRNAs, degradome predication and analysis revealed that only 12 had target genes. Moreover, 20 (39 transcripts) target genes cleaved by nine differentially expressed known miRNAs were classified as Class I from a total of 40 (76 transcripts) identified target genes of all differentially expressed known miRNAs ( Table 3; Supplementary Table S5). The miR_name is composed of the 1st known miR name in a cluster, a underscore, and a matching annotation. Such as: L-/+ n, means the miRNA_seq (detected) is n base less or more than known rep_miRSeq in the left side, respectively; R-/+ n means the miRNA_seq (detected) is n base less or more than known rep_miRSeq in the right side, respectively; 1ss6TC13TA means 1 substitutin (ss), which are T->C at position 6 and T->A at position 13 from 5 ′ of miRNA.
Sixty percent of these target genes were members of different families of transcription factors; namely, TCP, MYB, ARF, NAC, NF-YA, AP2, GRF, GRAS, and HD-ZIP (Table 3), while one gene, ARF, was found to be targeted by miR160 as verified by RNA ligase-mediated 5 ′ -RACE (Supplementary Figure S2). Some of the differentially expressed known miRNAs targeted multiple gene loci; for example miR159 ( Figure 5). Expression of eight miRNA target genes of differentially expressed known and novel miRNAs was subsequently detected using qRT-PCR analysis. As expected, expression of all target genes was negatively correlated with expression levels of the corresponding miRNA in both the IFC and CFC ovaries ( Figures 3A,B, 5), This finding is also in accordance with the gene cleavage function of the miRNAs.

Identification of Potential Novel miRNAs
To identify novel miRNAs, sequencing reads that did not match any of the known miRNAs were further analyzed with reference to the criteria for annotation of plant miRNAs (Meyers et al., 2008). A total of 43 miRNA candidates, present in at least one of the four libraries, were revealed and designated predicated candidates (PC)-5P/3P. In addition, 10 miRNA new isoforms in known miRNA loci were also identified (Supplementary Table  S6). Secondary structures of eight of the selected novel miRNAs are presented in Supplementary Figure S2. Both miRNAs of predicated candidates and new isoforms were defined as novel miRNAs. Of these, a total of 10 were differentially expressed between the IFC and CFC ovaries (Supplementary Table S7).
Potential targets of the novel miRNAs were also predicted via degradome sequencing. A total of 11 target genes (14 transcripts) of 9 novel miRNAs were identified, all belonging to Classes II and III (Supplementary Table S8). However, only seven targets of four differentially expressed novel miRNAs were revealed (Supplementary Table S9). These results clearly differ from the identified targets of known miRNAs (Table 3; Supplementary  Table S9; Figure 4). The four differentially expressed novel miRNAs and their target genes were subsequently selected for qRT-PCR analysis (Figure 6A), revealing a negative correlation between all miRNAs and their corresponding target genes except for PC-5p-189225_15 ( Figure 6B). This exception was possibly due to its low expression level as observed in the sRNA sequencing results (Supplementary  Table S7).

Dynamic Changes in Phytohormones during IFC and CFC Ovary Development
To examine the changes in endogenous phytohormone accumulation during deficient carpel fusion, we compared the contents of various endogenous phytohormones in IFC and CFC ovaries at silking (when the deficiency could be observed in serial observations) and 1, 2. and 3 days after silking, respectively. In general, the endogenous IAA content of the IFC ovaries was significantly higher than that of the CFC ovaries at each developmental stage except the initial observation at silking. In FIGURE 3 | qRT-PCR analysis of the identified differentially expressed known miRNAs and their targets in IFC and CFC ovaries. (A) The copy number of miRNAs was normalized by comparison with maize U6. Relative expression levels of each miRNA were normalized by comparison with expression of the CFC ovaries at developmental stage "0," which was set as 1. (B) The copy number of the corresponding target mRNAs was normalized by comparison with the maize gene EF1a (gene ID: GRMZM2G153541). Relative expression levels of each target were then normalized by comparison with expression of the CFC ovaries at developmental stage of "3," which was set as 1. A primer pair spanning the cleavage site was used to quantify the expression of uncleaved target mRNA. Developmental stage "0" represents silking, when carpel fusion deficiency was initially observed. Developmental stages "1," "2" and "3" represent the number of days after silking. Error bars represent the standard error. Samples used at developmental stage "0" were the same as those used for small RNA sequencing. Lowercase a-d, f, j, k, a-e and g after the miRNA name refer to highly homologous miRNAs. IFC, incompletely fused carpels; CFC, completely fused carpels; *P < 0.05; **P < 0.01.
contrast, the contents of all other endogenous phytohormones, including zeatin riboside and isopentenyl adenine (ZR + iPA), GA, brassinosteroids (BR), jasmonic acid (JA) and abscisic acid (ABA), were lower in the IFC ovaries. The pattern of IAA accumulation was also similar between the IFC and CFC ovaries, both showing a gradual increase; however, the changes in ZR + iPA and GA contents showed different patterns. In the CFC ovaries, ZR + iPA and GA increased with ovary growth, whereas a reverse trend was observed in the IFC ovaries. Changes in endogenous BR, JA and ABA contents showed consistent tendencies, with a decrease in BR and JA and an increase in ABA in both the IFC and CFC ovaries (Figure 7).

Correlation Analysis of Phytohormone Contents, miRNAs and their Targets
Target gene expression was examined using qRT-PCR during ovary formation, prior to pollination. We then focused our attention primarily on miRNAs and their target genes in relation to phytohormone synthesis and metabolism, since they were deemed most likely to play important roles in carpel fusion development. To further investigate the potential relationship between phytohormone signals and the expression of miRNAs and their target genes during IFC ovary development, we selected eight differentially expressed known miRNAs and their targets for correlation analysis. As a result, expression FIGURE 4 | Representative target plots (t-plot) depicting categories of the cleaved mRNAs confirmed by degradome sequencing, reflecting the reliability of the miRNA targets. As previously reported (Karlova et al., 2013), cleavage products can be categorized into three classes: Class I t-plots, category 0 and 1 peaks with a P ≤ 0.05; Class II, category 0 and 1 peaks with a P > 0.05, and category 2 peaks; and Class III, category 3 and 4 peaks. Representative examples of (A,B) target genes belonging to Class I, (C,D) target genes belonging to Class II, and (E,F) target genes belonging to Class III are shown.
patterns of all 8 miRNAs and 8 corresponding targets were found to be significantly correlated with the phytohormone signals (Table 4). Furthermore, these miRNAs and their target genes were subsequently divided into two groups. In the first, the miRNAs and their corresponding targets showed opposite correlation relationships with the phytohormone signals. That is, they all (miR159 and its target with IAA; miR166, miR169, miR393, miR396, and their targets with ZR + iPA; miR159, miR396 and their targets with GA; miR396 and their targets with BR; miR159, miR396 and their targets with JA; miR160, miR164, miR166, miR169, miR172, and miR396 and their targets with ABA) showed a significant positive and/or negative correlation with the phytohormone levels. In contrast, the second group included only those miRNAs or targets that were significantly correlated with the phytohormone signals. For example, miR160 was significantly correlated with IAA, but its target gene was not.
Taking into account both the miRNAs and their targets, it should be noted that either the selected miRNA and / or its target had a significant correlation with at least one phytohormone signal ( Table 4).

DISCUSSION
With the development of sequencing technologies, sRNA sequencing has become an efficient and economical technique for estimating expression profiles of miRNA genes. Recent studies have reported the role of miRNAs in determining gene expression during maize ear and kernel development Xin et al., 2015). However, few studies have examined miRNA expression during ovary growth, in particular, during carpel fusion. In this study, we conducted deep sequencing of sRNAs in IFC and CFC ovaries at the moment of initial  carpel fusion deficiency. As a result, 20 differentially expressed known miRNAs were identified, all of which cleaved target genes related mainly to TFs. Furthermore, a large number of these differentially expressed miRNAs and their targets showed a significant correlation with endogenous hormones during ovary development. These findings suggest that miRNAs, especially those differentially expressed between IFC and CFC ovaries, may play crucial roles in regulating ovary development.

Analysis of the Identified miRNAs and their Target Genes
miRNAs, a class of sRNAs, were previously shown to play an important role in controlling gene expression via cleavage to a target gene(s), many of which are members of TF families (Bartel, 2004;Kidner and Martienssen, 2005). In this study, we obtained an average of 2.7 and 3.6 million unique valid reads for 18-26 nt sequences in the IFC and CFC libraries, respectively, substantially increasing the available data on maize ovary sRNAs. Analysis of these unique sRNAs by mapping to miRBase revealed that most of the identified miRNAs were previously unannotated, suggesting that many more remain to be explored.
In comparison, a total of 30 differentially expressed miRNAs were revealed, of which 20 were known and 10 were novel. These findings suggest that these known miRNAs may play a role in the completion of carpel fusion during ovary development, as further suggested by their annotations in degradome sequencing (Table 3). However, some of these newly identified miRNAs probably had a specific role in carpel fusion or had no action, as implied by the degradome analysis (Supplementary Table S9), this might due to the lower expression of these novel miRNAs and/or the function of these novel miRNAs' targets on ovary development were not fully studied. As mentioned above, some of the identified novel miRNAs showed relatively low expression, consistent with previous observations suggesting that novel miRNAs are often expressed at lower levels than conserved miRNAs (Pantaleo et al., 2010;Liu et al., 2014). miRNAs showing low levels of expression are generally thought to have limited biological function, and therefore, their effects on the target genes would, in most cases, be weak or negligible (Hackenberg et al., 2015). However, despite this, we cannot rule out the possibility that these miRNAs participate in carpel fusion development in maize ovaries.  The copy number of the miRNAs was normalized by comparison with maize U6. Relative expression levels of each miRNA were normalized by comparison with expression of the CFC ovaries at developmental stage of "0," which was set as 1. (B) The copy number of corresponding target mRNAs was normalized by comparison with the maize gene EF1a (gene ID: GRMZM2G153541). Relative expression levels of each target were then normalized by comparison with expression of the CFC ovaries at developmental stage "3," which was set as 1. A primer pair spanning the cleavage site was used to quantify the expression of uncleaved target mRNAs. Developmental stage "0" represents silking, when carpel fusion deficiency was initially observed. Developmental stages "1," "2," and "3" represent the number of days after silking. Error bars represent the standard error. Samples used at developmental stage of "0" were the same as those used for small RNA sequencing. IFC, incompletely fused carpels; CFC, completely fused carpels; *P < 0.05; **P < 0.01. FIGURE 7 | Phytohormone contents at different developmental stages during ovary development in IFC and CFC ovaries. Developmental stage "0" represents silking, when carpel fusion deficiency was initially observed. Developmental stages "1," "2," and "3" represent the number of days after silking. Error bars represent the standard error. Samples used at developmental stage of "0" were the same as those used for small RNA sequencing. IFC, incompletely fused carpels; CFC, completely fused carpels; *P < 0.05; **P < 0.01. Degradome sequencing is a relatively powerful tool for identifying potential targets, and was used here to identify the functions of candidate miRNAs using samples from IFC and CFC ovaries. As the result, a total of 181 transcripts of 46 target genes were identified as Class I miRNAs (Supplementary Table S4). The remaining transcripts all belonged to Classes II and III, the latter of which are considered low confidence miRNA-target pairs (Karlova et al., 2013;Baksa et al., 2015). These findings were possibly due to the complex processes of transcriptional regulation, or the limitations of target identification with degradome sequencing.

TFs Targeted by Differentially-Expressed known miRNAs might be Involved in Regulation of Carpel Fusion
Most gynoecia require complete carpel fusion to ensure ovary formation. miRNAs are important regulators of carpel fusion during gynoecium development (Sieber et al., 2007;Larue et al., 2009). For example, a previous study showed that maize ts4, which encodes miR172, plays a key role in carpel fusion in pistils of tassel flowers (Chuck et al., 2007). Moreover, the TF AP2, one identified target of miR172, was found to be involved in carpel fusion during Arabidopsis gynoecium development (Ripoll et al., 2011). It was also revealed that miR156 controls the initial steps of fleshy fruit development in tomato, playing an important role in ovary and fruit development (Ferreira et al., 2014). Moreover, two targets of miR156 were found to positively regulate miR172 expression by binding their sequence to the regulatory region of miR172 (Wu et al., 2009). In our study, expression of miR172 was down regulated, while that of miR156 was up regulated in IFC compared to CFC ovaries ( Table 2), suggesting a similar interaction between miR172 and miR156 may exist in maize ovary as well. Overall, the degradome and qRT-PCR analyses suggest that the differentially expressed miR172 might also regulate carpel fusion during maize ovary development by targeting AP2.
Carpel fusion requires coordination among various functional genes. A large proportion of the genes known to be involved in normal carpel fusion are TF genes, two of which, MYB and TCP, are known to play a key role in regulation of carpel fusion development in Arabidopsis (Shin et al., 2002;Koyama et al., 2007), were found to be targeted by miR159 (Table 3). Moreover, mRNA expression of both TFs was found to show an opposite trend compared to miR159 during ovary development in both IFC and CFC ovaries (Figures 3, 5), suggesting a consistent regulatory role of miR159 during ovary development in maize. miR164 and its target gene, NAC, have been described in relation to carpel marginal tissue development in Arabidopsis (Sieber et al., 2007). Furthermore, CUC2, a member of the NAC family, is known to be a key regulator of carpel fusion development in Arabidopsis (Nahar et al., 2012). In addition, miR169 and its targets, members of the NF-YA family of TFs, were previously found to exert homeotic control over the carpel identity gene AG in Petunia hybrida and Antirrhinum majus (Cartolano et al., 2007). Overexpression of miR169 was also found to cause changes in fruit shape and size in tomato (Teotia et al., 2015). In this study, AP2, MYB and TCP, and other TF family members including NAC, ARF, GRF, and NF-YA, were identified as Class I target genes of the differentially expressed known miRNAs (Table 3). Taken together, these findings suggest that these miRNAs may play a vital role in controlling carpel fusion development.

Differentially Expressed miRNAs may Control Ovary Development Through Regulation of Phytohormone Homeostasis
During formation and growth of the maize pistil, a period of 3 days after silking is sufficient for pollination of the maximum number of whole ear kernels (Culy et al., 1992). Thus, at this developmental stage, the pistil undergoes appropriate ovary development, with sufficient silk vitality and silk length, and develops sufficiently for fertilization (Culy et al., 1992). Therefore, understanding phytohormone levels at this stage as well as expression levels of miRNAs and their targets is important, as is determining their interaction and regulation in IFC and CFC ovary development.
In cereals, carpel fusion and gynoecium development is reportedly regulated by plant hormones (Barazesh and Mcsteen, 2008;Larsson et al., 2013;Reyes-Olalde et al., 2013;Marsch-Martínez and de Folter, 2016). Accumulating evidence suggests that miRNAs affect plant hormone signaling and hormone synthesis, In addition, sensing and transport of these hormones relies on the activity of proteins, with most genes involved in this regulatory network regulated by miRNAs (Curaba et al., 2014;Marsch-Martínez and de Folter, 2016).
Auxin, which stimulates cell differentiation, is extremely important to plant development, particularly gynoecium development (Hawkins and Liu, 2014;Sehra and Franks, 2015). Experimental evidence suggests that auxin also plays a crucial role in carpel marginal tissue development in Arabidopsis (Reyes-Olalde et al., 2013). In rice, ARF6 controls inflorescence development through coordinated activation of auxin biosynthesis and auxin response factors (Gao et al., 2015). ARF8 and ARF6, important signaling components in the auxin signaling pathway, are closely related TFs both targeted by miR167, and involved in carpel maturation and development in Arabidopsis (Nagpal et al., 2005). Our results showed no difference in IAA content when carpel fusion deficiency was first observed; however, soon after, the content became significantly higher in IFC than CFC ovaries (Figure 7). This is consistent with previous studies suggesting that a high concentration of auxin is essential for carpel fusion (Larsson et al., 2013;Sehra and Franks, 2015). Gene expression of auxin response factor (ARF), a target gene of miR160, was also higher in IFC ovaries (Table 2; Figure 3), indicating a major role in auxin homeostasis and carpel fusion. This is consistent with a previous study showing that mutation in the ARF gene can lead to alterations in carpel and gynoecium development, and patterning in Arabidopsis (Nagpal et al., 2005;Crawford and Yanofsky, 2011).
Cytokinins (CKs) are a structurally diverse species of N 6substituted purine derivatives that stimulate mitotic divisions, and also known to participate in carpel marginal tissue development (Reyes-Olalde et al., 2013;Marsch-Martínez and de Folter, 2016). OsGRF4, a target of miR396, reportedly regulates a CK dehydrogenase precursor gene and controls rice grain shape (Sun et al., 2016), while TCP14 and TCP15 TFs, which mediate CK responses, were found to cause excessive proliferation of the boundaries of the Arabidopsis gynoecium replum (Takei et al., 2001). Isopentenyl adenine (iPA) and zeatin riboside (ZR) are the major CK species in maize (Takei et al., 2001). Our data showed lower iPA and ZR contents in IFC ovaries except at the initial stage of incomplete carpel fusion (Figure 7), which is possibly related to the gradual increase in morphological differences between the IFC and CFC ovaries. Our data also revealed that expression of TCP, the target gene of miR159, was higher in IFC than CFC ovaries (Table 2; Figure 3), suggesting a regulatory function of miR159 in carpel fusion development.
Gibberellin (GA) is thought to play diverse roles in plant growth and development, including flowering time, with overexpression of miR156 reducing GA responses during flowering (Yu et al., 2012). The GA3ox1 gibberellin biosynthesis gene is a direct target of INDEHISCENT (IND) TF; however, ind mutant plants, which have low GA levels, show abnormal carpel valve margin development in the Arabidopsis gynoecium (Arnaud et al., 2010). The role of brassinosteroids (BR) in gynoecia development was also recently examined. For example, OsGRF4, a target of miR396, was found to positively regulate BR content through direct interaction with GSK2, the central negative regulator of brassinosteroid signaling in rice (Che et al., 2015). The enzyme CYP85A2 is also known to participate in brassinolide biosynthesis, with the double mutant seu cyp85A2 causing carpel fusion defects in the gynoecial apex of Arabidopsis (Nole-Wilson et al., 2010). Furthermore, decreased production of the plant hormone JA reportedly accounts for a subset of arf6 arf8 mutant phenotypes, including developmental defects in carpels, (Nagpal et al., 2005). Moreover, in Arabidopsis, it is reported that ARF6 and ARF8 are directly regulated by miR167 at the post-transcription level (Wu et al., 2006). miR159 also targets genes encoding the MYB TF, which is involved in regulation of ABA (Reyes and Chua, 2007). ABA was also shown to have a positive effect on flowering initiation and a negative effect on flower development (Wang et al., 2002). Therefore, these miRNAs and their targets as well as phytohormones all together may contribute to a balanced manner during ovary formation and development, with interruption of any one of them might result in an abnormal ovary development.
In conclusion, in this study, we performed genome-wide identification of known and novel miRNAs and experimentally predicated their targets through degradome analysis. Furthermore, the changes in endogenous phytohormone contents between IFC and CFC ovaries were also examined, and the correlation with miRNAs and their targets determined. The presence of differentially expressed miRNAs during transcription regulation and phytohormone signaling processes in IFC and CFC ovaries suggests the significant roles of miRNAs in carpel fusion development. This study provides valuable information on carpel fusion development that will be beneficial in breeding programs aimed at enhanced seed vigor and quality in maize. Further genetic studies are now needed to determine the nature of these regulatory interactions.

Plant Materials
Zea mays (maize) inbred line Yu-A474 was planted in the farm of Henan Agricultural University, Zhengzhou city (Henan Province, China), under non-stress conditions during June-October 2014. This inbred line is a female of the maize hybrid Yudan 603, some ovaries of which show various phenotypes of incompletely fused carpels. Before silk emergence, ears were bagged to prevent pollination. To determine the timing of initial emergence of the IFC phenotype, we carried out continuous and repeated microscopic observations from the stage of floret primordium differentiation to floret organ differentiation during ear development in August 2013 (Supplementary Figure S3). Ovary samples were collected at silking, the time at which carpel wall fusion deficiency was initially observed. Briefly, ovaries were manually collected at the base using forceps after removal of the glumes, lemma and palea. To eliminate inconsistencies from sampling different parts of the ear, IFC ovaries were sampled close to the CFC ovaries, in the middle part of the same ear. All collected ovaries were immediately frozen in liquid nitrogen and stored at −80 • C for RNA isolation. Two replicates from two different ear rows were collected, respectively, for the two ovary phenotypes. In total, four samples, including two biological replicates from CFC and IFC, respectively, were collected and prepared for total RNA extraction.
To confirm expression patterns of the 12 differentially expressed miRNAs and their target genes, IFC and CFC ovaries were collected on the day of silking then 1, 2, and 3 days after silking, prior to pollination. This ensured that the period of ovary formation and growth was covered as well as meeting the conditions of fertilization.

RNA Isolation, Small RNA Library Construction and Sequencing
Total RNA for sRNA sequencing, qRT-PCR and degradome analyses was extracted from IFC and CFC ovaries using Trizol reagent (Invitrogen) according to the manufacturer's instructions. Quality and purity of the total RNA was analyzed using Bioanalyzer 2100 (Agilent) and an RNA 6000 Nano LabChip Kit (Agilent). From each sample, 1 µg of total RNA was ligated to RNA-DNA chimeric oligonucleotide adaptors then converted to cDNA by RT-PCR. The resulting cDNA was amplified by PCR and gel-purified to produce sequencing libraries. Finally, sRNA sequencing (single-end, 50 bp) was performed on an Illumina Hiseq2500 (Illumina) platform according to the manufacturer's recommended protocol.

Identification of known miRNAs
To identify known miRNAs, the following were performed. Raw reads were subjected to the Illumina pipeline filter (Solexa 0.3) then the dataset processed using ACGT101-v4.2-miR (LC Sciences, Houston, Texas, USA) to remove adapter sequences, junk reads, reads less than 15 nt, common RNA families (rRNA, tRNA, snRNA, and snoRNA), repeats and sRNA reads assigned to exon regions. The remaining small RNAs were classified by alignment to mRNA, RFam and Repbase databases then filtered.
The last remaining unique sequences, which were 18-25 nt in length, were mapped to maize precursors in miRBase 21.0 (ftp:// mirbase.org/pub/mirbase/CURRENT/) using a BLAST search to identify known miRNAs and new isoforms. Length variation at the 3 ′ and 5 ′ ends and one mismatch inside the sequence were allowed in the alignment. Unique sequences mapped to mature maize miRNAs in hairpin arms were classified as known miRNAs. Unique sequences mapped to the other arm of known maize precursor hairpins, opposite the annotated mature miRNA-containing arm, or to different positions on the same arm were considered new isoforms. Remaining sequences were mapped using bowtie, an alignment tool in the proprietary pipeline script Houston,Texas,USA), to precursors of other selected species (with the exclusion of maize) to identify known and new isoforms (one mismatch inside these sequences was allowed with a seed length of 16 nt). The identified precursors were further mapped to the maize genome using Megablast, another alignment tool in ACGT101-v4.2-miR (LC Sciences, Houston, Texas, USA), to determine their genomic locations. Only those showing a similarity rate of more than 90% were selected.

Identification of Potential Novel miRNAs
After identification of known and new isoforms, remaining sequencing reads that did not match any known miRNA precursors were subjected to "ACGT101-v4.2-miR" to further determine novel miRNAs. Criteria were mainly those of Meyers and Lee (Meyers et al., 2008;Peng et al., 2013). Parameters for detailed identification of secondary structures were as follows: (1) number of nucleotides in one bulge in the stem region <13; (2) number of base pairs in the stem region of the predicted hairpin >15; (3) cutoff of free energy during the formation of secondary structures < −15 kcal/mol; (4) length of the hairpin (up and down stems + terminal loop) >49; (5) length of the hairpin loop <351 nt; (6) number of nucleotides in one bulge in a mature region <5; (7) number of biased errors in one bulge in a mature region <3; (8) number of biased bulges in a mature region <3; (9) number of mismatch errors in a mature region <5; (10) number of base pairs in a mature region of the predicted hairpin >11; and (11) percentage of mature base numbers in the stem >80%. Furthermore, unmapped sequences were BLASTed against the maize genome, and the mapped sequences containing hairpin RNA structures predicated from flank 120-nt sequences using mfold software (part of ACGT101-v4.2-miR).

Analysis of Differentially Expressed miRNAs
To compare differentially expressed miRNAs between IFC and CFC, expression abundances were normalized to obtain the expression of transcripts per 1,000,000 using the following formulae: Normalized expression (NE) = (Actual miRNA reads count/Total count of clean reads) × 1,000,000.
Before identification of differentially expressed miRNAs, we conducted reproducibility analysis of data from two replicate IFC and CFC ovaries using SCC analysis (Zhan et al., 2015). Log 2transformed NE values [log 2 (NE + 1)] of the expressed miRNAs were used as input for the SCC analysis. Differential expression was analyzed using a t-test, based on normalized counts with a significance threshold set as a fold change of NE >1.5 and a Pvalue < 0.05 (Wang et al., 2012;Boke et al., 2014;Chang et al., 2015).

Degradome Library Construction and Bioinformatics Analysis
Two degradome cDNA libraries were constructed from the same ovary samples used for sRNA analysis, mixing the two biological samples into one for IFC and CFC, respectively. Single-end sequencing (50 bp) was performed on an Illumina Hiseq2500 (Illumina) according to the method of German et al. (2008) with some modifications. Briefly, the extracted poly (A) RNA was ligated to a 5 ′ adapter containing a MmeI site at its 3 ′ end. The ligated products were then used for cDNA production and amplified by PCR for 5 cycles. The PCR products were purified and digested with MmeI, and the resulting fragments ligated to a second double-stranded DNA oligonucleotide. The ligation products were further purified and amplified for another 10 PCR cycles, and the final product purified and subjected to high throughput sequencing.
The publicly available CleaveLand pipeline version 3.0.1 software package (Addo-Quaye et al., 2009) and Target Finder program (http://targetfinder.org/) (Kiełbasa et al., 2010) were used to detect potentially sliced targets of the known and novel miRNAs identified in sRNA sequencing. To account for inaccurate target cleavage or variations in miRNA 5 ′ ends, the pipeline was modified to recognize targets cleaved at the 9th, 11th, and 10th positons. All targets were classified as t-plot peaks according to 5 categories (0-4) based on the abundance of the resulting mRNA tags relative to the overall profile of the degradome reads matching the target (Addo-Quaye et al., 2008). Classification was as follows: peaks in categories 0-3, >1 read per peak; category 0, peaks representing a single maximum in a particular transcript; category 1, peaks equal to the maximum, with more than one maximum per transcript; category 2, peaks lower than the maximum but higher than the median of a transcript; and category 3, peaks with an equal or less than median number of reads. Category 4 peaks had only 1 read. The statistical significance of an observed peak-miRNA match was represented by a P-value < 0.05.

Quantitative Real-Time RT-PCR (qRT-PCR)
Expression levels of the differentially expressed miRNAs and the predicated targets were validated by qRT-PCR analysis. Eight known miRNAs (miR159, miR160, miR164, miR166, miR169, miR172, miR393, and miR396), four newly identified miRNAs, and 12 target genes were selected for qRT-PCR validation. qRT-PCR was performed using the SYBR Green PCR Master Mix (Applied Biosystems) and an ABI PRISM_7900 Sequence Detection System (ABI, USA) following the manufacturer's instructions. To detect the level of non-cleaved mRNAs, all primers designed for the target genes spanned the miRNA cleavage site. Primer sequences of the miRNAs and their targets are listed in Supplementary Table S10. U6 RNA and EF1a (gene ID: GRMZM2G153541) were used as an internal reference for miRNA and their target genes, respectively. All reactions were performed in triplicate for technical and biological repetition of the IFC and CFC ovaries, respectively, and the generated real-time data analyzed using the comparative 2 − Ct method.

RNA Ligase-Mediated 5 ′ Race
To determine the cleavage sites of the target transcripts, we performed RNA ligase-mediated rapid-amplification of 5 ′ complementary DNA ends (5 ′ -RLM-RACE) with a GeneRacer kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions with slight modifications. Briefly, total RNA from mixtures of IFC and CFC ovaries were ligated directly to a 5 ′ RACE RNA adapter (5 ′ -GCTGATGGCGATGA ATGAACACTGCGTTTGCTGGCTTTGATGAAA-3 ′ ) followed by reverse transcription with the GeneRacer (dT) primer. The reverse transcription product was used as template for PCR, with GeneRacer TM 5 ′ Primer (5 ′ -GCTGATGGCGATGAATGAACA CTG-3 ′ ), GeneRacer TM 5 ′ Nested Primer (5 ′ -CGCGGATCCG AACACTGCGTTTGCTGGCTTTGATG-3 ′ ), and two genespecific reverse primers used in each RACE (Supplementary Table S10). The RACE products were gel-purified, cloned to the pGEM-T Easy vector (Promega, Madison, WI, USA), and 10 independent clones from each reaction sequenced.

Detection of Phytohormones
Samples used for phytohormone measurement were the same as those used for qRT-PCR validation. Approximately 0.2 g of frozen samples were ground and homogenized in 2 ml of 80% methanol extraction medium (containing 40 mg l −l butylated hydroxytoluene as an antioxidant). The extract was incubated at 4 • C for 24 h then centrifuged at 4,000 revolutions per minute for 15 min at 4 • C. The supernatant was passed through Chromosep C18 Sep-Pak cartridges (Waters Corp., Millford, MA, USA) then prewashed with 10 ml of 100% (v/v) methanol followed by 10 ml 80% (v/v) methanol. The eluate was dried under pure N 2 and the extracts dissolved in 2.0 mL phosphate-buffer saline (PBS) (pH 7.5) containing 0.1% (v/v) Tween 20 and 0.1% (w/v) gelatin to examine free IAA, ZR, iPA, GA, BR, ABA, and JA via indirect competitive enzyme-linked immunosorbent assay (icELISA) according to Yang et al. (2001), Zhao et al. (2006), and Deng et al. (2008). Mouse monoclonal antigen and antibodies against free IAA, ZR, iPA, GA, BR, ABA, and JA were produced according to Weiler et al. (1981) at the Center of Crop Chemical Control, China Agricultural University, China. JA was derivatized into methyl jasmonate using JA in its free-acid form, and used for ELISA analysis. Anti-ZR antibody was used to detect ZR-type CKs, and the anti-iPA antibody was used to detect iP and iPR (iPA-type CKs). Calculation of the ELISA data was performed as described in Weiler et al. (1981). Percentage recoveries obtained using internal standards during extraction and analysis were all above 90%.

Statistical Analysis
A one-tailed t-test was used to compare the significance of differences in gene expression data from RT-PCR analysis.
One-way ANOVA was used to compare the differences in phytohormone contents at different stages before pollination. All values are reported as means ± standard error (SE). Differences were considered significant at P < 0.05. To evaluate the relationship between phytohormone contents and expression levels of miRNAs and their targets, the qRT-PCR results and phytohormone levels of both the IFC and CFC ovaries were used for correlation analysis. All data analysis was performed using SPSS 19.0 (SPSS Inc., Chicago, IL, USA).

Data Access
The raw data of sRNA sequencing and degradome supporting the results of this study are available in the National Center for Biotechnology Informations (NCBI) Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), under accession GSE80998. The following link was created to allow review of record GSE80998: http://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?token=khehakekjfmfjsd&acc=GSE80998

AUTHOR CONTRIBUTIONS
CL participated in design of the study and wrote the manuscript. HL performed the experiments and data analysis, and wrote the draft manuscript. TP contributed to the data analysis and revised the manuscript. YW and QW participated in data analysis. JC and MZ participated in data analysis and field management. GT edited the manuscript and gave valuable advice on chart display. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
This work was supported by the China Agriculture Research System (No. CARS-02-19) and the Special Fund for Agroscientific Research in the Public Interest (No. 201203100). We are grateful to Dr. Dong Ding for his help with the initial data analysis.