Genome-wide identification and characterization of novel microRNAs in seed development of soybean

ABSTRACT

MicroRNAs (miRNAs) are important and ubiquitous regulators of gene expression in eukaryotes. However, the information about miRNAs population and their regulatory functions involving in soybean seed development remains incomplete. Base on the Dicer-like1-mediated cleavage signals during miRNA processing could be employed for novel miRNA discovery, a genome-wide search for miRNA candidates involved in seed development was carried out. As a result, 17 novel miRNAs, 14 isoforms of miRNA (isomiRs) and 31 previously validated miRNAs were discovered. These novel miRNAs and isomiRs represented tissue-specific expression and the isomiRs showed significantly higher abundance than that of their miRNA counterparts in different tissues. After target prediction and degradome sequencing data-based validation, 13 novel miRNA–target pairs were further identified. Besides, five targets of 22-nt iso-gma-miR393h were found to be triggered to produce secondary trans-acting siRNA (ta-siRNAs). Summarily, our results could expand the repertoire of miRNAs with potentially important functions in soybean.

Graphical Abstract

Identification of novel microRNAs and their targets in soybean.

KEYWORDS:

• miRNA
• degradome sequencing
• ta-siRNA
• Glycine max

MicroRNAs (miRNAs) are a class of ~ 21nt non-coding RNAs, which are well known as important regulators ubiquitously existing in both plants and animals. There are increasing evidences that they take part in a wide range of biological processes in eukaryotes, such as organogenesis [1–3], signal transduction [4] and stress responses [1,5–8]. In plant, a primary miRNA (pri-miRNA) is recognized and processed by DICER-LIKE1 (DCL1) to produce a miRNA precursor, which is further cropped by DCL1 to yield a single small RNA (sRNA) duplex, miRNA-miRNA*, with two-nucleotide 3ʹ overhangs [9]. One strand of the duplex is the mature miRNA, it loaded onto an Argonaute (AGO) protein to form a RNA-induced silencing complex (RISC) to exert the cleavage of target transcripts. The other strand, miRNA*, is normally rapidly degraded [10]. Some miRNAs have been revealed capable of triggering phased siRNA (phasiRNAs) biogenesis from their targets. The products of miRNA-guided cleavage can further serve as the substrate for RDR6, generating a dsRNA precursor which is further processed by DCL4 to produce phasiRNAs [11]. Some phasiRNAs have been revealed function in trans regulation of target genes via cleavage at complementary target sites, just like many plant miRNAs. These phasiRNAs were named as trans-acting siRNAs (ta-siRNAs) [12–14].

Soybean is the most economically important legume in the world. It produces seeds that are enrich with proteins and oils which provide sources of energy and nutrition for both humans and animals. A lot of nutritional contents are accumulated in the soybean seeds during their development. Thus, understanding the regulation of metabolic pathways during seed development is one of the prerequisites for crop improvement [15].

MiRNAs have been revealed to play an essential role in seed development [15,16]. In recent years, increasing investigations focused on miRNA discovery and their regulatory function analysis contributed greatly to the quality improvement of soybean seeds, such as increasing oil production by controlling fatty acid pathway [17], and surviving in abiotic stress environment [18,19].The recent advance in next-generation small RNA sequencing technology coupled with developments of bioinformatics tools provided a powerful way for genome-wide discovery of considerable miRNAs with important regulatory functions in plants [15,20–25]. To date, 639 mature miRNAs and 573 miRNA precursors of soybean have been registered in the miRBase (release 21). However, there are increasing evidences that the information of miRNAs population and their regulatory functions in soybean seed development remains incomplete [26–28].The investigation of discovering novel miRNA with potential regulatory functions during soybean seed development period is necessary to continue, which will benefits the fully understand of miRNA regulatory network, seed quality controlling and improvement of soybean.

