Identification of soybean seed developmental stage-specific and tissue-specific miRNA targets by degradome sequencing

Background MicroRNAs (miRNAs) regulate the expression of target genes by mediating gene silencing in both plants and animals. The miRNA targets have been extensively investigated in Arabidopsis and rice using computational prediction, experimental validation by overexpression in transgenic plants, and by degradome or PARE (parallel analysis of RNA ends) sequencing. However, miRNA targets mostly remain unknown in soybean (Glycine max). More specifically miRNA mediated gene regulation at different seed developmental stages in soybean is largely unexplored. In order to dissect miRNA guided gene regulation in soybean developing seeds, we performed a transcriptome-wide experimental method using degradome sequencing to directly detect cleaved miRNA targets. Results In this study, degradome libraries were separately prepared from immature soybean cotyledons representing three stages of development and from seed coats of two stages. Sequencing and analysis of 10 to 40 million reads from each library resulted in identification of 183 different targets for 53 known soybean miRNAs. Among these, some were found only in the cotyledons representing cleavage by 25 miRNAs and others were found only in the seed coats reflecting cleavage by 12 miRNAs. A large number of targets for 16 miRNAs families were identified in both tissues irrespective of the stage. Interestingly, we identified more miRNA targets in the desiccating cotyledons of late seed maturation than in immature seed. We validated four different auxin response factor genes as targets for gma-miR160 via RNA ligase mediated 5’ rapid amplification of cDNA ends (RLM-5’RACE). Gene Ontology (GO) analysis indicated the involvement of miRNA target genes in various cellular processes during seed development. Conclusions The miRNA targets in both the cotyledons and seed coats of several stages of soybean seed development have been elucidated by experimental evidence from comprehensive, high throughput sequencing of the enriched fragments resulting from miRNA-guided cleavage of messenger RNAs. Nearly 50% of the miRNA targets were transcription factors in pathways that are likely important in setting or maintaining the developmental program leading to high quality soybean seeds that are one of the dominant sources of protein and oil in world markets.


Background
MicroRNAs (miRNAs) are endogenous noncoding small RNAs which play significant roles in the regulation of gene expression. Post-transcriptional gene regulation by miRNAs constitutes one of the most conserved and well characterized gene regulatory mechanisms. It is important for growth, development, stress responses and numerous other biological processes in eukaryotes [1][2][3][4][5]. Therefore, identification of miRNAs and their targets in diverse species has been a major focus in recent years [6,7]. In higher plants, miRNAs play significant roles in different developmental stages by regulating gene expression at transcriptional and post-transcriptional levels [8][9][10][11][12].
Most plant miRNAs facilitate the degradation of their mRNA targets by slicing precisely between the tenth and eleventh nucleotides (nt) from the 5' end of the miRNA. As a result, the 3' fragment of the target mRNA possesses a monophosphate at its 5' end. This important property has been used to validate miRNA targets [8]. Isolation of such fragments is one of the critical steps for validating miRNA guided cleavage of target mRNA. A major limitation of this procedure is that every single predicted gene has to be verified separately. So, one-ata-time isolation of cleaved RNA fragments is laborious, time-consuming and expensive. Recently, highthroughput sequencing methods, known as degradome analysis or PARE (parallel analysis of RNA ends) that can globally identify small RNA targets have been developed to overcome such limitations [13,14].
Soybean is one of the most important crops cultivated all over the world. It is a good source of vegetable protein and oil. However, the role of miRNAs in soybean seed development is mostly unknown. So it is important to identify the seed developmental stage-specific and tissues-specific miRNAs and their potential target genes. Identification of the consequences of miRNA-guided target degradation that occurs in a developmental and tissue specific manner could help to elucidate how lipid and protein metabolic pathways operated during seed development. The soybean genome (cv. Williams82) was decoded a year ago [15], and this information has accelerated molecular research on soybeans. Although many soybean miRNAs were identified in previous research [16][17][18][19][20], the number of miRNAs known in soybean is still very small and considerably lower than that in Arabidopsis or rice. High-throughput sequencing technologies such as massively parallel signature sequencing (MPSS), 454 and sequencing-by-synthesis (SBS) have enabled the identification of miRNAs in soybean. The extent of miRNA-directed post-transcriptional gene regulation in any organism can only be fully realized by identifying not only the miRNA component but also the set of their RNA targets.
Recently, miRNA targets have been reported for one of the many stages of soybean seed development, namely very early at 15 days after flowering, and without dissection of the maternal seed coats from the cotyledons which develop from the zygote [21]. To comprehensively investigate small RNA targets and provide basic information for further understanding of the miRNAmediated post-transcriptional regulation during different soybean seed developmental stages, we constructed five separate degradome libraries derived from seed coats and cotyledons of different developmental stages representing the early, mid, and late maturation stages of seed development. The libraries were sequenced using SBS sequencing technology. The degradome dataset for the five different libraries was computationally analyzed. The majority of these reads mapped to the soybean transcriptome. A total of 183 target genes were confirmed as miRNA targets, which included both conserved and non-conserved miRNAs. Additionally, we have identified targets for 25 cotyledon-specific miRNAs, as well as 12 miRNAs and their potential targets found only in the seed coats. We found 16 miRNA families and their large number of targets that are found in both tissues. Moreover, we have validated Auxin Response Factors (ARFs) to be targets of gma-miR160, as verified by RNA ligasemediated 5' rapid amplification of cDNA ends (RLM 5'-RACE). The identification of developmental stagespecific and tissue-specific miRNA targets including many transcription factors advance our understanding of the miRNA-mediated post-transcriptional gene regulation during soybean seed development.

