A Multi-Level Iterative Bi-Clustering Method for Discovering miRNA Co-regulation Network of Abiotic Stress Tolerance in Soybeans

Abstract

Although growing evidence shows that microRNA (miRNA) regulates plant growth and development, miRNA regulatory networks in plants are not well understood. Current experimental studies cannot characterize miRNA regulatory networks on a large scale. This information gap provides an excellent opportunity to employ computational methods for global analysis and generate valuable models and hypotheses. To address this opportunity, we collected miRNA-target interactions (MTIs) and used MTIs from Arabidopsis thaliana and Medicago truncatula to predict homologous MTIs in soybeans, resulting in 80,235 soybean MTIs in total. A multi-level iterative bi-clustering method was developed to identify 483 soybean miRNA-target regulatory modules (MTRMs). Furthermore, we collected soybean miRNA expression data and corresponding gene expression data in response to abiotic stresses. By clustering these data, 37 MTRMs related to abiotic stresses were identified, including stress-specific MTRMs and shared MTRMs. These MTRMs have gene ontology (GO) enrichment in resistance response, iron transport, positive growth regulation, etc. Our study predicts soybean MTRMs and miRNA-GO networks under different stresses, and provides miRNA targeting hypotheses for experimental analyses. The method can be applied to other biological processes and other plants to elucidate miRNA co-regulation mechanisms.

Keywords: abiotic stress tolerance in soybeans; bi-clustering; homology expansion; miRNA co-regulation; miRNA–target.

FIGURE 1

Flowchart of the authors’ research method.

FIGURE 2

miRNA-target regulatory modules (MTRM) iterative merge algorithm flowchart used to derive a gene interaction network.

FIGURE 3

Gene ontology (GO) term analysis of MTRM genes under various abiotic stresses (A–C). (A) GO semantic correlation analysis of abiotic stress, (B) drought stress, and (C) salt stress, and (D) the GO BP regulatory network of cooperative miRNAs under abiotic stresses. Triangles represent different miRNAs and circles represent different GOs. The size of the circle is determined by the number of genes contained in the GO in this article. The color of the circle depends on the representative GO. The areas with different colors show the modules obtained by our method.

FIGURE 4

Multi-level iterative biclustering results of soybean MTRMs. (A) Results under different iteration times at the 2 × 2 scale, (B) results under different iteration times at the 6 × 2 scale, and (C) the boxplot of the cluster score is calculated based on the gene ontology (GO) under the two scales when converging to a stable level, where based on the overall distribution, the results at the 6 × 2 scale are better; (D) shows the MTRM bicluster at level 1 before the 6 × 2 scale fusion, and (E) shows the corresponding MTRM bicluster at level 7 after the 6 × 2 scale fusion.

FIGURE 5

Gene ontology (GO) analysis of soybean MTRMs. (A) Semantic relevance of GO terms wherein the GO pathway has a certain concentration. (B) GO annotation enriched with 483 soybean MTRMs with enrichment results, which mainly involve positive regulation of development, heterochronic, chalcone biosynthesis, defense responses, and mitochondrial mRNA modification. (C) GO enrichment of the top five soybean MTRMs. The listed GO terms were enriched with significant p-values <0.00001.

FIGURE 6

Collected miRNA data on soybeans involved in various abiotic stress responses based on the data statistics from the literature. (A) The distribution of stress types in miRNAs where each vertical line represents one miRNA, and red is marked as relevant, (B) UpSet diagram (Lex et al., 2014) of modular genes under various abiotic stresses within the horizontal correspondence, where dots are used to refer to the corresponding cold stress, acid stress, heat stress, drought stress, and salt stress on the left. The point-to-point connection is realized longitudinally to indicate the intersection between the corresponding data sets, and the upper bar graph shows the number of genes in the intersection. In panel (C), the differentially expressed genes in each MTRM under abiotic stress are shown after screening. We used three indicators to filter the candidate clusters. According to the p-value, the related miRNA purity and the cluster score of each MTRM gene are placed under the corresponding stress. We selected the corresponding threshold, obtained the stress-related MTRMs with higher reliability, and marked them as red dots in Panel (C). Supplementary Figure 2 shows MTRMs under other types of stress.

FIGURE 7

Soybean MTRMs under various abiotic stresses. (A) Venn diagram of 37 kinds of soybean MTRMs under various abiotic stresses, including 14 drought-specific MTRMs, seven salt-specific MTRMs, two heat-specific MTRMs, and two shared MTRMs, (B) the miRNAs in the two shared MTRMs 31 and 493, and (C) a GO Treemap of 37 MTRMs under abiotic stress.