In previous study, we found DCL1-mediated cleavage signals for the processing of the miRNA precursors could be detected by degradome sequencing. Besides, their positions are highly correlated with the boundaries of the miRNAs (−5p)/miRNAs*(−3p). So, we constructed a bio-software, miRNA Digger (Figure 1), for genome-wide mining of miRNAs by using the degradome-supported Dicer cropping signals as the markers to detect potential miRNA loci [29].

Figure 1. Workflow of miRNA Digger. miRNA Digger is free for download from http://www.bioinfolab.cn/miRNA_Digger/index.html [29].

In this study, 12 sRNA high-throughput sequencing (HTS) data sets from six tissues of cotyledon stage soybean seed were input into miRNA Digger for miRNAs mining, followed by the expression specificity analysis and target identification. A number of novel soybean miRNAs and their targets were identified and some miRNAs registered in miRBase (Release 21) were reevaluated.

Materials and methods

Data sources

The Williams 82 soybean genome sequences, transcripts sequences and annotation information were retrieved from Phytozome. The sRNA HTS data (Supplementary Table S1a) and the degradome sequencing data libraries (GSM647200,GSM825574, GSM848963,GSM848964,GSM848965,GSM848966 and GSM848967) employed for miRNA screen were public available data sets from GEO (Gene Expression Omnibusncbi.nlm.nih.gov/geo/).

Pre-treatment of high-throughput sequencing data

In order to enable the cross-comparison between different libraries, all the HTS data sets were employed for pre-treatment. The pre-treatment procedure discarded the short sequences containing ‘N’s or with ‘0ʹ read counts, and normalized the read count of each short sequence with reads per million (PRM) as described previously [29].

Genome-scale screen of miRNA-miRNA* duplexes and their corresponding precursors

miRNA Digger reports miRNA-miRNA* duplexes and their corresponding precursors based on the following five filtering procedures: (1) The genomic sequences were subdivided into 10,000-nt fragments with 500-nt overlaps. (2) The degradome tags are then mapped onto both sense and antisense strands of those 10,000-nt fragments to obtain the degradome signatures and their locations. The averaged read count (in RPM) of a degradome signature which mapped to a potential cleavage site from one library should be ‘five times or more [30,31] than the surrounding degradome tags, and it also should be among the top 12 most abundant degradome tags that perfectly matched to the corresponding 10,000-nt fragment in the same library [32]. (3) The length of the longest pre-miRNA (the parameter “L”) was set as 700nt. And based on a single potential cleavage locus, two pre-miRNA candidates of 730 nt in length were obtained as pre-miRNA candidates, as described previously [29]. (4) The pre-miRNA candidates were further employed for secondary structure prediction and followed by base complementarity analysis. Two 30nt possible “Stem” structure regions of a single miRNA precursor candidates around the cleavage locus were analyzed, and the miRNA precursor candidates with hairpin-like structures and possess complementarities higher than 70% were retained. (5) The sRNAs from HTS data sets were then mapped onto the miRNA precursor candidates and retained the potential pre-miRNAs with detectable miRNA and miRNA* clusters. The most abundant sRNA within a miRNA cluster was treated as mature miRNA, and the corresponding miRNA* was extracted based on the structurally featured 2-nt 3ʹ overhangs of the miRNA-miRNA* duplex [31]. The most abundant isomiR* was treated as the counterpart of mature miRNA, if no canonical miRNA* is detected within the miRNA* cluster, which might causes by the rapid degradation of miRNA* after DCL1-mediated processing.

Expression analysis of miRNAs across different tissues

Twelve sRNA HTS data sets (Supplementary Table S1a) were employed for expression analysis. There were two duplicated data sets for every tissue sample, therefore, the average abundant (PRM) of a single miRNA was calculated and used for analysis. The log2 transformed values of abundance were displayed as heat maps. A previous reported toolkit: Heml, was utilized for heatmap profiling the differential expressions of miRNAs [33].

Target prediction and identification