Results
Seed developmental stage-specific library construction, sequencing and sequence analysis In higher plants, most miRNAs regulate their targets via cleavage, which normally occurs between the tenth and eleventh nucleotides of the complementary region between the miRNA and the mRNA target [22]. The 3' cleavage fragments contain both a free 5' monophosphate and a 3' polyA tail. So, these cleavage products can be successfully ligated with RNA ligase, whereas full length cDNAs with a 5' cap structure or other RNAs lacking the 5' monophosphate group are not compatible for ligation [8] and thus will be unavailable for subsequent amplification and sequencing reactions.
Five different degradome libraries, which capture the cleaved mRNAs, were constructed from cotyledons and seed coats from different seed developmental stages. These represented early maturation seed (25-50 mg fresh weight, green seed) and mid-maturation (100-200 mg fresh weight, green seed), the stages when the biosynthetic capacity of the seed is maximal and proteins and oils are accumulated at a high rate. In addition, we constructed a library from the yellow cotyledons (300-400 mg fresh weight) that are undergoing dehydration and desiccation.
SBS sequencing of these libraries produced raw reads from 10 million to 45 million (Table 1). After removal of low quality sequences and adapter removal, 95% of these reads lengths were 20 or 21 nt in length as expected from the cloning procedure. More than 97% of reads mapped to the soybean genome available at the Phytozome data base [23]. We also used the computationally predicted cDNA transcripts from the soybean genome sequence consisting of 78,773 high and low confidence gene models (Glyma models) for mapping degradome reads and found that more than 95% of reads matched to the Glyma models. These data indicate our degradome libraries to be of high quality and efficiency in recovering degraded mRNA targets that should contain the sequence profile resulting from miRNA directed cleavage.

