- Open Access
Differential microRNA expression, microRNA arm switching, and microRNA:long noncoding RNA interaction in response to salinity stress in soybean
BMC Genomics volume 23, Article number: 65 (2022)
Soybean is a major legume crop with high nutritional and environmental values suitable for sustainable agriculture. Noncoding RNAs (ncRNAs), including microRNAs (miRNAs) and long noncoding RNAs (lncRNAs), are important regulators of gene functions in eukaryotes. However, the interactions between these two types of ncRNAs in the context of plant physiology, especially in response to salinity stress, are poorly understood.
Here, we challenged a cultivated soybean accession (C08) and a wild one (W05) with salt treatment and obtained their small RNA transcriptomes at six time points from both root and leaf tissues. In addition to thoroughly analyzing the differentially expressed miRNAs, we also documented the first case of miRNA arm-switching (miR166m), the swapping of dominant miRNA arm expression, in soybean in different tissues. Two arms of miR166m target different genes related to salinity stress (chloroplastic beta-amylase 1 targeted by miR166m-5p and calcium-dependent protein kinase 1 targeted by miR166m-3p), suggesting arm-switching of miR166m play roles in soybean in response to salinity stress. Furthermore, two pairs of miRNA:lncRNA interacting partners (miR166i-5p and lncRNA Gmax_MSTRG.35921.1; and miR394a-3p and lncRNA Gmax_MSTRG.18616.1) were also discovered in reaction to salinity stress.
This study demonstrates how ncRNA involves in salinity stress responses in soybean by miRNA arm switching and miRNA:lncRNA interactions. The behaviors of ncRNAs revealed in this study will shed new light on molecular regulatory mechanisms of stress responses in plants, and hence provide potential new strategies for crop improvement.
Soybean is an economically and environmentally important crop, providing a major source of dietary protein and oil for human food and animal feed. Cultivated soybean was domesticated from its wild relatives approximately 5000 years ago in China before being introduced to other parts of Asia 2000 years ago [1, 2]. While Asian countries such as China are major consumers of soybean, it has now become an important cash crop in North and South America for export. Unfortunately, soybean is generally sensitive to saline soil, (often a result of long-term cultivation and over-fertilization of the land), resulting in the reduction of its yield . Earlier studies on salt tolerance of soybean have mainly focused on identifying the protein-coding genes and related mechanisms involved .
Besides protein-coding genes, noncoding RNAs (ncRNAs) are equally important regulators of gene expressions in plants, ranging from moderating development to coping with abiotic and biotic stresses. One type of ncRNAs, microRNAs (miRNAs) are typically 21-23 nucleotides long and can bind to the 3′ untranslated regions (3’UTRs) of messenger RNAs (mRNAs), resulting in the downregulation of their target genes. Soybean mutants lacking the miRNA biogenesis component, Dicer-like 1, have reduced seed size and defective seedling development . Many other miRNAs were also found to be associated with the development of soybean seeds [6,7,8], roots , and flowers . For instance, the overexpression of miR156b negatively regulated squamosa promoter-binding protein-like genes (GmSPLs) and postponed flowering time [11, 12]. In addition to the developmental roles, miRNAs are also associated with the responses toward abiotic stresses in soybean [3, 13,14,15,16,17], such as the involvements of gma-miR169c and gma-miR394a in drought tolerance and salt sensitivity [18,19,20], and miR399a and miR172c in root development under salinity stress [21, 22].
Despite a better understanding of the different roles of miRNAs in plants including soybean, several aspects remain poorly known, one of which is the phenomenon of miRNA arm-switching ([23, 24]; Fig. 1). In the biogenesis of miRNA, the dominant strand could change depending on cellular contents. Since the complementary strands of mature miRNAs will target different sets of mRNAs, alternative strand selection has been implicated in various biological mechanisms in animals (e.g. [26, 28, 29]). In plants, there are only a handful of studies related to miRNA arm-switching. For instance, miR393 and miR399 could undergo arm-switching in Arabidopsis under pathogen infection and phosphate deprivation [30, 31]. The analyses of 29 rice small RNA libraries also revealed miRNA arm-switching in different tissues and under abiotic stresses . Nevertheless, whether miRNA arm-switching could potentially contribute to soybean abiotic stress responses remains unexplored.
Another aspect that requires further investigation is the interactions between miRNAs and other ncRNAs, such as long noncoding RNAs (lncRNAs). Various roles of lncRNAs in gene regulation have already been demonstrated in plants, including transcriptional regulation [33,34,35,36] and alternative splicing . LncRNAs could also augment the repertoires of RNA molecules . For instance, the lncRNA IPS1 could interact with the miRNA miR399 to play a role in phosphate homeostasis by regulating its target gene PHO2 in Arabidopsis [38, 39] and Medicago ; and the lncRNA osa-eTM160 could interact with osa-miRNA160 to regulate rice development . In soybean, there has only been a genome-wide prediction of lncRNAs and miRNAs  but no systematic study on their interactions has been carried out yet. Therefore, the interactions between miRNAs and lncRNAs in relation to physiological responses toward abiotic stresses in soybean have remained elusive.
Here, we generated miRNA sequencing data from soybean leaf and root tissues from two soybean accessions (a cultivated Gylcine max C08 and a wild Glycine soja W05 accession as representatives) treated with NaCl at six time points with three biological replicates, in order to shed light on (i) differential miRNA expression, (ii) miRNA arm-switching, (iii) interactions between miRNAs and lncRNAs; and (iv) conserved and accession-specific responses between cultivated and wild accessions in response to salinity stress in soybean.
Differential miRNA expression in soybean under salinity stress
The setting of the experiment is described in materials and methods and a total of 72 small RNA datasets were generated for analysis (Fig. 2). Details of the datasets are listed in Supplementary Table 2. The small RNA sequencing data generated above were then used to identify miRNA candidates responsive to salinity stress. By comparing samples across different treatments, 35 upregulation and 77 downregulation events were found in the leaf and root between the different time points (Fig. 3).
Among the upregulated miRNA candidates (Fig. 3A), certain members within the same miRNA family sub-functionalized and were differentially expressed at different salinity conditions, including miR166 and miR390. In the case of the miR166 family, the expression of miR166m was found to be upregulated in the root of both W05 and C08 from 0 to 1 h upon salinity treatment, while the expression of miR166i-5p was upregulated in the leaf of both W05 and C08 from 4 h to 24 h. Similar sub-functionalization of expression patterns was also observed in the downregulated miRNAs, for instance, the miR156 family (Fig. 3B).
In addition, divergent expression trends of miRNAs upon salinity stress were also revealed within or between the tissues of the cultivated and wild soybean, including miR156b, miR156f, miR160a, miR166i, miR390a, miR390e, miR390f, miR390g, miR394a, miR4413a, miR4416c, and miR5225 (Fig. 3). For example, the expression of miR156b was upregulated in the root of C08 from 4 h to 24 h, while its expression was downregulated in the leaf of W05 from 4 h to 24 h.
On the other hand, some other miRNAs demonstrated consistent expression responses under salinity stress in the same tissue of both accessions, including miR156r, miR166i, miR3522, miR394a, miR397a, miR397b, miR5225, miR5671a, miR5761a, miR5761b, miR5770a, and miR5786 (Fig. 3). For example, miR156r was downregulated from 4 h to 24 h in the root of both C08 and W05. These miRNAs could play potentially important roles in the response to salinity stress and therefore deserve further investigations (Fig. 3).
MiRNA arm expression and switching of dominance
In addition to the conventional comparisons of differential expressions in miRNA family members upon salinity stress, we also investigated the expression patterns between the two arms (5p and 3p) generated from the same miRNA locus. Conserved expression patterns between miRNA arms were observed in both soybean accessions, an example of which is miR169l, where both the 5p and 3p arms had increased expression under salinity treatment between 4 to 24 h in the root samples (Fig. 3A).
However, divergent expression patterns between the two miRNA arms were also observed in the different tissues of the cultivated and wild soybean. For instance, decreased expression of the 3p arm was revealed in miR408a in the leaf tissues of W05 between 4 to 24 h, while downregulation of the 5p arm was observed between 4 to 24 h in the root tissues of C08 (Fig. 3B). Another example is miR408c which represents the same expression pattern as miR408a (Fig. 3B). Next, we explored whether miRNA arm-switching has occurred, using the stringent criteria as described in Materials and Methods, and found that miR166m underwent arm-switching in both soybean accessions (Fig. 4). Specifically, in all three biological replicates of the C08 root samples challenged with NaCl for 48 h, 5p arm dominance was recorded for miR166m (Fig. 4A; Supplementary Fig. 1), and yet at the same time point in the C08 leaf tissue, all three biological replicates exhibited 3p arm dominance for the same miRNA (Fig. 4B; Supplementary Fig. 1). The same pattern was also observed in W05 at 4 h upon salinity treatment, where miR166m exhibited different arm dominance in the root versus the leaf (Fig. 4C & D; Supplementary Fig. 1). In silico prediction of target genes of miR166m-5p and miR166m-3p in C08 and W05 were carried out and summarised in Supplementary Table 5. Conserved targets of the two arms were identified, such as the chloroplastic beta-amylase 1 targeted by miR166m-5p and calcium-dependent protein kinase 1 targeted by miR166m-3p. Predicted target genes regulated by miR-166 m and were previously reported relating to salinity stress were further selected for validation by dual luciferase reporter assay (Supplementary Fig. 3), suggesting arm-switching of miR166m play roles in soybean in response to salinity stress. In any case, this is the first time a robust and confident case of miRNA arm-switching was demonstrated in both cultivated and wild soybean accessions.
MiRNA and lncRNA interactions in soybean under salinity conditions
The potential interaction between miRNAs and lncRNAs was then investigated by combining the data generated in this study with our previous studies on soybean transcriptome reprograming in response to salinity stress [43,44,45]; Supplementary Table 4). Two criteria were used for predicting interactive pairs of miRNA and lncRNA, including: (1) sequence complementarity between the miRNA and the lncRNA, and (2) opposite trends in the expression patterns between the miRNA and the lncRNA upon salinity stress (Fig. 5A & B). Predicted miRNA:lncRNA pairs meeting these two criteria could suggest their functional relevance in salinity stress responses, and a total of four possible miRNA:lncRNA pairs were identified, with miR166i and miR394a each targeting two lncRNAs (Fig. 5C).
Specifically, in C08, miR166i-5p decreased in expression in leaf between 4 and 24 h upon salt treatment, while its predicted target lncRNAs, Gmax_MSTRG.11852.1 and Gmax_MSTRG.35921.1, increased in expression levels (Fig. 5C). On the other hand, in the root, upon salinity treatment between 2 and 4 h, miR166i-5p increased in expression level while the expressions of Gmax_MSTRG.11852.1 and Gmax_MSTRG.35921.1 decreased (Fig. 5C). Using the LAMP assay followed by RT-PCR, the interactions between miR166i and lncRNA Gmax_MSTRG.35921.1 was validated (Fig. 5D). This miRNA:lncRNA pair showed tissue-specific variations in their antagonistic expression patterns during soybean salinity stress.
On the other hand, time-specific variations in expression patterns were observed for miR394a and its predicted target lncRNAs in the root samples of C08, where different antagonistic expression patterns were observed between miR394a and its two predicted targets, Gmax_MSTRG.12568.1 and Gmax_MSTRG.18616.1, at different time points after salt treatment (Fig. 5C). However, only the interaction between miR394a and Gmax_MSTRG.18616.1 could be validated by the LAMP assay and RT-PCR (Fig. 5D). Furthermore, utilising the dual luciferase reporter assay, we validated the potential interactions of these two microRNAs with lncRNAs and genes related to salinity stress, suggesting the potential regulation of these two miRNAs under salinity stress (Supplementary Table 6 and Supplementary Fig. 3).
Previous studies have explored and established the potential roles of soybean miRNAs under salt treatment, ranging from elucidating individual miRNA functions [19, 21] to the genome-wide profiling of miRNAs via transcriptome sequencing [13, 48, 49]. This study further advances the knowledge of miRNAs in several aspects.
In terms of miRNA sequencing data, previous studies have usually focused on a single tissue/time point/accession. For instance, information was provided on a salt-sensitive soybean inbred line “HJ-1” based on a single tissue type at a single time point, i.e. root at 48 h ; effects of salt treatment on the leaf and root samples of “William 82” at various early time points (0, 1, 3, 6, 9 and 12 h; ); and longer-term effects on the root tissues of “William 82” at 15d after salt treatment . On the other hand, the data generated in this study covers two tissue types (leaf and root) at a large range of time points (0, 1, 2, 4, 24 and 48 h) from two accessions (G. soja “W05” and G. max “C08”), providing an unprecedented opportunity to understand the various dynamic contributions of miRNAs to soybean physiology in coping with salinity.
In addition to providing a list of differentially expressed miRNAs, this study also uncovered new aspects of miRNA responses to abiotic stresses in plants. Similar to other studies, we found that in this study, members of the same miRNA family exhibited different expression patterns upon salinity challenge, but because we could compare the data between two accessions in the same tissues and time points, we were then able to reveal the conserved and divergent trends of expression patterns of miRNAs, which in turn allowed us to narrow down those miRNAs which are potentially being negatively selected in the salinity responses in soybean, which can then lead us to their target genes that could be the important players in these responses.
This study has also shed new light on how the different arms of the miRNA could contribute to plant physiology. Here we have identified conserved and divergent expression patterns between the two arms across the two accessions. Specifically, the two arms of miR169l exhibited the same expression profiles/response upon salinity challenge. This could potentially be used as another strategy to identify miRNAs which could play important roles in soybean physiology, by screening for those miRNAs with similar expression profiles under abiotic stresses.
Furthermore, in animal models, miRNA arm-switching is now well-known as a mechanism for controlling various biological processes, but there are still only limited studies on Arabidopsis and rice to understand how such miRNA arm selection could be a potential stress response strategy in plants [30,31,32]. This study also provided the first documented case of a miRNA (miR166m) undergoing arm-switching in different tissues in soybean, in both G. max and G. soja, when under salinity stress. Last but not least, we have also revealed for the first time miRNA and lncRNA interactions in soybean under salinity stress (Fig. 6). Previous studies in other plant species have revealed potential miRNA:lncRNA interactions under different conditions, such as phosphate starvation, pathogen infection, and heat stress [47, 50, 51]. Two pairs of miRNA:lncRNA interacting partners (miR166i-5p and lncRNAs Gmax_MSTRG.35921.1; and miR394a-3p and lncRNA Gmax_MSTRG.18616.1) could be identified and validated in soybean under salinity stress. Intriguingly, the interactions between miR166i-5p and its lncRNA interacting partner showed tissue-specific expression patterns, whereas that between miR394a and its lncRNA partner exhibited time-dependent patterns after salt treatment. It has previously been proposed that down-regulated miRNAs may target genes involved in stress responses, while upregulated miRNAs probably target genes involved in development . Validation of the predicted targets for these selected miRNAs indicates the potential regulation of these miRNAs in response to salinity stress given these targets such as SOS2-like protein kinase , calcium-dependent protein kinase 1 , nudix hydrolase 2  were previously shown to be related to salinity stress. Nonetheless, whether and how these different kinds of noncoding RNAs interaction under salinity stress spatially/temporally, and whether these interactions could contribute to the improvement on salinity tolerance of soybean, remain to be tested.
This study profiled the time-series differential expression patterns of a group of miRNAs in the root and leaf of two soybean accessions under salinity stress. We also discovered the first evidence of miRNA arm-switching in soybean, indicating that different tissues may have different miRNA arm preference that contributes to different biological functions. Two validated pairs of miRNA:lncRNA interacting partners involved in the response to salinity treatment in soybean were also identified, setting up the foundation for further investigations into the roles of different types of non-coding RNAs in plant physiology.
Soybean treatments and sample collection
Cultivated (C08) and wild (W05) soybean plants were grown and harvested as described in one of our previous publications . Salt treatment was carried out by transferring to half-strength Hoagland solution supplemented with 0.9% (w/v; ~ 150 mM) NaCl when the primary leaves were fully opened. Three independent sets of primary leaf and root samples of each cultivar and treatment were harvested at 0, 1, 2, 4, 24, and 48 h post-treatment. A total of 72 samples were generated.
Small RNA sequencing
MicroRNA (miRNA) was isolated using mirVana™ miRNA Isolation Kit (Thermo Scientific) with Acid-Phenol:Chloroform (Thermo Scientific), following the manufacturer’s protocol. After the quality of miRNA was examined via bio-analysis, small RNA library construction and deep sequencing were performed at BGI (Shenzhen, Guangdong, China). Briefly, after PAGE purification, small RNA molecules between 16 and 30 bases were collected and ligated with a pair of Solexa adaptors to their 5′ and 3′ ends respectively. The small RNAs were then amplified using the adaptor primers for 17 cycles. Amplicons of around 90 bp (small RNA + adaptors) were isolated from the agarose gel and were used directly for cluster generation and sequencing analyses using Illumina Hiseq 2000 according to the manufacturer’s instructions.
MiRNA arm-switching detection and differential miRNA expression
FastQC was run as quality control  for the 72 small RNA datasets. Adaptor sequences were trimmed from small RNA sequencing reads, and reads with the Phred quality score less than 20 were removed. The expression levels of the 5p and 3p arms of soybean miRNAs from miRBase (Release 22.1) were calculated based on the number of sequencing reads mapped to the respective arm region in the miRNA hairpin using bowtie/miRDeep2 [57, 58]. The miRNAs having either arm with absolute counts > 50 were included in the arm-switching analysis. The formula, ω = 5p/(5p + 3p), where 5p and 3p refer to the number of predicted 5p arm and 3p arm, respectively, was adopted to calculate the arm selection value, ω, which ranged from 0 to 1. Smaller ω values indicate higher tendencies of 3p preference and larger values indicate higher tendencies of 5p preference. We adopted ω < 0.3 as the indicator of 3p dominance and ω > 0.7 as the indicator of 5p dominance. The shift of arm dominance is defined as arm switching . The expression raw read table, calculated as previously described, was submitted to Degust  for the comparative analyses of differential miRNA gene expression using the edgeR method . (Degust link: https://degust.erc.monash.edu/degust/compare.html?code=708756e3e3e0ad71f68186a9e83f19c2#/.)
Interaction prediction between miRNAs and long noncoding RNAs (lncRNAs)
The soybean lncRNA annotation was retrieved from a previously published soybean lncRNA catalog . The published time-series transcriptome data upon salt treatment [43, 45] were processed as described . Briefly, adapter and quality trimming were performed using Trimmomatic 0.36 . Then, the clean reads of C08 and W05 were mapped to the reference genomes of Williams 82 (G. max)  and W05 (G. soja) , respectively, using TopHat v2.1.1 . With the read mapping results, Cufflinks [64, 65] was used to perform annotation-free transcriptome assembly. The assembled transcripts were merged using the Cufflinks script, “cuffmerge”, to produce unified transcript sets separately in both genomes. These transcripts were then compared against the protein-coding transcripts, as well as small non-coding transcripts as predicted by RNAmmer v1.2 , tRNAscan-SE v2.0  and Infernal v1.1.2  with Rfam v13.0 database , and searched against the nonredundant protein database from NCBI using BLASTx  to identify potential lncRNAs. Trinity  was used for de novo transcriptome assembly with the clean reads, separately for the C08 and W05 data. All Trinity-assembled transcripts were aligned to the two reference genomes using GMAP v2018-03-25 . The above-mentioned potential lncRNAs were used for downstream analyses only when they were supported by both the Trinity-assembled transcripts and the lncRNA catalog. The differential gene expression analysis was performed using the Cufflinks script, “cuffdiff”, with the criteria set at |log2FC| > 1 and q-value < 0.05. Target prediction between miRNA and lncRNA was performed by RNAhybrid , based on the seed region pair-matching with default parameters, and by LncMirNet  based on deep-learning of ribonucleic acid sequences, respectively. For the comparison across the time-series of salinity treatment, parameters were set to screen for differentially expressed miRNAs and lncRNAs to investigate potential interactions between them. For lncRNAs, |log2FC| ≥ 1 and q-value ≤0.05 were used as screening criteria for differential expression. For miRNAs, counts per million (CPM) ≥ 50 in at least one sample, |log2FC| ≥ 1, and false discovery rate (FDR) cut-off ≤0.05 are the criteria used for differential expression screening.
Target gene prediction of miRNAs
Target predictions between miRNAs and the 5′ untranslated regions (UTR), coding sequence (CDS) and 3’UTR of protein-coding genes were performed by psRNATarget . The functional term annotations were performed using eggNOG  with default parameters and taxon restricted to Fabids (Taxon ID:91835). Genes were assigned with Gene Ontology (GO), EuKaryotic Orthologous Groups (KOG), Kyoto Encyclopedia of Genes and Genomes (KEGG), and KEGG Orthology (KO) terms. Functional enrich of mirna target gene was tested using function ‘compareCluster()’ in R package ‘clusterProfiler’ v.3.16.1  under the environment of R 4.0.4 . Significantly enriched terms were determined with pvalueCutoff = 0.05, pAdjustMethod = “BH”, and qvalueCutoff = 0.2. Data was visualised using R packages ‘ggplot2’ .
Labeled miRNA pull-down (LAMP) assay
The predicted interactions between miRNAs and lncRNAs were validated by biotin-labeled miRNA pull-down (LAMP) assay  using Dynabeads®M-280 Streptavidin (Invitrogen™) and biotin-labeled gma-miR166i-5p and gma-miR394a-3p (Integrated DNA Technologies [IDT]), followed by RT-PCR. The primers of four target lncRNAs for PCR are listed in Supplementary Table 1. RNA from the leaf of C08 under salt treatment for 24 h was the starting material and the RNA pulled-down products were reverse-transcribed using the iScript cDNA Synthesis Kit (BIO-RAD) following the manufacturer’s protocol to obtain the cDNAs. Those originating from the four target lncRNAs were amplified by PCR and separated on agarose gel. The target bands were excised and purified by QIAquick gel extraction kit (Qiagen) following manufacturer’s protocol and sent for Sanger sequencing for confirmation.
Dual luciferase reporter (DLR) assay
Genomic DNA and cDNA were used to amplify microRNA hairpins with flanking sequences, target genes and lncRNAs. Primers information were listed in Supplementary Table 1. Sequences used in DLR assay were listed in Supplementary Table 7. MicroRNAs were cloned into pAC5.1 vector, while the target genes and lncRNAs were cloned into psicheck-2 vector. Dual luciferase reporter assay was conducted in Drosophila S2 cells as previously described  with following modifications: Drosophila S2 cells were cultured in Shields and Sang M3 Insect Medium (Sigma) supplemented with 10% (v/v) heat-inactivated fetal bovine serum (Gibco, Life Technologies) and 1: 100 penicillin-streptomycin (Gibco, Life Technologies). Measurement was performed at 36-48 h post-transfection by Tecan Spark 10 M Microplate Reader (Eastwin International Trading Ltd) with two technical replicates four biological replicates.
Availability of data and materials
The raw reads of the 72 small RNA samples generated in this study have been deposited to the NCBI database under the BioProject accessions: PRJNA720229.
Li Y, Guan R, Liu Z, Ma Y, Wang L, Li L, et al. Genetic structure and diversity of cultivated soybean (Glycine max (L.) Merr.) landraces in China. Theor Appl Genet. 2008;117:857–71. https://doi.org/10.1007/s00122-008-0825-0.
Wilson RF. Soybean: market driven research needs. In: Stacey G, editor. Genetics and genomics of soybean. Plant genetics and genomics: crops and models. New York: Springer; 2008. p. 2. https://doi.org/10.1007/978-0-387-72299-3_1.
Dong Z, Shi L, Wang Y, Chen L, Cai Z, Wang Y, et al. Identification and dynamic regulation of microRNAs involved in salt stress responses in functional soybean nodules by high-throughput sequencing. Int J Mol Sci. 2013;14:2717–38. https://doi.org/10.3390/ijms14022717.
Phang T-H, Shao G, Lam H-M. Salt tolerance in soybean. J Integr Plant Biol. 2008;50:1196–212. https://doi.org/10.1111/j.1744-7909.2008.00760.x.
Curtin SJ, Michno J-M, Campbell BW, Gil-Humanes J, Mathioni SM, Hammond R, et al. MicroRNA maturation and MicroRNA target gene expression regulation are severely disrupted in soybean dicer-like1 double mutants. G3 (Bethesda). 2015;6:423–33. https://doi.org/10.1534/g3.115.022137.
Gupta OP, Dahuja A, Sachdev A, Kumari S, Jain PK, Vinutha T, et al. Conserved miRNAs modulate the expression of potential transcription factors of isoflavonoid biosynthetic pathway in soybean seeds. Mol Biol Rep. 2019;46:3713–30. https://doi.org/10.1007/s11033-019-04814-7.
Wang T, Sun M-Y, Wang X-S, Li W-B, Li Y-G. Over-expression of GmGIa-regulated soybean miR172a confers early flowering in transgenic Arabidopsis thaliana. Int J Mol Sci. 2016;17:645. https://doi.org/10.3390/ijms17050645.
Yu J-Y, Zhang Z-G, Huang S-Y, Han X, Wang X-Y, Pan W-J, et al. Analysis of miRNAs targeted storage regulatory genes during soybean seed development based on transcriptome sequencing. Genes (Basel). 2019;10:408. https://doi.org/10.3390/genes10060408.
Noon JB, Hewezi T, Baum TJ. Homeostasis in the soybean miRNA396– GRF network is essential for productive soybean cyst nematode infections. J Exp Bot. 2019;70:1653–68. https://doi.org/10.1093/jxb/erz022.
Kulcheski FR, Molina LG, da Fonseca GC, de Morais GL, de Oliveira LFV, Margis R. Novel and conserved microRNAs in soybean floral whorls. Gene. 2016;575(2 Pt 1):213–23. https://doi.org/10.1016/j.gene.2015.08.061.
Cao D, Li Y, Wang J, Nan H, Wang Y, Lu S, et al. GmmiR156b overexpression delays flowering time in soybean. Plant Mol Biol. 2015;89:353–63. https://doi.org/10.1007/s11103-015-0371-5.
Ding X, Zhang H, Ruan H, Li Y, Chen L, Wang T, et al. Exploration of miRNA-mediated fertility regulation network of cytoplasmic male sterility during flower bud development in soybean. 3 Biotech. 2019;9:22. https://doi.org/10.1007/s13205-018-1543-1.
Li H, Dong Y, Yin H, Wang N, Yang J, Liu X, et al. Characterization of the stress associated microRNAs in Glycine max by deep sequencing. BMC Plant Biol. 2011;11:170. https://doi.org/10.1186/1471-2229-11-170.
Xu F, Liu Q, Chen L, Kuang J, Walk T, Wang J, et al. Genome-wide identification of soybean microRNAs and their targets reveals their organ-specificity and responses to phosphate starvation. BMC Genomics. 2013;14:66. https://doi.org/10.1186/1471-2164-14-66.
Zeng HQ, Zhu YY, Huang SQ, Yang ZM. Analysis of phosphorus-deficient responsive miRNAs and cis-elements from soybean (Glycine max L.). J Plant Physiol. 2010;167:1289–97. https://doi.org/10.1016/j.jplph.2010.04.017.
Zhang S, Wang Y, Li K, Zou Y, Chen L, Li X. Identification of cold-responsive miRNAs and their target genes in nitrogen-fixing nodules of soybean. Int J Mol Sci. 2014;15:13596–614. https://doi.org/10.3390/ijms150813596.
Zheng Y, Hivrale V, Zhang X, Valliyodan B, Lelandais-Brière C, Farmer AD, et al. Small RNA profiles in soybean primary root tips under water deficit. BMC Syst Biol. 2016;10:126. https://doi.org/10.1186/s12918-016-0374-0.
Ni Z, Hu Z, Jiang Q, Zhang H. Overexpression of gma-MIR394a confers tolerance to drought in transgenic Arabidopsis thaliana. Biochem Biophys Res Commun. 2012;427:330–5. https://doi.org/10.1016/j.bbrc.2012.09.055.
Ni Z, Hu Z, Jiang Q, Zhang H. GmNFYA3, a target gene of miR169, is a positive regulator of plant tolerance to drought stress. Plant Mol Biol. 2013;82:113–29. https://doi.org/10.1007/s11103-013-0040-5.
Yu Y, Ni Z, Wang Y, Wan H, Hu Z, Jiang Q, et al. Overexpression of soybean miR169c confers increased drought stress sensitivity in transgenic Arabidopsis thaliana. Plant Sci. 2019;285:68–78. https://doi.org/10.1016/j.plantsci.2019.05.003.
Sahito ZA, Wang L, Sun Z, Yan Q, Zhang X, Jiang Q, et al. The miR172c-NNC1 module modulates root plastic development in response to salt in soybean. BMC Plant Biol. 2017;17:229. https://doi.org/10.1186/s12870-017-1161-9.
Sun Z, Wang Y, Mou F, Tian Y, Chen L, Zhang S, et al. Genome-wide small RNA analysis of soybean reveals auxin-responsive microRNAs that are differentially expressed in response to salt stress in root apex. Front Plant Sci. 2016;6. https://doi.org/10.3389/fpls.2015.01273.
Griffiths-Jones S, Hui JHL, Marco A, Ronshaugen M. MicroRNA evolution by arm switching. EMBO Rep. 2011;12:172–7. https://doi.org/10.1038/embor.2010.191.
Marco A, Hui JHL, Ronshaugen M, Griffiths-Jones S. Functional shifts in insect microRNA evolution. Genome Biol Evol. 2010;2:686–96. https://doi.org/10.1093/gbe/evq053.
Li C, Wong A, Wang S, Jia Q, Chuang W-P, Bendena W, et al. miRNA-mediated interactions in and between plants and insects. Int J Mol Sci. 2018;19:3239. https://doi.org/10.3390/ijms19103239.
Medley JC, Panzade G, Zinovyeva AY. microRNA strand selection: unwinding the rules. WIREs RNA. 2020. https://doi.org/10.1002/wrna.1627.
Budak H, Kaya SB, Cagirici HB. Long non-coding RNA in plants in the era of reference sequences. Front Plant Sci. 2020;11:276. https://doi.org/10.3389/fpls.2020.00276.
Hui JHL, Marco A, Hunt S, Melling J, Griffiths-Jones S, Ronshaugen M. Structure, evolution and function of the bi-directionally transcribed iab-4/iab-8 microRNA locus in arthropods. Nucleic Acids Res. 2013;41:3352–61. https://doi.org/10.1093/nar/gks1445.
Nong W, Cao J, Li Y, Qu Z, Sun J, Swale T, et al. Jellyfish genomes reveal distinct homeobox gene clusters and conservation of small RNA processing. Nat Commun. 2020;11:3051. https://doi.org/10.1038/s41467-020-16801-9.
Hsieh L-C, Lin S-I, Shih AC-C, Chen J-W, Lin W-Y, Tseng C-Y, et al. Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing. Plant Physiol. 2009;151:2120–32. https://doi.org/10.1104/pp.109.147280.
Zhang X, Zhao H, Gao S, Wang W-C, Katiyar-Agarwal S, Huang H-D, et al. Arabidopsis Argonaute 2 regulates innate immunity via miRNA393(∗)-mediated silencing of a Golgi-localized SNARE gene, MEMB12. Mol Cell. 2011;42:356–66. https://doi.org/10.1016/j.molcel.2011.04.010.
Hu W, Wang T, Yue E, Zheng S, Xu J-H. Flexible microRNA arm selection in rice. Biochem Biophys Res Commun. 2014;447:526–30. https://doi.org/10.1016/j.bbrc.2014.04.036.
Ariel F, Jegu T, Latrasse D, Romero-Barrios N, Christ A, Benhamed M, et al. Noncoding transcription by alternative RNA polymerases dynamically regulates an auxin-driven chromatin loop. Mol Cell. 2014;55:383–96. https://doi.org/10.1016/j.molcel.2014.06.011.
Ariel F, Lucero L, Christ A, Mammarella MF, Jegu T, Veluchamy A, et al. R-loop mediated trans action of the APOLO long noncoding RNA. Mol Cell. 2020;77:1055–1065.e4. https://doi.org/10.1016/j.molcel.2019.12.015.
Csorba T, Questa JI, Sun Q, Dean C. Antisense COOLAIR mediates the coordinated switching of chromatin states at FLC during vernalization. Proc Natl Acad Sci. 2014;111:16160–5. https://doi.org/10.1073/pnas.1419030111.
Heo JB, Sung S. Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science (80- ). 2011;331:76–9. https://doi.org/10.1126/science.1197349.
Bardou F, Ariel F, Simpson CG, Romero-Barrios N, Laporte P, Balzergue S, et al. Long noncoding RNA modulates alternative splicing regulators in Arabidopsis. Dev Cell. 2014;30:166–76. https://doi.org/10.1016/j.devcel.2014.06.017.
Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, et al. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39(8):1033–7.
Lin S-I, Chiang S-F, Lin W-Y, Chen J-W, Tseng C-Y, Wu P-C, et al. Regulatory network of MicroRNA399 and PHO2 by systemic signaling. Plant Physiol. 2008;147:732–46. https://doi.org/10.1104/pp.108.116269.
Wang T, Zhao M, Zhang X, Liu M, Yang C, Chen Y, et al. Novel phosphate deficiency-responsive long non-coding RNAs in the legume model plant Medicago truncatula. J Exp Bot. 2017;68:5937–48. https://doi.org/10.1093/jxb/erx384.
Wang M, Wu H-J, Fang J, Chu C, Wang X-J. A long noncoding RNA involved in rice reproductive development by negatively regulating osa-miR160. Sci Bull. 2017;62:470–5. https://doi.org/10.1016/j.scib.2017.03.013.
Ye C-Y, Xu H, Shen E, Liu Y, Wang Y, Shen Y, et al. Genome-wide identification of non-coding RNAs interacted with microRNAs in soybean. Front Plant Sci. 2014;5. https://doi.org/10.3389/fpls.2014.00743.
Liu A, Xiao Z, Li M-W, Wong F-L, Yung W-S, Ku Y-S, et al. Transcriptomic reprogramming in soybean seedlings under salt stress. Plant Cell Environ. 2019;42:98–114. https://doi.org/10.1111/pce.13186.
Lin X, Lin W, Ku Y-S, Wong F-L, Li M-W, Lam H-M, et al. Analysis of soybean long non-coding RNAs reveals a subset of small peptide-coding transcripts. Plant Physiol. 2020;182:1359–74. https://doi.org/10.1104/pp.19.01324.
Xie M, Chung CY-L, Li M-W, Wong F-L, Wang X, Liu A, et al. A reference-grade wild soybean genome. Nat Commun. 2019;10:1216. https://doi.org/10.1038/s41467-019-09142-9.
Fernandes J, Acuña S, Aoki J, Floeter-Winter L, Muxel S. Long non-coding RNAs in the regulation of gene expression: physiology and disease. Non-Coding RNA. 2019;5:17. https://doi.org/10.3390/ncrna5010017.
Zhang X, Wang W, Zhu W, Dong J, Cheng Y, Yin Z, et al. Mechanisms and functions of long non-coding RNAs at multiple regulatory levels. Int J Mol Sci. 2019;20:5573. https://doi.org/10.3390/ijms20225573.
Liu W, Deng Y, Zhou Y, Chen H, Dong Y, Wang N, et al. Normalization for relative quantification of mRNA and microRNA in soybean exposed to various abiotic stresses. PLoS One. 2016;11:e0155606. https://doi.org/10.1371/journal.pone.0155606.
Wang Q, Yang Y, Lu G, Sun X, Feng Y, Yan S, et al. Genome-wide identification of microRNAs and phased siRNAs in soybean roots under long-term salt stress. Genes Genomics. 2020;42:1239–49. https://doi.org/10.1007/s13258-020-00990-0.
Fan C, Hao Z, Yan J, Li G. Genome-wide identification and functional analysis of lincRNAs acting as miRNA targets or decoys in maize. BMC Genomics. 2015;16:793. https://doi.org/10.1186/s12864-015-2024-0.
He X, Guo S, Wang Y, Wang L, Shu S, Sun J. Systematic identification and analysis of heat-stress-responsive lncRNAs, circRNAs and miRNAs with associated co-expression and ceRNA networks in cucumber ( Cucumis sativus L.). Physiol Plant. 2020;168:736–54. https://doi.org/10.1111/ppl.12997.
Sanz-Carbonell A, Marques MC, Bustamante A, Fares MA, Rodrigo G, Gomez G. Inferring the regulatory network of the miRNA-mediated response to biotic and abiotic stress in melon. BMC Plant Biol. 2019;19:78. https://doi.org/10.1186/s12870-019-1679-0.
Liu J, Ishitani M, Halfter U, Kim CS, Zhu JK. The Arabidopsis thaliana SOS2 gene encodes a protein kinase that is required for salt tolerance. Proc Natl Acad Sci U S A. 2000;97(7):3730–4. https://doi.org/10.1073/pnas.060034197.
Vivek PJ, Resmi MS, Sreekumar S, Sivakumar KC, Tuteja N, Soniya EV. Calcium-dependent protein kinase in ginger binds with importin-α through its junction domain for nuclear localization, and further interacts with NAC transcription factor. Front Plant Sci. 2017;7:1909. https://doi.org/10.3389/fpls.2016.01909.
Huang H, Cao H, Niu Y, et al. Expression analysis of Nudix hydrolase genes in Chrysanthemum lavandulifolium. Plant Mol Biol Rep. 2012;30:973–82. https://doi.org/10.1007/s11105-011-0401-7.
Andrews S. FastQC - a quality control tool for high throughput sequence data: Babraham Bioinformatics; 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.
Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.
Powell D. Degust: visualize, explore and appreciate RNA-seq differential gene-expression data. In: COMBINE RNA-seq workshop; 2015.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40. https://doi.org/10.1093/bioinformatics/btp616.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20. https://doi.org/10.1093/bioinformatics/btu170.
Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, et al. Genome sequence of the palaeopolyploid soybean. Nature. 2010;463:178–83. https://doi.org/10.1038/nature08670.
Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25:1105–11. https://doi.org/10.1093/bioinformatics/btp120.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7:562–78. https://doi.org/10.1038/nprot.2012.016.
Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol. 2013;31:46–53. https://doi.org/10.1038/nbt.2450.
Lagesen K, Hallin P, Rødland EA, Stærfeldt H-H, Rognes T, Ussery DW. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 2007;35:3100–8. https://doi.org/10.1093/nar/gkm160.
Schattner P, Brooks AN, Lowe TM. The tRNAscan-SE, snoscan and snoGPS web servers for the detection of tRNAs and snoRNAs. Nucleic Acids Res. 2005;33(Web Server):W686–9. https://doi.org/10.1093/nar/gki366.
Nawrocki EP, Eddy SR. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics. 2013;29:2933–5. https://doi.org/10.1093/bioinformatics/btt509.
Kalvari I, Argasinska J, Quinones-Olvera N, Nawrocki EP, Rivas E, Eddy SR, et al. Rfam 13.0: shifting to a genome-centric resource for non-coding RNA families. Nucleic Acids Res. 2018;46(D1):D335–42.
Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52. https://doi.org/10.1038/nbt.1883.
Wu TD, Watanabe CK. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics. 2005;21:1859–75. https://doi.org/10.1093/bioinformatics/bti310.
Krüger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 2006;(34 Web Server):W451–4. https://doi.org/10.1093/nar/gkl243.
Yang S, Wang Y, Lin Y, Shao D, He K, Huang L. LncMirNet: predicting LncRNA–miRNA interaction based on deep learning of ribonucleic acid sequences. Molecules. 2020;25:4372. https://doi.org/10.3390/molecules25194372.
Dai X, Zhuang Z, Zhao PX. psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic Acids Res. 2018;46(W1):W49–54. https://doi.org/10.1093/nar/gky316.
Huerta-Cepas J, Szklarczyk D, Heller D, Hernández-Plaza A, Forslund SK, Cook H, et al. eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 2018;47:309–14. https://doi.org/10.1093/nar/gky1085.
Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7. https://doi.org/10.1089/omi.2011.0118.
R Core Team. R: a language and environment for statistical computing. 2021. https://www.r-project.org/.
Wickham H. ggplot2: elegant graphics for data analysis: Springer-Verlag; 2016. https://ggplot2.tidyverse.org
Hsu R-J, Yang H-J, Tsai H-J. Labeled microRNA pull-down assay system: an experimental approach for high-throughput identification of microRNA-target mRNAs. Nucleic Acids Res. 2009;37:e77. https://doi.org/10.1093/nar/gkp274.
Qu Z, Bendena WG, Nong W, Siggens KW, Noriega FG, Kai ZP, et al. MicroRNAs regulate the sesquiterpenoid hormonal pathway in Drosophila and other arthropods. Proc Biol Sci. 2017;284(1869):20171827. https://doi.org/10.1098/rspb.2017.1827.
Ms. Jee Yan Chu copy-edited this manuscript. Any opinions, findings, conclusions or recommendations expressed in this publication do not reflect the views of the Government of the Hong Kong Special Administrative Region or the Innovation and Technology Commission.
C.L. is supported by postgraduate studentship at The Chinese University of Hong Kong. This work was supported by Hong Kong Research Grant Council Area of Excellence scheme (AoE/M-403/16) awarded to H.-M.L., T.-F.C., and J.H.L.H.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Mapped reads of miR166m in individual samples. Supplementary Figure 2.1. PHRED Score Distribution. Supplementary Figure 2.2. Read Length Distribution. Supplementary Figure 2.3. Quality Control Statistics. Supplementary Figure 2.4. RNA Type. Supplementary Figure 2.5. miRNA Complexity. Supplementary Figure 2.6. Contamination. Supplementary Figure 3. Dual luciferase reporter assay.
Oligos and primer sequences used in this study. Supplementary Table 2. Number of raw reads generated in the datasets in this study. Supplementary Table 3. Differential miRNA expression in different samples as shown in Fig. 3. Supplementary Table 4. Differential lncRNA expression in different samples as shown in Fig. 5B and C. Supplementary Table 5. Target predictions of miR166m-5p and miR166m-3p in C08 and W05. Supplementary Table 6. Genes selected for validation. Supplementary Table 7. Sequences used in the DLR assay.
About this article
Cite this article
Li, C., Nong, W., Zhao, S. et al. Differential microRNA expression, microRNA arm switching, and microRNA:long noncoding RNA interaction in response to salinity stress in soybean. BMC Genomics 23, 65 (2022). https://doi.org/10.1186/s12864-022-08308-y
- microRNA (miRNA)
- Long noncoding RNA (lncRNA)
- microRNA arm switching
- miRNA:lncRNA interactions
- Salt stress