The potential targets of miRNAs and ta-siRNAs were predicted by using miRU algorithm, which could thoroughly search for potential complementary sites on cDNA sequences with mismatches tolerable when a miRNA/ta-siRNA sequence was given, with default parameters [34]. Then, in order to get the genuine regulation of “miRNAs/ta-siRNA–target”, the degradome sequencing data (or Parallel Analysis of RNA Ends (PARE) libraries) were employed for verification of miRNA/ta-siRNA mediated cleavage sites on the predicted targets. Only when the cleavage signal resides within the middle (10 to 11nt away from the 5ʹ-end) of miRNA/ta-siRNA binding site, can the target be considered as the real target of the miRNA/ta-siRNA. Finally, target plots (T-plots) were performed to visualize the match results of the degradome data on targets [35–38].

PhasiRNA identification

Phased duplexes with 2-nt overhang at 3ʹ-end and 21–22 nt length of the target transcripts were computationally predicted by “phase processing” in double strand form from the miRNA-mediated cleavage site (termed as register 1) and other forward registers. Then each of these phased duplexes was separated into a plus siRNA and a minus siRNA, and followed by read count searching in the 12 sRNA HTS data sets. PhasiRNA searching was performed using a five-cycle sliding window, the phasiRNA production region should contain at least 5 tandem “processing” duplexes, and each of these duplexes should contain at least one phasiRNA [39].

Results

Identification of novel miRNAs in soybean seeds

In order to discover miRNAs, the soybean genome (version 1.1) and 12 sRNA HTS libraries from six tissues: Embryo Proper (EP), Embryo Proper Cotyledon(EPC), Seed Coat Epidermis (SCE), Seed Coat Hilum (SCH), Seed Coat Inner Integument (SCII) and Seed Coat Outer Integument (SCOI) of soybean at cotelydon-stage (Supplemental Table S1a) were input into our previously constructed bio-software, miRNA Digger [29]. Five filtering steps were automatically implemented to screen hairpin structure miRNA precursors containing both degradome-supported Dicer cropping sites and detectable miRNA-miRNA* duplexes (see details in Materials and methods). As a result, miRNA Digger reported totally 62 miRNA-miRNA* duplexes and their corresponding precursors (Supplemental Table S1b,c,d). Among them, 16 novel miRNAs were discovered from 16 novel precursors (Figure 2 and Supplemental Table S1b). All these novel pre-miRNAs were classic “stem-loop”, and each of them possessed at least one degradome-supported DCL1 cleavage signal at the boundary of the miRNA-miRNA* duplex loci. Interestingly, another novel miRNA, gam-miR482a-5p.2, was identified from different locus on the same precursor of gam-miR482a-5p and gam-miR482a-3p which had been registered in miRBase (Figure 3 and Supplemental Table S1b). Recently, Zabala et al. reported that gma-miR482a-3p represented higher expression levels than gma-miR482a-5p in various tissues of soybean [40]. Besides, gma-miR482a-3p was found to be conservatively expressed in several legume and non-legume dicots, to target resistance gene receptor kinases related to defense responses and to play important role during initial stages of nodule organogenesis [41,42]. However, in our study, the expression analysis showed that the novel gam-miR482a-5p.2 represented a significantly high expression count in all six different tissues of soybean seed while gam-miR482a-5p and gam-miR482a-3p were almost unexpressed in all the six tissues (Figure 3). Besides, 31 known miRNAs (Supplemental Table S1d) and 14 isomiRs, the variants of microRNAs，(Figure 4 and Supplemental Table S1c) were also identified by miRNA Digger in this study.

Figure 2. Novel miRNA-miRNA* duplexes discovered in soybean. A-P: The newly identified miRNAs and miRNA*s (marked with lines) were mapped to the pre-miRNAs. The degradome-supported DCL1-mediated cleavage signals were marked with arrows.

Figure 3. Novel miRNA, gam-miR482-5p.2, uncovered from known pre-miRNA. A: The newly identified gam-miR482-5p.2/gam-miR482-5p.2* (marked with lines) and miRBase registered gam-miR482-5p/gam-miR482-3p (marked with dotted lines) were mapped to the pre-miRNA. The degradome-supported DCL1-mediated cleavage signal was marked with arrow. B: Expression of gam-miR482-5p.2, gam-miR482-5p.2*, gam-miR482-5p and gam-miR482-3p in six tissues.

