Construction and analysis of degradome-dependent microRNA regulatory networks in soybean

Background Usually the microRNA (miRNA)-mediated gene regulatory network (GRN) is constructed from the investigation of miRNA expression profiling and target predictions. However, the higher/lower expression level of miRNAs does not always indicate the higher/lower level of cleavages and such analysis, thus, sometimes ignores the crucial cleavage events. In the present work, the degradome sequencing data were employed to construct the complete miRNA-mediated gene regulatory network in soybean, unlike the traditional approach starting with small RNA sequencing data. Results We constructed the root-, cotyledon-, leaf- and seed-specific miRNA regulatory networks with the degradome sequencing data and the forthcoming verification of miRNA profiling analysis. As a result, we identified 205 conserved miRNA-target interactions (MTIs) involved with 6 conserved gma-miRNA families and 365 tissue-specific MTIs containing 24 root-specific, 45 leaf-specific, 63 cotyledon-specific and 225 seed-specific MTIs. We found a total of 156 miRNAs in tissue-specific MTIs including 18 tissue-specific miRNAs, however, only 3 miRNAs have consistent tissue-specific expression. Our study showed the degradome-dependent miRNA regulatory networks (DDNs) in four soybean tissues and explored their conservations and specificities. Conclusions The construction of DDNs may provide the complete miRNA-Target interactions in certain plant tissues, leading to the identification of the conserved and tissue-specific MTIs and sub-networks. Our work provides a basis for further investigation of the roles and mechanisms of miRNA-mediated regulation of tissue-specific growth and development in soybean. Electronic supplementary material The online version of this article (10.1186/s12864-019-5879-7) contains supplementary material, which is available to authorized users.


Background
MiRNAs, a class of~21 nt non-coding endogenous small RNAs, play crucial roles in soybean growth, development and stress responses [1][2][3][4][5][6] by pairing to the target mRNAs to direct their post-transcriptional repression [7]. MiRNA itself does not have an open reading frame (ORF) or encode any protein, but it has a high degree of evolutionary conservatism among all kind of species and obvious different expression along with time and space transmission, showing its important statute in various biological processes of plants.
Since the first miRNA Line4 was discovered in larval development of nematode [8,9], the discoveries of novel miR-NAs were explosive with small RNA sequencing technique, providing influx of information for researchers. Then, the scholars started to explore the interactions and regulatory networks medicated by two or more miRNAs with the prosperity of next generation sequencing, all kinds of small RNA databases and bioinformatics tools, for instance, miRBase [10], psRNATarget service [11], bowtie [12] and Cytoscape 8 [13].
With such abundant small RNA datasets and tools, it is allowed to apply different patterns of biosynthesis into small RNA quantification and their expression detection. MiRNAs are generated from stem-loop pre-miRNAs while siRNAs are mainly exacted from double-stranded precursors. According to the latest criteria of miRNA annotation updated by Axetell and Meyers on plant cell [14], its length is limited to 20~24 nt, and mature 21-nt miRNA discomposes from double-strand RNA and form miRISC (miRNA-induced silencing complex) with Argonaute protein(AGO) family to repress translation of mRNAs or degrade mRNAs in the seed region [7,15,16]. The deep sequencing results of cleaved targets, generated latter, which is known as degradome, are used to validate the authenticity of the predicated MTIs.
Since the technique of degradome sequencing was performed to identify the miRNA-target relations in combination with rapid amplification of cDNA ends (RACE), high-throughput sequencing techniques and bioinformatics analysis [17], and mature protocols of Gene Ontology (GO) and Pathway analysis, we can understand thoroughly the regulation by miRNAs and it supports the results of bioinformatics and experimental supplement.
Simultaneously, some miRNAs have been reported with their specific regulation on different target genes in various biological pathways in plants. However, they tended to focus on miRNA expression performance and just with the predicted miRNA-Target relations, which usually results in incomplete regulatory relations. MiR172c was found to involve the repression of Auxin protein 2 (AP2) and to modulate the nod-related transcription factors NNC1 signaling pathway associated with nodule initiation in soybean [18]. Xu constructed miRNA-mediated gene regulatory networks in soybean cyst nematode (SCN) within 32 miRNA families and found six of them regulate the formation of SCN in root [19]. Even some genome-wide miRNA investigations, their explorations are also limited to expression verification. For instance, Chen did an investigation in soybean beginning with sRNA differential expression analysis under SMV infection, to find process-specific miRNAs, but most of the differentially expressed miRNAs did not have the high-efficiency cleavages based on degradome verification [20], which possibly dropped vital processes. The similar phenomenon showed up in Volkdin's experiments about seed development [21].
MicroRNAs mediated gene regulatory networks generally started with investigating the miRNA's expression profiling and target predictions so far. However, the level of miRNA's expression cannot represent the level of cleavage, thus, such analysis focusing on expression level often ignores the crucial cleavage events corresponding with low-level expressed miR-NAs. In the present work, reversely, the degradome sequencing data were employed to study the complete miRNA mediated regulatory networks in soybean. Based on degradome-dependent miRNA regulatory network (DDN) analysis, we found and validated conserved and specific MTIs, using tissue-specific MTIs to construct the root-, cotyledon-, leave-and seed-specific miRNA regulatory networks based on the degradome sequencing data, annotated with miRNA profiling analysis and discussed some co-interacted miRNA families mediated GRNs. All the results above indicate that the construction of DDN is a useful approach to describe a comprehensive regulatory network mediated by miRNAs, and will help botanist to learn further about the roles of miRNAs in soybean growth and development in the view of miRNA and target interactions.