Systematic identification of miRNA targets in soybean
Systematic identification of miRNA targets was accomplished using previously described methods by analyzing the 20 and 21 nt reads with the CleaveLand pipeline for miRNA target identification [24] using all Glycine max miRNAs from miRBase [25,26]. The identified targets were grouped into five categories by the program based on the relative abundance of the number of reads mapping to the predicted miRNA target site relative to other sites in the gene model (see Methods for details). Those in category 0 clearly have the majority of tags located at the miRNA-guided cleavage site.
The identified miRNA targets using degradome sequencing are presented in the form of target plots (t-plots) that plot the abundance of the signatures relative to their position in the transcript [14]. Representative t-plots are shown, one from each of four different degradome libraries such as cotyledon (25-50 mg), seed coat (25-50 mg), cotyledon (100-200 mg) and seed coat (100-200 mg) ( Figure 1). In each of the four Glyma models, a clear peak for the absolute number of tags is found at the predicted cleavage site for gma-miR159, gma-miR4406, gma-miR167, or gma-miR164.
The location of the cleavage site within the target gene is another important aspect of miRNA-mediated gene silencing. In soybean, the cleavage site of the miRNA was usually located in the CDS (coding sequence) of the target genes (Tables 2, 3 and 4 and Additional file 1). Since the soybean genome at Phytozome [23] used computational predictions of gene models (known as Glyma models), some are likely deficient at the 5' and 3' UTRs (untranslated regions). Due to the some gene models being incomplete in the UTRs, there are likely other genes targeted by miRNA-guided cleavage in the UTR regions that may not be detected in our alignment analyses. In addition, miRNAs that function through translational repression, as opposed to cleavage of the target mRNA, will also not be identified by degradome or PARE sequencing techniques.
The full complement of targets found in each of the five degradome libraries is presented in Additional file 1. In total, 183 targets representing 53 different miRNAs families were identified. Of those 133 targets were found representing the putative action of 16 different miRNAs in common between both tissues. Table 2 presents a subset of those that are found in at least one stage of development for both seed coats and cotyledons. The Clea-veLand program predicts any gene family members that have a splice site matching the degradome data. Some miRNA family members residing at different genomic locations have very similar, if not identical mature miRNA sequences. Thus, the predictions from analysis of degradome data do not necessarily mean that the particular miRNA family member revealed from degradome data is the one expressed in that tissue. Direct sequencing of the small RNA population is required to verify the presence of a particular gene family member. Inspection of small RNA sequencing data from seed coats and cotyledons of Williams [27] shows the presence of various miRNA family members for gma-miR156, 159, 160, 164, 166, and 167, thus confirming that these miRNAs are present during seed development.  Identification of miRNA targets specific to either seed coat or cotyledons during seed development Tissue specific miRNA and target identification is very important for understanding the regulation of gene expression in a spatial manner. In this study, we constructed cotyledon and seed coat libraries separately to identify miRNA targets both at younger (25-50 mg fresh weight seed) and older stages (100-200 mg and 300-400 mg fresh weight seed) of soybean seed development. Tissue specific siRNAs generated from a cluster of inverted repeat chalcone synthase (CHS) genes that downregulate CHS mRNAs and lead to lack of pigment Cot25 and SC25 are libraries from cotyledons and seed coats dissected from early-maturation green seed of 25-50 mg fresh wt.; Cot100 and SC100 are cotyledons and seed coats dissected from mid-maturation green seed of 100-200 mg fresh weight range; Cot300 are cotyledons dissected from late maturation yellow seeds of 300-400 mg weight range. * and bold are miRNAs whose targets shown in the target plots (t-plots) (Figure 1). # and bold are miRNAs whose RLM-5'RACE validated targets are shown in the target plots (t-plots) ( Figure 2).
on soybean seed coats have been described [27,28], but very little is known about the miRNAs and their targets in developing seed tissues. We analyzed the degradome data from seed coats versus cotyledons and identified 25 miRNAs and their 32 different targets that were found only in the cotyledons and not the seed coats (Table 3). Likewise, 12 miRNAs and their 18 targets are associated with the seed coats only ( Table 4).