Figure 4. isomiRs uncovered in known precursors. The miRNA Digger discovered isomiRs/isomiRs* (marked with lines) and miRBase registered miRNAs counterparts (marked with dotted lines)were mapped to the pre-miRNAs. The degradome-supported DCL1-processing sites were marked with arrow.

Expression analysis of novel miRNAs across different tissues

An assessment of differential expression was performed for the novel miRNAs across six tissues. The results indicated that, 16 novel miRNAs showing tissue-preferential accumulation (Figure 5 and Supplemental Table S2). Most of them represent enrichment in more than two tissues; i.e., gma-miRNA-15 was observed strongly present in EP and EPC but showed low expression levels in other tissues. Gma-miRNA-5 was found strongly present in all tissues, however, it presented extremely high abundance (2900.81 RPM) in SCOI, which is 1.3 to 4.2 times of that in other tissues.

Figure 5. Expression profiling of novel miRNAs. The expression pattern of novel miRNAs in cotyledon stage across six different tissues, including: Embryo Proper (EP), Embryo Proper Cotyledon (EPC), Seed Coat Epidermis (SCE), Seed Coat Hilum (SCH), Seed Coat Inner Integument (SCII) and Seed Coat Outer Integument (SCOI) were analyzed and the log2 transformed values of abundance are displayed as a heat map.

As miRNA Digger only extracted the most abundant sRNA sequence in the miRNA cluster as the mature miRNA and 14 isomiRs of known miRNAs (miRBase release 21) were reported as mature miRNAs in this study, the enrichment abundance of isomiRs and their miRNA counterparts in different tissues were also analyzed. As shown in Figure 6 and Supplementary Table S3, all these isomiRs represented tissue specific expression and showed significantly higher accumulation than their miRNA counterparts in different tissues, i.e., iso-gma-miR1510a-5p, iso-gma-miR1510b-5p, iso-gma-miR2118a-3p and iso-gma-miR4996 showed significantly high enrichment in all six tissues, about 6 to 1100 times of that of their miRNA counterparts (Supplementary Table S3a). Besides, in other two stages of seed development, globular and heart stage, almost all of these isomiRs present higher abundance than their miRNA counterparts (Supplementary table S3b and c). Some annotated miRNAs in miRBase showed very low enrichment in seed development stages, e.g. gma-miR393h, gma-miR482a, gma-miR1507c-3p, gma-miR5037a and gma-miR5037c (Figure 6 and Supplementary Table S3). As they have high sequence similarity, isomiRs and miRNA are usually predicted to target the same mRNAs. So, the isomiRs identified in this study which had much higher expression levels may play primary regulatory roles in soybean than their miRNA counterparts.

Figure 6. Expression profiling of isomiRs and their miRNA counterparts. The expression pattern of newly discovered isomiRs and their miRBase-registered mature miRNA counterparts in cotyledon stage across six different tissues, including: Embryo Proper (EP), Embryo Proper Cotyledon (EPC), Seed Coat Epidermis (SCE), Seed Coat Hilum (SCH), Seed Coat Inner Integument (SCII) and Seed Coat Outer Integument (SCOI) were analyzed and the log2 transformed values of abundance are displayed as a heat map.

miRNA target verification using PARE libraries

All the novel miRNAs were employed for target prediction using miRU algorithm, followed by degradome-based verification using public PARE data and T-plot profiling. The principle of target identification were described in the “materials and methods”. As a result, 2 novel miRNAs involving in 13 novel miRNA-target interaction were identified. The novel miRNAs cleave their targets all in the position of 10-11nt, which fits the characteristic of canonical miRNA very well.