Construction and description of global DDNs
We took raw tables of degradome-dependent networks as fundamental information to classify the regulatory types of different miRNA mediated subnetworks, and studied further about conserved MTIs and tissue-specific MTIs. Based on the data predicted with psRNATarget server and relations validated by degradome data, we constructed a degradome-dependent miRNA regulatory network (Fig. 1). The grey-green circles represent target, the grey-green squares represent target gene-encoding protein family (indicating biological function), blue nodes are expressionnon-changed miRNAs and orange ones are tissue-specific miRNAs. The width of edges is gradient from 1 to 7, representing the number of degradome libraries, which verified the same MTIs. The color of edges represents the category of cleavage of highest degradome count (red is Cat_1 level, and blue is Cat_2). It is obvious to observe that the combination of blue node and red edges accounting for large part of the verified networks. In the all-relation clusters, 9518 predicted relations were involved with 611 miRNAs, while 1804 of predicted MTIs were validated, involving 225 miRNAs.
Combining with identification of expression profile, we diagram occupation of differentially expressed miRNA (DEM) regulated MTIs in each tissue (Fig. 3): of 2561 MTIs verified in seed degradome, 60 are center with seed-specific DEMs; 13.6% MTIs (196 DEM-centered ones in all 1446 MTIs) verified by cotyledon degradome are DEM-center; 11.53% in leaf are mediated by leaf-specific DEMs, and none of the center miRNAs is root-specific in MTIs verified by root degradome. Difference between expression-nonchanged miRNAs and DEM, distribution difference among different tissues can both indicate the different activity of miRNA in soybean. Detailed workflow is present in Additional file 5: Figure S1, and results of expression verifications were shown as a heatmap in Additional file 6: Figure S2.
Tissue-specific miRNA-mediated regulatory networks Tissue-specific MTIs are particularly crucial to conclude the regulatory roles of these miRNAs in soybean. This approach only considers whether the MTIs are specific not the expression specificity of miRNAs. If degradomespecific MTI's degradome read count is over 10 TP10M and miRNAs can be detected in certain tissue, we described them as tissue-specific MTIs. While those miR-NAs whose expression in certain tissue is over twice of that in rest tissues (for each specific miRNA (SSR > Mean Square), we applied one-way ANOVA, when compared with the rest other tissues, if we observe log 2 (Expspecific /Exp other ) > 2, then this miRNA is considered as tissue-specific miRNA), were considered as tissuespecific miRNAs. Among 365 degradome-specific MTIs, there are 24 root-specific, 45 leaf-specific, 63 cotyledonspecific and 225 seed-specific MTIs (Additional file 3: Table S3), however, in 22 miRNAs from cotyledonspecific MTIs, only miR156f is cotyledon-specific; miR4996 is the only leaf-specific one of all 33 miRNAs from leaf-specific MTIs and miR156b is the only seedspecific miRNA among 112 miRNAs from seed-specific MTIs. There is no consistently specific-expression miR-NAs in root (Table 3). There are most abundant kinds of miRNAs from seed-specific MTIs, while distribution of miRNAs is similar in rest tissues.
We picked out tissue-specific subnetworks in every tissue--miR5674a and miR2109 co-mediated leaf-specific subnetworks, miR396 mediated cotyledon-specific subnetworks, miR164 mediated seed-specific subnetworks and root-specific networks, which were not tissuespecific miRNAs in certain tissue but had high efficiency of cleavage on targets (Fig. 6). Of these tissue-specific MTIs (listed in Additional file 3: Table S3), 45 are leafspecific, 63 are cotyledon-specific, 24 are root-specific and 225 are seed-specific.
In seed-specific networks, there are 225 seed-specific MTIs, including 3 seed-specific miRNA-center MTI-miR156b regulating Squamosa promoter-binding protein- like (SBP domain) transcription factor family protein (SPL), which is a typical miRNA-Transcription factor pattern. The rest seed-specific MTIs were mediated by nonspecific miRNA but had high degraded read count, and of these MTIs, We found some published miRNA-target pairs, such as miRNA2119-ADH (alcohol dehydrogenase 1), a flooding-stress-resistant gene [24](alcohol dehydrogenase 1), miR164-NAC (No Apical Meristem domain transcriptional regulator superfamily protein), droughtresponse transcriptional regulator [25], and miR1515dicer-like2 pairs. Seed-specific MTIs contain more abundant and high cleavage efficient regulatory relations between miRNA and target genes (Additional file 3: Table S3, Additional file 4: Table S4). Detailed information of tissue-specific MTIs is present in Additional file 3: Table S3, miRNAs expression level in Additional file 4: Table S4.