Validation of miRNA targets
We report here that many targets were captured by the degradome analysis, which provided experimental evidence to support previous computational predictions. Because of its polyploid genome, many soybean genes are present in multiple copies. As a result, some of the reads align to multiple members of the same gene family. To further confirm the degradome data for some of the family members, a RLM-5' RACE experiment was performed to examine which family members were targeted by the miRNA for degradation. For gma-miR160 in the cotyledon (100-200 mg) degradome library (Table 2 and Additional file 1), we have identified five targets annotated as Auxin Response Factors (ARFs). Four of the five, namely Glyma12g08110.1, Glyma12g29720.1, Glyma14g33730.1 and Glyma11g20490.1, were also verified by RLM-5'RACE to be subjected to cleavage guided by gma-miR160 ( Figure 2).

GO analysis of miRNA target genes in soybean seed developmental stages
The identified targets for miRNAs in the three cotyledon degradome libraries were classified by their gene ontology (GO) using the AgriGO toolkit [29] (Figure 3). Higher percentages of these targets were found to be involved in developmental, reproductive, and regulatory and metabolic processes with respect to their proportions within the GO classification of all soybean cDNAs. The same general pattern is found for the targets predicted with the seed coats (Additional file 2). The enrichment of the genes involved in developmental and regulatory processes may be consistent with the fact that the degradome libraries were constructed from different stages of developing soybean seeds. For the developing seeds, it is of utmost important to accumulate proteins and lipids that are subsequently used as the source of energy and amino acids for the germinating seedling. The corresponding miRNAs may regulate the expression of these target genes during different seed developmental stages in soybean through affecting various transcription factors that induce or shut off specific metabolic networks during the course of seed development.
Interestingly, we identified more miRNA targets in the cotyledons of late seed maturation than earlier stages with a total of 92 different targets in the 300-400 desiccating, yellow seeds compared to 60 and 53 total in the early and mid-maturation, immature green seed respectively.