As shown in Table 1 and Figure 7, gma-miRNA-1 targets three SQUAMOSA PROMOTER BINDING PROTEIN-LIKE6 (SPL6) genes (Glyma02g13370.1, Glyma05g00200.1 and Glyma17g08840.1), four SPL9 genes (Glyma02g30670.1, Glyma03g29900.1, Glyma18g36960.1 and Glyma19g32800.1), one SPL2 gene (Glyma11g36980.1) and one SPL3 gene (Glyma07g31880.1). SPL2 has been revealed controlling proper development of lateral organs in association with shoot maturation in the reproductive phase [43], SPL3 functions in the control of flowering time and phase change [44], SPL6 has been revealed essential for the resistance of pathogens and virus in Arabidopsis and Tobacco, respectively [45] and SPL9 has strong effects on leaf initiation rate [44]. These results imply that gma-miRNA-1 is involved in the regulation of soybean development and pathogen resistance. Gma-miRNA-7 was detected targeting two genes (Glyma14g16730.1 and Glyma17g29790.1) encoding DUF581 containing proteins. DUF581 domain has been revealed encoding plant specific FCS-Like Zinc Fingers and involving in protein-protein interaction. Many DUF581 family proteins have been identified in A. thaliana, and one of them has been revealed involving in ABA-mediated seed dormancy and induced during senescence [46–49].

Table 1. Targets of novel miRNAs.

miRNA Sequence	cDNA ID	miRU start-ending	PARE cut sites	Gene annotation
gma-miRNA-1	Glyma02g13370.1	1209–1229	1220	squamosa promoter-binding protein-like 6 (SPL6)
	Glyma02g30670.1	758–778	769	SPL9; transcription factor
	Glyma03g29900.1	1138–1158	1149	SPL9; transcription factor
	Glyma05g00200.1	1192–1212	1203	squamosa promoter-binding protein-like 6 (SPL6)
	Glyma07g31880.1	601–621	612	SPL3; DNA binding/transcription factor
	Glyma11g36980.1	1213–1233	1224	SPL2; DNA binding/transcription factor
	Glyma16g05900.1	315–335	326	unknown
	Glyma17g08840.1	1001–1021	1012	squamosa promoter-binding protein-like 6 (SPL6)
	Glyma18g00890.1	239–259	250	unknown
	Glyma18g36960.1	713–733	724	SPL9; transcription factor
	Glyma19g32800.1	737–757	748	SPL9; transcription factor
gma-miRNA-7	Glyma14g16730.1	583–604	595	DUF581 family protein
	Glyma17g29790.1	738–759	750	DUF581 family protein

Figure 7. Degradome sequencing data-based validation of the miRNA-target interactions in soybean. Five libraries of degradome sequencing data libraries (GSM647200,GSM848963,GSM848964,GSM848965 and GSM848967) were recruited for T-plot profiling. The IDs of the target transcripts and the corresponding miRNAs are listed on the top. The y axes measure the normalized reads (in RMP, reads per million) of the degradome signals, and the x axes represent the position of the cleavage signals on the target transcripts. The binding sites of the phased siRNA on their target transcripts were denoted by horizontal lines, and the dominant cleavage signals were marked by arrows.

According to the previous reports, 22-nt miRNAs have the ability to trigger the ta-siRNAs production from the cleaved targets [50]. So the twelve sRNA HTS data sets were employed for 22-nt novel miRNAs/isomiRs triggered secondary ta-siRNAs identification. The results showed that only five targets of iso-gma-miR393h generated secondary phasiRNAs. Based on the modified ta-siRNA biogenesis model which we have reported previously, 135 novel phasiRNAs were discovered (Supplementary Table S4). All these phasiRNAs were then employed for target identification, and 3 novel ta-siRNA-target interactions were identified (Table 2 and Figure 7).

Table 2. Targets of miRNA-triggered secondary ta-siRNAs.

Ta-siRNA ID	Ta-siRNA sequence	Target cDNA ID	miRU Start-ending	PARE cut sites	Gene annotation
Glyma02g17170.1(1742)21 3ʹD11(+) Glyma10g02630.1(2167)21 3ʹD10(+)	CUGCAUUAAGGAUUUCUUGAG	Glyma05g03930.1	564–584	575	Class I peptide chain release factor
Glyma19g39420.1(1752)21 3ʹD8(+)	GGUUGGACAUGCCUGGUUUUG	Glyma01g00720.1	629–649	640	Regulator of Vps4 activity in the MVB pathway protein