Interacted networks among different miRNA families in different tissues
By observing our constructed DDNs, we found an interfamilies regulatory network involving miR319 family and miR4996 (Fig. 7), since our DDNs were mainly separated as every single miRNA-family regulatory group. They shared the target Glyma.13 g219900, which is homologous to TEOSINTE BRANCHED 1, cycloidea and PCF (TCP) transcription factor 2 in Arabidopsis, where miR319regulated TCP has reported to control flowering time [26] . But their cleavage belongs to Cat_2 level, we hypothesized there might be a RNA edition of miRNA for higher probable cleavage on targets. So we found out the MTIs of miR319f/l and miR4996's isoforms, which still showed no Cat_1 level cleavage on Glyma.13 g219900. It indicated that Glyma.13 g219900 probably is regulated by both miRNA and an uncharacterized small RNA.

Difference analysis of networks regulated by the miRNAs from the same family
We observed that some miRNAs share totally different targets although they come from the same miRNA family. Different from the most miRNA family we detected, the members of miR167 family, as a typical example, do not work as a unit like others, but far apart into two general (Fig. 8), targets of gma-miR167k are mainly encoding auxin response factors, a regulator of plant growth and development from embryogenesis to A B C D Fig. 7 Interacted networks between different gma-miRNA families. a Co-interacted network of gma-miR 319 family and gma-miR4996, sub-network of Fig.3a; b Verified t-plot of gma-miR4996/Glyma.13 g219900, Y-aix represents degradome abundance(TP10M), X-aix represents nucleotide position (nt), the aligned region of microRNA on mRNA. c Verified t-plot of gma-miR319l/Glyma.13 g219900. d Verified t-plot of gma-miR319f/Glyma.13 g219900 senescence, GmARFs are reported to regulate the growth of soybean lateral root [27], other miRNA167 members repress zinc finger (C3HC4-type RING finger) family protein, nuclear protein X1 and HSP20-like chaperones superfamily protein, relating to abiotic stress response. To further dig out the possible molecular proof causing the huge functional difference, we compared the mature sequences of miRNA167 members and found that miR167k differentiate with other members in miR167 family from the 13rd base which is close to its seed region.

Discussion
Usually the miRNA mediated GRN is constructed by investigating the miRNA expression profiling and target predictions. However, miRNA's abundance sometimes cannot reflect whether these miRNAs, with high or low expression levels, play the crucial roles by cleaving the target genes in/ under certain tissues or treatments. Specifically, the higher/ lower level of miRNA does not mean the higher/lower level of cleavage. The same miRNAs may target different mRNAs in different tissues. Therefore, such analysis sometimes filters the crucial cleavage events where the corresponding miRNAs were lowly expressed. The degradome sequencing allows us to find out whether the cleavage events really occur and the cleavage abundances. Furthermore, whether the miRNA is expressed could be used to validate the degradome-dependent miRNA-target interaction in/under the certain tissues or treatments. This reverse approach may provide the reliable and complete miRNA-target interactions for network construction.