Discussion
Regulation of gene expression by miRNAs has been comprehensively investigated in animals and plants [3,4,22,30,31]. In the case of higher plants, Arabidopsis and rice miRNA targets have been widely studied by high throughput sequencing [13,[32][33][34]. Soybean is a polyploid crop plant having a complex and large genome compared to Arabidopsis and rice. The number of identified miRNAs and their potential targets in soybean is limited. To date, degradome sequencing has been reported for only one soybean tissue, namely the very young whole seed extracted 15 days after flowering from the cultivar Heinong44 [21]. In order to study the regulation of gene expression during soybean seed development, we constructed and sequenced five distinct degradome libraries using cotyledons and seed coats as a tissue source from different seed developmental stages of the cultivar Williams. The seed coat is primarily maternal tissue while the cotyledons represent the embryo * and bold is miRNA whose target is shown in the target plot (t-plot) (Figure 1). of the next generation. Both tissues have the same genotype in inbred soybean lines. After computational analysis using the Cleaveland pipeline [24], we identified a total of 183 potential targets of 53 miRNA families in five different seed developmental stage specific degradome libraries. Subsequent analysis and identification of cotyledon and seed coat specific miRNAs and their targets give us a better understanding of the regulation of gene expression in a spatial manner during soybean seed development at later stages of seed development in a widely used Maturity Group III cultivar "Williams". The soybean genome sequence and predicted gene models used in this study are derived from a closely related isoline Williams82 [15]. Validation of four Auxin Response Factor (ARF) genes as targets of gma-miR160 indicates degradome sequencing as an efficient strategy to identify miRNA targets in plants.
Comparison of the data in our report with the first soybean degradome data reported by Song et al. [21] from a very early stage of soybean seed development in the cultivar Heinong44 revealed a number of miRNA targets in common among the data sets including many transcription factors such as ARF, MYB, TCP, NF-Y, Growth Regulatory Factor, HD-ZIP, PPR, SBP and NAC family protein. In contrast, we found some other miRNA targets such as Permease family protein, LRR domain containing protein, transmembrane protein 14 C, Serine-Threonine kinase, BRE expressed protein and transcription factor TFIID only in our degradome libraries that represented later stages of seed development of the cultivar Williams. Overall there were 65 Glyma model targets in common between both of the data sets, 80 found only in the Heinong44 data set, and 118 found only in the Williams data sets.
The 183 identified targets from our degradome analyses belonged to 57 different annotation groupings (Additional file 1). Comparing those to the Song et al. paper, 11 annotation groups are in common, 18 are unique to the Heinong44 very young seed, and 46 are unique to the five libraries of the Williams data sets ( Table 5). Many of the identified soybean miRNA targets belong to diverse gene families of transcription factors such as ARFs, MYBs, TCPs, NACs, HD-ZIPs and NF-Y subunits (Tables 2, 3 and 4 and Additional file 1). Many of these transcription factors are known to regulate diverse aspects of plant growth and development. Patterning and outgrowth of lateral organs in plants depend on the expression of HD-ZIP transcription factors that specify adaxial/upper cell fate [35,36]. We identified a number of HD-ZIP transcription factors as targets for gma-miR166 in our degradome libraries. MYB family members in rice, which are targeted by miR159, appear to play an important role in response to the presence  of abscisic acid (ABA) during plant embryonic development, suggesting their roles in seed development [37,38]. In this study, we detected a number of MYB family transcription factors regulated by gma-miR159 tlsb 01w(Additional file 1. The gma-miRNA156 family members target sites in numerous proteins containing the Squamosa Promoter Binding (SBP) domain. SBP and SBP-LIKE (SPLs) proteins play multiple roles in plant development [39,40]. In Arabidopsis, rice and in some other plants, miR156 regulates leaf development by targeting Squamosa-Promoter Binding protein-like (SBP) transcription factors [31,39]. The identification of SBP as a target of gma-miR160 may indicate the additional level of regulation for SBP during soybean seed development.
Our analyses of the early, mid and late maturation developmental stages of soybean show a number of targets similar to those found by Song et al., [21] in the very young seeds of the cultivar Heinong44 including the SPB transcription factors. One notable difference was the absence in our degradome data of miRNA172 targets which include members of the AP2 transcription factor family. From inspection of sequenced small RNA populations from the 50-75 mg seed coats and cotyledons of Williams [27], we find only a few occurrences of the miR172 family (less than 30 occurrences per million reads) while some family members of the miR156 family are highly abundant (99,000 per million) in the cotyledons. We speculate that miR172 and/or its targets may be more abundant in the very young seed used by the Song et al. group [21] and not prevalent in the midmaturation seed that we have examined. In Arabidopsis, miR172 has been reported to be involved in the regulation of flowering time and floral development [40]. Alternatively, the AP2 factors may not be detected as targets in the degradome data if translational repression by miR172 is operative as has been shown in Arabidopsis flower development [41]. Nuclear Factor Y (NF-Y) was shown to control a variety of agronomically important traits, including drought tolerance, flowering time, and seed development [42]. We detected seven NF transcription factor YA subunit mRNAs specifically in seed coats that are targets of miR169 family members that occur in both seed coats and cotyledons (Additional file 1). These targets may indicate some specific regulation of NF-YA transcription factors during soybean seed development.
To obtain a deeper understanding of soybean seed development, we investigated tissue specific miRNA target identification in the cotyledons and seed coats at different seed developmental stages. Based on the degradome data, we identified some miRNAs that may act differentially in the cotyledons versus the seed coats to degrade their targets (Tables 3 and 4). F-box proteins involved in auxin-stimulated protein degradation (TIR1-like) were among the identified targets specifically found in soybean cotyledon ( Table 3). The LRR kinases have been reported to play important roles in plant development and brassinosteroid and ABA signaling [43]. We identified several LRR domain containing proteins as targets for gma-miR393, gma-miR1523 and gma-miR2109 (Tables 2, 3 and 4). The presence of these miRNA targets implies their regulation during soybean seed development. As a storage organ, the soybean seed contains significant amounts of lipid and protein. Thus the regulation of energy metabolism is very important during seed development. We identified a number of targets in the soybean cotyledon such as NADP/FAD oxidoreductase, ribose-5-phosphate isomerase, GTPase activating proteins and ferredoxin related proteins which are related to energy metabolism. Both in the soybean cotyledon and seed coat, we found pentatricopeptide repeat (PPR) proteins as targets of miR1520 (Table 2 and Additional file 1) which regulates gene expression in the mitochondria and chloroplasts [44,45]. Since we constructed our degradome libraries using cotyledons and seed coats from different seed developmental stages, we identified targets of miRNAs during a broad range of soybean seed development.
Auxin is an important phytohormone in higher plants. It acts as a key player in plant development [46]. As the transducer of auxin signaling, ARFs play vital roles in plant development, including shoot, root and flower formation [47,48]. In Arabidopsis, miR160 and miR167 are involved in auxin signaling via regulation of ARF genes [1]. In rice, a number of ARF encoding genes have been identified which are regulated by osa-miR160 and osa-miR167, respectively [33,34]. In our study, we identified a large number of ARF genes as targets for different miRNAs such as gma-miR160 and miR167. In the cotyledon (100-200 mg) degradome library (Additional file 1), we identified five targets annotated as Auxin Response Factors (ARFs) for gma-miR160, and four Williams data sets from this study were generated using five degradome libraries of early, mid and late soybean seed developmental stage. The Heinong44 data set represented an earlier seed stage [21]. of these, Glyma12g08110.1, Glyma12g29720.1, Gly-ma14g33730.1 and Glyma11g20490.1, were validated by RLM-5'RACE showing precise cleavage as expected ( Figure 2). These results suggested that gma-miR160 could participate in auxin signaling via down-regulation of ARFs during soybean seed developmental stages. The cleaved mRNAs captured by the degradome procedure indicate that the levels of the ARF mRNA targets are likely to be decreased, but qRTPCR or RNA sequencing data would be needed to directly confirm the effect on mRNA levels for a particular ARF target gene.
In our degradome libraries, miRNA targets are involved in major transitions between each stage of seed development and transcription factors account for approximately half of these targets. Of 183 identified targets in our soybean degradome libraries, GO analysis for biological function indicates that these genes are mainly involved in developmental and metabolic processes ( Figure 3 and Additional file 2). Enrichment of developmentally related genes as target miRNAs suggests the high level of regulation of gene expression during soybean seed development. The larger number of targets found in the 300-400 mg desiccating, yellow cotyledons of late maturation implies that post-transcriptional regulation by miRNAs may aid in shifting the developmental program of the immature soybean cotyledons from biosynthesis of storage reserves to a catabolic role in utilization of those reserves during seed germination and growth. The miRNA targets verified by degradome sequencing will provide useful information for understanding and revealing significant roles of miRNAs during soybean seed development.

Conclusion
Degradome sequencing is a valuable tool for the experimental confirmation of miRNA targets in higher plants. This method can reveal additional targets which are difficult to identify by computational prediction alone and confirm that the targets genes have been cleaved in specific tissues. Five degradome libraries from three different developmental stages identified 183 miRNA targets. Identification of soybean seed coat and cotyledon specific miRNA targets gives better understanding of tissue specific miRNA targets during seed development. In summary, the current study has confirmed a large set of targets that are subjected to miRNA guided degradation including many transcription factors and a surprisingly large number of targets in the late stages of cotyledon development. The data provides an avenue to explore more details about developmental stage specific miRNA targets that play critical roles in each of the important tissues during seed development.

Plant materials
Soybean (Glycine max cv. Williams) plants were grown in a greenhouse and seeds were collected at different developmental stages including early maturation, green 25-50 mg fresh weight seed, mid-maturation green 100-200 mg, and late maturation yellow 300-400 mg fresh weight seed. Immediately, cotyledons and seed coats were separated by dissecting whole seeds and then frozen in liquid nitrogen. Subsequently the tissue was freeze dried and stored at −80°C.