The ta-siRNA, Glyma19g39420.1 (1752)21 3ʹD8(+), was found targeting a gene (Glyma01g00720.1) which encoding a regulator of Vps4 activity in the Multivesicular Body (MVB) pathway. In recent years, the MVBs have been revealed playing important roles in mediating protein traffic and plant immune responses in plants [51]. Although the regulatory mechanism of Vps4 and MVBs in soybean are not well studied, however, Shahriari et al. revealed that the homologous of Vps4 in Arabidopsis, AtSKD1, related to the cell surviving during seed development [52]. Therefore, we speculated that Glyma19g39420.1(1752)21 3ʹD8(+) might involving in regulation of mediating protein traffic, pathogen response and cell surviving during seed development of soybean.

Another two ta-siRNAs, Glyma10g02630.1 (2167)21 3ʹD10(+) and Glyma02g17170.1(1742)21 3ʹD11(+), were found targeting a Class I peptide chain release factor gene (Glyma05g03930.1), indicating that they are both play important roles for cell viability [53].

Discussion

Soybean is an important economic crop. In recent years, a number of miRNAs in soybean have been revealed playing crucial regulatory roles in organ development [53], nutrient signaling [45], biotic and abiotic stress responses. In this study, a previously constructed bio-software, miRNA Digger, was used for miRNA mining from 12 public available sRNA data sets. As a result, 17 novel miRNAs were discovered. Two of them were further identified involving in 13 novel miRNA-target interactions and might play important roles in regulation of protein–protein interaction, stress response, plant development, plant proliferation and differentiation.

In addition, 14 isomiRs rather than their miRNA counterparts were discovered in this study. Recently, literatures revealed that, a single pre-miRNA can generate a series of miRNA variants (termed as isomiRs) with various length and 5ʹ or 3ʹ ends [54,55], and different variants of the same miRNA accumulated divergently in different tissues and development stages. However, the isomiRs uncovered in this study represented significantly higher abundance than those registered in miRBase (Release 21) in all three stages of seed development (Supplementary Table S3), which indicated that these isomiRs might play primary regulatory roles than their miRNA counterparts in soybean. Although numerous miRNAs from variant species have been discovered and registered in miRBase, some incorrect annotation of miRNAs in miRBase have been reported [1,56,57]. So, using different methods to reevaluate the annotated miRNAs is necessary and important for fully understanding of the miRNA family and their functions.

The 22-nt size class miRNAs have been drawn a great attention in recent years because of their potential ability of triggering the production of secondary ta-siRNAs which have multiple biological functions [12,58,59]. Therefore, all the target transcripts of 22-nt miRNAs/isomiRs were employed for secondary ta-siRNAs identification. However, except those of iso-gma-miR393h, the other target transcripts could not enter into the synthesis pathway for secondary ta-siRNA generation. This anomaly indicates that the current ta-siRNA triggering models are imperfect and needs to be further investigated.

Author contributions

L.Y. accomplished the identification of novel miRNAs and phasiRNAs; R.G. accomplished miRNAs expression analysis and miRNAs/ta-siRNAs target verification. L.Y. and R.G. wrote the main manuscript. L.Y, R.G, Y.M. and C.S. designed the experiments. Y.J., X.Y. and Z.Y. contributed to collection and pre-treatment of the data, and preparation of supplemental files. Y.M. and C.S. revised the manuscript. All of the authors reviewed the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China, under Grants[numbers:31801102,31771457,31571349]; Zhejiang Provincial Natural Science Foundation of China, under Grants [numbers:LQ16C060001, LY17C190001], Public Welfare Technology Application Research Project of Zhejiang Province, under Grant[number 2016C33193], and Scientific Research Fund of Zhejiang Provincial Education Department, under Grant [number Y201533314].