Tissue-conserved MTIs and tissue-conserved miRNAs
By comparison of different networks constructed for different tissues/treatments, we can finally identify the crucial miRNAs and the miRNA-target interactions responsible for specific biological processes in plants. Besides getting rid of heavy laboratory work of identification, our method is more focused on the network level and discusses more about miRNAs' cleavage efficiency in their MTIs . Some conserved miRNAs we identified have been reported. Li's team identified conservation and diversification of the miR166 family in soybean and discussed their potential roles [28], which we described as tissue-conserved miRNA family. Successively, some conserved MTIs were set as a solid model for further discussion. Turner reported miR160 regulated auxin responsive factor [29], Wang's research discussed miR167-GmARF regulatory interactions' function in soybean root development [27], and Xu found A C B Fig. 8 The gma-miR167 family regulatory network. a Gma-miRNA167 family regulatory networks (including miR 167a/b/c/d/e/f/g/j); b Gma-miRNA167k regulatory network; c Mature gma-miRNA167 sequence variation. The bases in red are sequence variations, the grey box is homologous part of gma-miRNA167 family, and the green box is homologous part of gma-miRNA167a/b/c/d/e/f/g/j that miRNA167 in soybean nodules and lateral root growth can regulate auxin response factor, but we also detected gma-miR167 meditated MTIs showed up and high expressed in the soybean seed(Additional file 2: Table S2) and target on zinc finger (C3HC4-type RING finger) family protein, which was related to stress response in rice and hasn't been reported in soybean yet.
However, tissue-conserved MTIs were not often regulated by conserved miRNA, and conserved miRNA also did not always in tissue-conserved MTIs. On the other side, their targets functioned in soybean development and disease resistance, indicating the conservatism of MTIs and miRNA families. Besides, some tissue-conserved miRNA families are identified to regulate totally different target genes in tissue-specific MTIs. For instance, miR171 is a conserved miRNA family in soybean, even miRNA171-GRAS regulatory pattern is conserved among different species, but in cotyledon, it has specific MTI with Leucine-rich repeat protein kinase family LRR-RPK) gene, a vital gene in regulation of leaf senescence.

Tissue-specific MTIs and tissue-specific miRNAs
In this investigation, we found a total of 156 miRNAs in tissue-specific MTIs, however, only 3 miRNAs were detected consistently specific expression (Table 3). Of 365 tissue-specific MTIs, 40 are mediated by tissue-specific miRNA, while only 5 pairs are regulated by consistent specific miRNA, containing 3 seed-specific pairs, one leafspecific pair and a cotyledon-specific one (seed-specific MTI is center with seed-specific miRNA) (Additional file 3: Table S3).
Tissue-specific miRNAs mediated MTIs we detected in this research only account for 8.8% of all tissue-specific MTIs. Tissue-specific MTIs required the high efficient cleavage (degradome count at least over 10 TP10M) of miRNA with different targets, so specific expression of miRNA is not necessary and tissue-specific miRNAs may have different tissue-specific MTIs for their different targets. MiR159e-3p, a cotyledon-specific miRNA and it targeted on MYB domain genes Glyma.04g125700 and Glyma.06g312900 with Cat_1 level cleavage, while in leaf it regulated Glyma.12g228100 and Glyma.13g271900, and had seed-specific MTIs in regulation of HCO3-transporter family genes Glyma.03g222300, Glyma.19g219500. Besides, these tissue-specific miRNAs may also get involved in conserved regulatory patterns, for example, leaf-specific miR156l regulated squamosa promoter binding proteinlike 9 (SPL9) in both cotyledon and leaf, and in cotyledon miR156l detected higher cleaved efficiency (dependent on level of degraded target read count) than that in leaf, this pattern has reported in soybean to regulate the plant architecture [33]. Some tissue-specific miRNA may have different tissue-specific MTIs in other tissue, such as, leafspecific miR1507c-3p has root-specific cleavage on Glyma.15g226100, encoding a LRR and NB-ARC domainscontaining disease resistance protein.
On the other hand, those non-specific miRNAs can participate in some tissue-specific MTIs, related to some crucial biological and biochemistry processes because of their specific targets, for instance, gma-miR5674a, a nonchanged miRNA, has a leaf-specific MTI with Glyma.09175800, homologous to NOP56-like pre RNA processing ribonucleoprotein in Arabidopsis. Of these MTIs, literally miRNA408c-3p-plantacyanin pair has been reported in Arabidopsis [34].
Co-regulated networks by cross-family miRNAs and the specific networks regulated by the miRNAs from the same family Overall our results, the distribution of 225 gma-miRNAs and their regulatory relations in different tissues are typically individuated into three conditions: (1) MiRNAs in soybean usually regulate as a unit of family, and their relations often perform as one-to-one, one-to-many or many-tomany form for regulating one or more target genes. The co-regulation between miRNAs often occurs in the same family. (2)Additionally, different miRNAs from different families have co-interaction with the unique genes, such as miR319f / g / l and miR4996 (Fig. 7). (3) Some miRNAs are isolated from other members in the same families, and can be used to correlate with multiple family members. It is independent of the main regulatory function group mediated by other family members, such as miRNA167k (Fig. 8) and other members in miR167 family (Fig. 8) and the causes for the differences may be the single base mutant in the conserved region or key sites on miRNA. According to the intersected networks, these universal MTIs validated in all the investigated tissues, and at least one member performed high expression in all tissues. These miRNA families are often considered as conserved miRNA families, often playing a key statue in the process of evolution and regulating the necessary biological processes, and generally have a cocleavage of the same target transcript in several tissues, indicating that there are often over one members participating in the key metabolic processes to ensure the stability.
However, here is still boundedness of this method: 1) these validated MTIs are limited to the number of degradome libraries, and the expression verification still requires expression profiling data. 2) The expression of tissue-specific MTIs is hard to identify. The solution we adopted here is to combine small RNA expression data and degradome-verified cleavage data. With the increase of degradome sequencing data, it would be more convenient to study the miRNA-target interaction and build the miRNA mediated GRNs with the full branches of the interactions.

Conclusions
The miRNA-target interactions (MTIs) and networks, not the miRNA themselves, should draw our attention when studying plant miRNA regulation. DDN approach may provide a more complete and reliable miRNA-target interactions, especially in the analysis of tissue/treatment-specific MTIs, and to study miRNA regulation in plants. Whether the abundances of MTIs can be directly predicted based on the abundance of the corresponding degradome reads is a question requiring further discussions and validations. In the present work, we constructed DDNs for four soybean tissues and identified the conserved and tissue-specific MTIs/sub-networks, which provides a basis of further studies of miRNA regulation in soybean growth and development.

Processing of the degradome sequencing data
We employed the Cutadapt [37] to remove the adapters of the degradome sequencing reads, and then unify the forms of raw data sets. Normalization of the read count was Reads per 10 Million (RP10M), which was calculated by dividing the raw read counts with the total genome-mapped read counts and then multiplied by 10 7 . The processed degradome data are available in the DPMIND [38] repository (http://cbi.njau.edu.cn).

MiRNA target prediction and degradome-based validation of miRNA-target pairs
We uploaded all the mature soybean miRNAs from miR-Base v.22 into psRNATarget prediction server (http:// plantgrn.noble.org/psRNATarget) [11] with default parameters and soybean cDNA library so that we could get the predicted results which describe the potential interactions between miRNAs and their target transcripts, and take these predicted results into the use of network construction for the observation of complicated regulatory interactions between miRNAs and target genes. First, the normalized degradome data should be mapped to the putative target sequences from prediction results; with the bowtie program, then the predicted miRNAtarget pairs were retained if they met the following criteria: (1) there must be at least one degradome sequence with their 5'ends resided within 9~12 nt regions away from the 5'ends of the target binding sites; (2) read counts of at least one degradome sequence in above region should be more than 5RP10M; (3) read counts of one degradome sequence in above region should be the most abundant (Category 1) or higher than median (Category 2), among all the reads mapped to one target. Finally, the t-plot figures were generated, using our local developed python script.
Classification of tissue-conserved and tissue-specific miRNA-target interactions (MTIs) We imported degradome-verified results above into Cytoscape [13] (http://cytoscape.org) to construct raw regulatory networks of miRNA and their targets for different soybean tissues. To classify MTIs, we merged all seven tables of Degradome-dependent MiRNA Regulatory Networks (DDNs) in order to pick out overlapped MTIs in all four investigated tissues, which are treated as tissue-conserved MTIs, while those one-degradome-verified MTIs whose center miRNAs were detected express in certain tissue, are considered as tissue-specific MTIs. Simultaneously, we annotated the target genes with Phytozome.v11 Gmax_275_ Wm82.a2.v1 (https://phytozome.jgi.doe.gov/pz/portal.html) [39] to prepare the discussion of DDN's function (Additional file 5: Figure S1).

Expression analysis of MiRNAs from degradome-classified MTIs
A home-made python script was in-house developed and was used to quantify miRNA expression in different small RNA sequencing libraries (libraries' information are listed in Table 1). Then we merged all miRNA expression files in four tissues into an expression matrix and detected tissue-differential expressed miRNA with R package DESeq2 [40](Additional file 4: Table S4).
We applied one-way ANOVA to identify tissuespecifically expressed miRNAs: For each specific miRNA (SSR > Mean Square), when compared with each rest tissue, if we can observe: log 2 Exp specific =Exp other > 2;