Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Transcriptome-wide identification of novel circular RNAs in soybean in response to low-phosphorus stress

  • Lingling Lv,

    Roles Formal analysis, Investigation, Methodology, Validation, Writing – original draft

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

  • Kaiye Yu,

    Roles Investigation, Methodology

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

  • Haiyan Lü,

    Roles Funding acquisition, Software, Writing – review & editing

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

  • Xiangqian Zhang,

    Roles Investigation, Software

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

  • Xiaoqian Liu,

    Roles Investigation, Visualization

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

  • Chongyuan Sun,

    Roles Validation

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

  • Huanqing Xu,

    Roles Visualization

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

  • Jinyu Zhang,

    Roles Writing – review & editing

    Affiliation Collaborative Innovation Center of Modern Biological Breeding, Henan Institute of Science and Technology, Xinxiang, China

  • Xiaohui He ,

    Roles Resources, Writing – review & editing

    zhangd@henau.edu.cn (DZ); hexh@zzu.edu.cn (XH)

    Affiliation Smart City Institute, Zhengzhou University, Zhengzhou, China

  • Dan Zhang

    Roles Conceptualization, Data curation, Funding acquisition, Supervision, Writing – review & editing

    zhangd@henau.edu.cn (DZ); hexh@zzu.edu.cn (XH)

    Affiliation Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China

Abstract

Low-phosphorus (LP) stress is a major factor limiting the growth and yield of soybean. Circular RNAs (circRNAs) are novel noncoding RNAs that play a crucial role in plant responses to abiotic stress. However, how LP stress mediates the biogenesis of circRNAs in soybean remains unclear. Here, to explore the response mechanisms of circRNAs to LP stress, the roots of two representative soybean genotypes with different P-use efficiency, Bogao (a LP-sensitive genotype) and Nannong 94156 (a LP-tolerant genotype), were used for the construction of RNA sequencing (RNA-seq) libraries and circRNA identification. In total, 371 novel circRNA candidates, including 120 significantly differentially expressed (DE) circRNAs, were identified across different P levels and genotypes. More DE circRNAs were significantly regulated by LP stress in Bogao than in NN94156, suggesting that the tolerant genotype was less affected by LP stress than the sensitive genotype was; in other words, NN94156 may have a better ability to maintain P homeostasis under LP stress. Moreover, a positive correlation was observed between the expression patterns of P stress-induced circRNAs and their circRNA-host genes. Gene Ontology (GO) enrichment analysis of these circRNA-host genes and microRNA (miRNA)-targeted genes indicated that these DE circRNAs were involved mainly in defense responses, ADP binding, nucleoside binding, organic substance catabolic processes, oxidoreductase activity, and signal transduction. Together, our results revealed that LP stress can significantly alter the genome-wide profiles of circRNAs and indicated that the regulation of circRNAs was both genotype and environment specific in response to LP stress. LP-induced circRNAs might provide a rich resource for LP-responsive circRNA candidates for future studies.

Introduction

Phosphorus (P), a macronutrient present in all living cells, is essential for plant growth and development because of its functions in energy transfer, in DNA in terms of inheritance, in membrane components, in carbon metabolism and in regulation and signaling [1]. The low availability of soil P and the acquisition of P exclusively through plant roots result in P being a major factor limiting plant growth and development. The sustainability of the use of P fertilizers to optimize crop yields is in jeopardy due to the crisis of the global reserve of rock phosphate and P pollution. Compared with other crop species, soybean needs more P for growth and development; thus, low-P (LP) stress has become a major factor limiting soybean production [2]. Therefore, elucidating the response and adaptation mechanism of soybean to LP stress is of great urgency and importance.

Plant growth and crop productivity are greatly affected by various environmental factors. As such, plants have evolved sophisticated adaptation mechanisms for sensing and responding to suboptimal environmental stresses. Recent studies have revealed that a vast amount of noncoding RNAs (ncRNAs), for example, small interfering RNA (siRNA), microRNA (miRNA), and long noncoding RNA (lncRNA), are induced during abiotic stress, including LP stress [3], suggesting their indispensable role in regulating plant defense responses. In Arabidopsis, miRNA can regulate phosphate homeostasis [4]. Changes in miRNA in leaves and roots constitute an important mechanism in maize to adapt to low-phosphate environments [5]. Phosphate starvation also induces miRNA responses in Populus tomentosa [6]. Moreover, as a long-distance signal, miR399 can regulate phosphate homeostasis in plants [7] and can serve as a potential integrator of photoresponses and phosphate homeostasis [8, 9]. In Arabidopsis, a new class of small RNAs has been identified in phosphate-starved roots [10], and these small RNA mediate responses to P deficiency and regulate P homeostasis [11]. In yeast, lncRNA can dynamically regulate pho1 expression by recruiting RNAi and the exosome on the basis of phosphate levels [12]. Moreover, lncRNA not only governs the expression of the phosphate transporter gene Pho84 but also exerts linkage effects on prt lncRNA and pho1 genes in flanking regions [13]. Recent studies have shown that lncRNAs exert a non-cell-autonomous response to P deficiency and may act as systemic signaling agents at the whole-plant level to coordinate early P deficiency signals [14]. LncRNAs also play key roles in regulating mRNA levels of a large number of genes associated with P starvation responses [15]. Moreover, in the legume model plant Medicago truncatula, novel P deficiency-induced lncRNAs can regulate P transport and P-deficiency signaling, suggesting that lncRNAs play an important role in regulating plant responses to LP stress [16].

Recently, a newly characterized type of endogenous ncRNA, circular RNA (circRNA), has been intensively studied in a variety of species [17]. Unlike traditional linear RNAs, circRNAs are generated from head-to-tail back-splicing; thus, they do not have the terminal 5′ cap and the 3′ polyadenylated tail structures found in linear genes [18]. Increasing amounts of evidence have indicated that circRNAs are commonly produced by thousands of genes and have substantial functions in regulating gene expression at multiple levels in eukaryotic cells [19]. For example, circRNAs can function as transcriptional and posttranscriptional regulators, miRNA sponges, small nucleolar RNAs (snoRNAs), RNA processing intermediates, and transcriptional cis-regulators [1921]. Moreover, increasing amounts of evidence have indicated that circRNAs play a key role in plant responses to abiotic stress [20, 2226]. For example, circRNAs are involved in plant dehydration and immune responses [26], chilling stress [27], heat stress [20], dehydration [28], infection by pathogenic bacteria [19, 21], and LP stress [29]. These findings indicate that circRNAs are abundant in plants and have important functions in response to abiotic stresses. However, circRNAs have rarely been reported in soybean [30]; furthermore, no studies have been reported on the response of soybean circRNAs to LP stress.

In this study, to quantify circRNAs in soybean and explore their potential functions in regulating the response to LP stress, we first treated two representative soybean genotypes that present contrasting P tolerances, Bogao (B) and NN94156 (sensitive and tolerant genotypes to LP stress, respectively), with two different supplies of P: LP (-P, 5 μM) and high P (HP) (+P, control, 500 μM). Using genome-wide high-throughput RNA sequencing (RNA-seq) technology, we then identified all circRNAs in the roots of soybean seedlings and validated them by real-time PCR. Our study specifically (1) identified and characterized soybean root circRNAs that are responsive to LP stress; (2) identified and characterized the possible roles of the circRNA-host genes of differentially expressed (DE) circRNAs in regulating soybean tolerance to LP stress via Gene Ontology (GO) enrichment analysis; and (3) predicted and discussed the sponge action of DE circRNAs and miRNA-targeted genes.

Materials and methods

Plant materials and hydroponic experiments

The LP-tolerant genotype Nannong 94156 (NN94156) and LP-sensitive genotype B were grown hydroponically as described previously [2]. Briefly, the seeds were germinated and grown in an artificial intelligence climate chamber under 10 h light/14 h dark and 28/20°C. When the two cotyledons had fully expanded, the soybean seedlings were transplanted into modified half-strength Hoagland's nutrient solution (pH = 5.8, 500 μM Pi, KH2PO4, sufficient P, HP). After three days, 1/2 of the seedlings were transferred to Hoagland’s nutrient solution that contained a low supply of P (5 μM P, P deficiency, LP), with the other 1/2 half the seeds remaining in the P-deficient conditions as controls. The solution was replenished every 3 d, and the soybean plants were grown in a completely randomized block design. Root tissues were collected from the treatment and control plants at 10 d after the plants were transferred to the P-deficient conditions, after which the tissues were stored at −80°C. Three independent biological replicates per genotype, each from which 12 samples were collected (HP-NR-1, 2, and 3; LP-NR-1, 2, and 3; HP-BR-1, 2, and 3; and LP-BR-1, 2, and 3), were used.

Library construction and Illumina sequencing

Total RNA from the 12 samples was extracted by Trizol (Life Technologies Inc., USA), and the quality of the RNA was verified via 1% RNase-free agarose gel electrophoresis. A Nanodrop 2000 device (Thermo-Fisher Scientific, USA) was used to quantify the RNA purity, and the RNA integrity was determined by a 2100 Bioanalyzer (Agilent Technologies, USA). First, to obtain rRNA-free residue, an Epicentre Ribo-zero rRNA Kit (Epicentre, USA, cat: MRZSR116) was used to remove the rRNA. Second, according to the instruction manual provided by Gene Denovo Biotechnology Co. (Guangzhou, China), RNA libraries that were highly strand specific were generated by a NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, USA) in conjunction with 1 μg of rRNA-depleted RNA. In brief, divalent cations were used to fragment the rRNA-depleted RNA under elevated temperature in NEBNext Reaction Buffer. First-strand cDNA was then synthesized by reverse transcriptase and random hexamer primers. Second-strand cDNA was subsequently synthesized by DNA polymerase I, second-strand synthesis reaction buffer, dUTP, and RNase H. The double-stranded cDNA was then purified and collected via AMPure XP beads for the subsequent steps. Afterward, to prepare for hybridization, the 3′ ends of the DNA fragments were repaired, and the hairpin loop structures were ligated to the fragments. To select cDNA fragments that were 200 to 500 bp in length, the fragments in each of the libraries were purified by an AMPure XP bead system. The second-strand cDNA, which contained uracil (dUTP), was digested using 3 μL of NEBNext USER Enzyme at 37°C for 15 min followed by 30 s at 98°C before PCR. Quantitative polymerase chain reaction (qPCR) was subsequently conducted in conjunction with universal PCR primers, NEBNext High-Fidelity PCR Master Mix, and Index (X) primer. To assess the library quality, an Agilent Bioanalyzer 2100 system was used to purify the PCR products (AMPure XP system). The qualified libraries were then constructed and sequenced on an Illumina HiSeq 4000 platform, and long paired-end reads (150 bp) were generated.

RNA-seq and identification of circRNAs

In this study, the FastQC program was used to check the quality of the raw sequencing reads [31], and quality control was performed using the Trimmomatic v0.30 program [32]. The reads were filtered according to three criteria: (1) reads containing adapters were removed; (2) reads containing more than 10% unknown nucleotides (nt) were removed; and (3) reads with more than 50% low-quality (Q-value ≤ 20) bases were removed. These filtered reads were subsequently aligned to the soybean reference genome, Williams 82 Wm82.a2.v1, by TopHat2 [33]. The unmapped reads were then collected and analyzed with find_circ [17] to identify circRNAs, which were annotated by searching against the circBase according to the E value via BLAST [34]. With respect to database comparisons, an E value < e-10 was defined for the existing circRNAs and, conversely, for the novel predicted circRNAs; as a result, DE circRNAs were identified between two comparison groups when the fold change (FC) was ≥ 2 (P < 0.05). Enrichment analysis of GO functions of the parent genes of DE circRNAs was performed via the GO term enrichment online tool at SoyBase (https://soybase.org/). A heatmap of the results was created using the heatmap.2 function in the gplots package [35]. When the genomic DNA sequence of annotated genes was uniquely aligned with the circRNAs containing head-to-tail splicing sites, the circRNAs were determined to originate from the annotated genes, which we refer to as the circRNA-host genes of the circRNAs. In addition, if circRNAs originated from a fragment between two genes, we considered that the circRNAs had no circRNA-host genes and are defined as “NA”.

Quantitative real-time PCR (qRT-PCR)

In this study, qRT-PCR was used to further validate the circRNA expression levels identified by RNA-seq. Thirteen randomly DE circRNAs, including nine circRNAs predicted to act as miRNA sponges, were selected for experimental validation, and the sequences of the primers used are listed in S2 Table. To maintain the reliability of the experiment, we used the same RNA samples for the transcriptome sequencing analysis to verify the RNA-seq results. Briefly, each PCR mixture consisted of real-time PCR SYBR MIX (10 μL, Toyobo, USA), first-strand cDNA (50 ng), and gene-specific primers (0.5 μL of 10 μmol L-1). The PCR amplification procedure used in this study was as follows: 95°C for 5 min, followed by 40 cycles of 95°C for 15 s and then 60°C for 60 s. The soybean internal control gene TUBULIN (GenBank accession: AY907703) was used as a control, and the cDNA template was replaced by ddH2O as a negative control. The relative gene expression data were analyzed using the 2-△△CT method [36], and three biological and technical replicates were included.

Analysis of miRNA sponges

For novel circRNAs, the software Patmatch (v1.2) (https://academic.oup.com/nar/article/33/suppl_ 2/W262/2505454) was used to predict target miRNAs within plant samples and to predict mRNAs targeted by miRNA sponges and the interaction between the mRNAs and miRNA.

Results and discussion

Identification and characteristics of circRNAs responsive to different P levels

Recently, although the involvement of circRNAs in plant tolerance to abiotic stress has been reported in several species [1921, 26], their role has rarely been reported in soybean [30], especially under LP stress. In the present study, to identify and explore how circRNAs respond to LP stress in soybean, root samples of two representative soybean genotypes, NN94156 and Bogao, which present significant differences in their traits related to P-use efficiency, were used to construct RNA-seq libraries under LP (5 μM, -P) and HP (500 μM, +P, control) conditions. In total, 12 root samples, including four samples (the roots of NN94156 (NR) and the roots of Bogao (BR) plants under two P treatments, HP and LP) with three biological replicates, were collected and sequenced. Transcriptome sequencing generated approximately 1087 million raw paired-end reads from the two soybean genotypes and two different P levels, ranging from 75.5 to 100.7 million reads per sample. After quality control was performed, more than 90% of reads with high-quality scores (Q ≥ 20) were retained and uniquely mapped to the soybean reference genome sequence (Wm82.a2.v1). The unmapped reads were then collected and processed via find_circ [17] to identify circRNAs. After this two-step analysis, a total of 371 candidate circRNAs were found, which we named novel_circ_000001 to novel_circ_000371, in these 12 samples; these candidate circRNAs were selected on the basis of at least two unique back-spliced reads and at least two biological replicate samples (S1 Table). Among these circRNAs, 300, 324, 308, and 262 were identified for HP-NR, LP-NR, HP-BR and LP-BR, respectively, across both soybean genotypes under the two different P-supply levels (Fig 1A).

thumbnail
Fig 1. Identification and characterization of circRNAs in soybean roots.

(A) Histogram of circRNAs identified in each sample under different P levels. (B) Numbers of circRNAs in each type. (C) Numbers of circRNAs from different chromosomes.

https://doi.org/10.1371/journal.pone.0227243.g001

To analyze the types of the 371 circRNAs, we compared the circRNA sequence with the genome origination. The results revealed that six types could be classified: annot_exons, antisense, exon_intron, intergenic, intronic, and one_exon, which contained 47, 144, 112, 23, 5 and 40 circRNAs, respectively (Fig 1B, S1 Table). The 144 antisense circRNAs accounted for 38.8% of the total, and 199 (53.6%) circRNAs contained an exon. These results are consistent with the findings in other species, such as 50.5% of circRNAs containing exon sequences in Arabidopsis and 85.7% in rice [29], indicating that the exon_intron type is the most common in all the types and that back-spliced exons are important for the generation of circRNAs. In addition, 204 (55%) circRNAs, including four types, intronic, annot_exons, exon_intron, and one_exon, originated from annotated protein-coding genes (Fig 1B, S1 Table). Moreover, we found that 147 of the 225 circRNA-host genes produced two or more circRNAs (S1 Table), suggesting that the same gene can produce multiple circRNAs by alternative splicing. The data of previous studies, as well as those presented here, demonstrate that a close relationship exists between the origination mechanism of circRNAs and the splicing mechanism of precursor mRNAs [37, 38].

Previous studies have shown that circRNA-host genes are unevenly distributed across different chromosomes [19]. In the present study, we also found that the circRNA-host genes were unevenly distributed across chromosomes (Fig 1C); for example, chromosome 9 contained the most (58) circRNA-host genes, followed by chromosomes 13 and 16, which contained 41 and 30 circRNA-host genes, respectively (Fig 1C). The lengths of the circRNAs ranged from approximately 90 nt to approximately 100 thousand nt, and their host gene distribution across chromosomes was also uneven (S1 Table). The mean length of the one_exon, intronic, annot_exons, intergenic, exon_intron, and antisense circRNAs was 361, 418, 574, 1409, 24607 and 7939 nt, respectively (S1 Fig., S1 Table). With respect to the length distribution of circRNAs, we found that all intronic and one_exon, most annot_exons and some antisense circRNAs were shorter than 2000 nt, whereas the majority of exon_intron and a small number of antisense and intergenic circRNAs were longer than 2000 nt (S1 Fig., S1 Table).

The number of circRNAs (371) identified in soybean roots was much lower in our study than in previous studies of model plant species and in studies that analyzed publicly available RNA-seq data; for example, several thousand circRNAs were reported to be present in Arabidopsis, rice, cucumber and soybean [29, 30, 39]. The main reason for these differences may be attributed to developmental-, tissue-, or LP stress-specific expression patterns of circRNAs [40]. For example, researchers have recently identified 776, 3171 and 2165 circRNAs in soybean leaves, roots and stems, respectively [30]. Similarly, Zhu et al. [39] identified 2787 circRNAs in cucumber, including 827 in the leaves and 2420 in the roots in response to salt stress. Thus, some of the circRNAs identified in our study might be root- or P supply-specific, although additional studies are needed. Different sequencing platforms and bioinformatics approaches might also affect the number of circRNAs detected. In cotton, 280 DE circRNAs were isolated on the basis of RNA-seq under conditions of Verticillium wilt stress [21]. Eighty-eight circRNA candidates were identified in wheat seedling leaves, and only 62 were DE under dehydration stress conditions compared with well-watered control conditions [28]. In this study, we used the RNA Ribo-minus RNA-seq method. Our results included mRNA, lncRNA, and circRNA data; however, the expression of circRNA in the examined tissues was extremely low, which may also have caused the small number of identified circRNAs.

circRNAs are genotype specific and LP stress specific in soybean roots

circRNAs whose expression is genotype specific and environment specific have been detected in plants in response to different abiotic stresses [20, 21]. We performed a Venn diagram analysis to elucidate the relationships among the circRNAs identified. Despite variation in the number of circRNAs across the four samples, a high percentage of common circRNAs (190) was shared among the four conditions, e.g., 63.3% in HP-NR, 58.6% in LP-NR, 61.7% in HP-BR, and 52.5% in LP-BR (Fig 2A, S1 Table). These circRNAs were expressed regardless of the genotype and P-supply conditions, indicating that these common responsive circRNAs are universally present in soybean roots and might be essential to root development. Notably, more circRNAs were induced in NN94156 under LP stress conditions than under the control conditions (HP), the majority of whose expression decreased in the sensitive genotype (Bogao) (Fig 1B, S1 Table). Moreover, the remaining circRNAs were determined to be genotype specific or LP-responsive specific. Briefly, 16 and 40 circRNAs were expressed specifically in NN94156 under HP conditions, and 67 and 21 were expressed specifically in Bogao under LP conditions, suggesting that, in Bogao, more circRNAs were expressed specifically under HP conditions, while in NN94156, more circRNAs were expressed specifically under LP conditions (Fig 2B). These results indicate that circRNAs may play a crucial role in soybean growth and development. In our study, the common responsive circRNAs among the four conditions might be important to ensure normal root growth, and circRNAs whose expression can be induced in both genotypes under P stress are P stress specific and could be critical for basal defense against environmental stress. Genotype-specific circRNAs have been identified in other species in response to Verticillium wilt stress [21]; thus, the roles of the genotype-specific circRNAs identified in our study merit further exploration and experimental verification.

thumbnail
Fig 2. Comparison of soybean circRNAs in each sample under different P levels.

(A) Venn diagrams of circRNAs identified in the four samples, including the roots of two representative soybean genotypes under different P levels. (B) Venn diagrams comparing expressed circRNAs in each root sample under different P levels.

https://doi.org/10.1371/journal.pone.0227243.g002

LP stress enhances the accumulation of DE circRNAs

Soybean is an important oil-bearing legume with high nutritional value [41]. Compared with the other nonleguminous species, soybean needs greater amounts of P for its growth and development, which is reflected by the greater P content in soybean seeds compared with rice, maize and wheat seeds [2]. Therefore, P deficiency is a more serious problem than other nutrient deficiencies in soybean [42]. In this study, to understand thoroughly how these circRNAs respond to LP stress, we investigated their behavior in gene expression across different genotypes and P treatments and identified DE circRNAs underlying the LP stress-related circRNAs. On the basis of the threshold (|log2FC| > 1 and P < 0.05), 120 of the 371 circRNAs were expressed at significantly different levels in four groups (including HP-NR-vs-LP-NR and HP-BR-vs-LP-BR, which denote comparisons between different P levels in the same genotype, and HP-NR-vs-HP-BR and LP_NR-vs-LP_BR, which denote comparisons between different genotypes at the same P level) (Table 1, S2 Table).

thumbnail
Table 1. Statistical analysis of intergroup significantly DE circRNAs.

https://doi.org/10.1371/journal.pone.0227243.t001

The relative expression levels of the 120 identified DE circRNAs were analyzed in each sample according to the RPM method (Fig 3A). It is clear that there were significant differences between the two genotypes. The circRNAs in the middle of the heat map may play a fundamental role in the roots of soybean plants under LP stress, while the genotype-specific circRNAs on either end of the heat map may be very important for soybean plants to respond to LP stress, which requires further study. In addition, a Venn diagram and histogram were constructed to indicate changes in expression that were specific and common among the four groups (Fig 3B and 3D). Briefly, the greatest amount of DE circRNAs occurred in the comparisons of the two soybean genotypes, NN94156 and Bogao, under both P conditions, especially under the LP conditions (Fig 3B and 3D). These results indicate that circRNAs exhibit genotype-specific regulatory patterns, which may play potential roles in tolerance to P deficiency.

thumbnail
Fig 3. DE circRNAs in the roots of soybean plant under LP stress.

(A) Heatmap and cluster analysis of circRNA expression levels in each sample. The expression levels were calculated according to log2FC values. The red/blue color indicates greater/lower levels of circRNA expression, respectively. (B) Number of DE circRNAs in each comparison. The numbers above the columns show the number of upregulated (blue) and downregulated (green) circRNAs. HP-NR-vs-LP-NR and HP-BR-vs-LP-BR denote comparisons between different P levels in the same genotype; HP-NR-vs-HP-BR and LP_NR-vs-LP_BR denote comparisons between different genotypes at the same P level. (C) Real-time PCR-based verification of the DE circRNAs on the basis of the RNA-seq data. (D) Venn diagrams of DE circRNAs identified in the different comparisons.

https://doi.org/10.1371/journal.pone.0227243.g003

In contrast, in the roots of NN94156 plants under LP stress, the expression of only three DE circRNAs was upregulated, and that of one circRNA was downregulated, while in the roots of the sensitive B plants, the expression of 15 and 13 circRNAs was upregulated and downregulated, respectively (Fig 3B and 3D). These results showed that the expression of more DE circRNAs changed significantly under LP stress in the sensitive genotype (Bogao) than in the tolerant genotype (NN94156), suggesting that the tolerant genotype was slightly affected by LP stress; in other words, compared with Bogao, NN94156 may have a better ability to maintain P homeostasis under LP stress. In addition, the circRNAs uniquely DE in each genotype or under each P level were identified (S2 Table). These results revealed that the DE circRNAs may play an important role in the response to LP stress in soybean.

We randomly selected 13 circRNAs and designed divergent primers (S3 Table) to validate the results of the DE circRNAs via qPCR experiments. The qPCR results were highly correlated (R2 = 0.82) with the RNA-seq results (Fig 3C), indicating the robustness of our study in the identification of circRNAs in response to P stress. To clarify the presence of the circRNAs further, a representative circRNA, novel_circ_000237, was selected, and a pair of divergent primers were designed. Both cDNA and gDNA (negative control) were then used as templates for reverse transcription RT-PCR amplification. The presence of the circRNAs and the reliability of the data were subsequently successfully validated by Sanger sequencing experiments (Fig 4, S2 Fig).

thumbnail
Fig 4. Validation of a circRNAs (novel_circ_000237) by qPCR and Sanger sequencing.

M, DL500 marker and the red arrow represent 100 bp. The yellow arrows denote the divergent primers for PCR amplification orientation.

https://doi.org/10.1371/journal.pone.0227243.g004

Putative functions and expression patterns of DE circRNAs and their circRNA-host genes

Studies have shown that circRNAs play an important role in gene regulation through the cis-regulatory activity of its circRNA-host genes [43]. In this study, to analyze the potential possible roles of DE circRNAs between the sensitive genotype Bogao and the tolerant genotype NN94156 under LP stress, we first predicted and annotated the circRNA-host genes of those DE circRNAs (S2 Table). The annotated genes that encode circRNAs were considered the circRNA-host genes of circRNAs, and the intergenic-type circRNA fragments originating between two genes were designated as no circRNA-host genes (NAs). As a result, among the 120 DE circRNAs, 101 were identified originating from 84 circRNA-host genes, and 19 circRNAs had no circRNA-host genes (S1 Table). Previous studies [2] have shown that the plant response to LP stress involves a complex regulatory network with many different biological processes involved, such as energy production, nucleic acid (DNA, RNA) synthesis, photosynthesis, glycolysis, respiration, cell membrane formation, redox reactions, signal transduction and nitrogen fixation [44]. In this study, GO enrichment analysis of the circRNA-host genes revealed that the enriched GO terms were related mainly to the defense response, ADP binding, nucleotide binding, organic substance catabolic processes, oxidoreductase activity, and signal transduction, suggesting the involvement of these circRNAs in the responsiveness of soybean to P deficiency (Fig 5A, S4 Table). Moreover, six circRNA-host genes were related to phosphate metabolic processes, although the associated GO terms were not significantly enriched.

thumbnail
Fig 5. GO enrichment analysis of DE circRNAs and comparisons between circRNA expression patterns and their circRNA-host genes.

(A) GO enrichment analysis of circRNA-host genes of DE circRNAs in response to LP stress. The red and green boxes highlight the gene clusters involved in response to abiotic stress and LP stress, respectively. (B) Similar expression patterns of circRNAs and their circRNA-host genes were defined as “positive”. In detail, if the expression of a circRNA and that of its circRNA-host gene are both upregulated and downregulated after LP stress, the pattern was defined as “positive”; if not, “negative”.

https://doi.org/10.1371/journal.pone.0227243.g005

Previous studies have shown that weak correlations exist between the expression of circRNAs and that of linear RNAs in both animal and human cells [38, 45]. In plants, however, studies have shown that the expression of circRNAs is positively correlated with that of their circRNA-host genes [29]. In the present study, we observed that some of the LP-responsive circRNAs were encoded by circRNA-host genes, which may play important roles under LP stress and exhibited differential expression patterns under LP stress (S2 Table). These observations revealed that the expression of circRNAs induced by LP stressed might be similar to the expression of their circRNA-host genes. We thus compared the differences in expression patterns between circRNAs and their circRNA-host genes by selecting circRNAs whose expression was significantly differentially up- or downregulated in each genotype after LP stress. If the expression patterns of a circRNA and its host gene indicated consistent upregulation or downregulation after LP stress, the circRNA was defined as a “positive” responsive circRNA; otherwise, it was defined as “negative”. Interestingly, with respect to the 32 DE circRNAs (four for HP-NR-vs-LP-NR and 28 for HP-BR-vs-LP-BR) (Table 1), a similar pattern was observed between the majority of circRNAs (28, 88%) and their circRNA-host genes under LP stress (Fig 5B). This finding suggested that the expression of most circRNAs was positively correlated with their host gene expression after LP stress in soybean.

Putative functions of circRNA-mediated miRNA networks

By acting as miRNA sponges and competing for endogenous mRNAs, circRNAs play crucial roles in regulating functional gene expression and transcription [46]. Many miRNAs have been reported to participate in the response to P deficiency in various plant species by regulating the expression of target genes, including miR399 [4750], miR319 [51], miR156 [52], and miR159 [53]. For example, a recent study in maize revealed that the PILNCR1-miR399 regulatory module is important for LP tolerance [47]. Similarly, Tian et al. [8] reported that miR399 functions as a potential integrator of phosphate homeostasis. The regulatory protein pho2 results from a nonsense mutation in the target gene of miR399, which can cause excessive phosphate accumulation [50]. Overexpression of GmMIR319 in tobacco improved tolerance to P deficiency [51]. Thus, in the present study, to reveal the functional importance of circRNAs further, we first identified all circRNA-originating target miRNAs (including circRNA-targeted miRNAs and miRNAs targeted by mRNAs). The candidate miRNA pairs were then screened on the basis of the DE circRNAs (S2 Table). We ultimately found that 70 out of 120 DE circRNAs contain miRNA-binding sites, which may function as miRNA sponges in response to P deficiency in soybean. Target prediction indicated that 570 miRNAs were targeted by these 70 circRNAs (S5 Table). Interestingly, we found that miR399, miR319, miR156 and miR159 were sponged by several DE circRNAs. Moreover, we also detected various novel DE circRNAs expressed specifically in the LP-sensitive genotype (Bogao) between the different P levels (HP-BR and LP-BR). For example, the novel circRNAs novel_circ_000013, novel_circ_000349, novel_circ_000351, and novel_circ_000277 acted as sponges of many miRNAs, suggesting that these novel DE circRNAs may play an important role in the response to LP stress (Fig 6, S6 Table), similar to how miR838 reportedly responds to abiotic stress [54]. These circRNAs may regulate their target genes in soybean in response to P deficiency by sponging these miRNA family members.

thumbnail
Fig 6. Prediction of circRNA-associated-miRNA networks in response to LP stress.

The circRNA-miRNA and miRNA-gene interaction networks were constructed on the basis of the HP-BR vs LP-BR comparison. The circles denote circRNAs, the large circular areas denote miRNAs associated with a circRNA, the rhombuses denote miRNAs, and the purple rhombuses denote different miRNA-gene interactions.

https://doi.org/10.1371/journal.pone.0227243.g006

miRNAs are involved in various plant physiological processes in response to environmental stresses by regulating the expression of functional genes [55]. To reveal the potential roles of soybean circRNAs in response to LP stress further, the target mRNAs of miRNAs were predicted by bioinformatic analysis. As a result, 95 mRNAs were identified out of the 570 putative circRNA-binding miRNAs (S7 Table). To determine the gene functions associated with the DE circRNAs, 95 targeted mRNAs were subjected to GO analysis. Interestingly, similar to that revealed by the GO analysis of the circRNA-host genes, ADP binding (GO: 0043531), defense response (GO: 0006952), nucleotide binding (GO: 0001883), signal transduction (GO: 0007165), phosphorylation (GO: 0016310), and phosphate metabolic processes (GO: 0006796) were specifically enriched (S8 Table). The GO analysis results of the predicted mRNAs revealed that the targets of the DE circRNAs in response to LP stress are involved in various functions associated with different biological processes, cellular components, and molecular functions. We speculate that the expression of circRNAs may change in response to LP stress because the circRNAs act as negative regulators of their circRNA-host genes.

In addition, to analyze the function of these DE circRNAs further, we used the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation method to analyze the predicted target mRNAs of circRNAs, and 21 pathways were revealed. These KEGG pathways are involved mainly in the seleno-compound and thiamine metabolic pathways, plant-pathogen interactions, secondary metabolism, phenylpropanoid biosynthesis, and so forth (S9 Table). Previous studies have shown that plants absorb P, sulfur, and selenium in the form of anions from the soil. A deficiency in P may affect the absorption and accumulation of sulfur and selenium by plants [56, 57]. The nucleotide-binding site-leucine-rich repeat (NBS-LRR) gene family plays a crucial role in the disease resistance and various resistance responses of plants [58]. In the present study, 13 of the 95 target genes were NBS-LRR family genes (S7 Table), suggesting that NBS-LRR family genes in soybean may play crucial roles in the response to LP stress and may be regulated by circRNAs involved in P metabolism. In addition, we found that many of the predicted targets of miRNAs encode transcription factors (TFs), including bHLH, WRKY, and nuclear TFs (S7 Table), which play crucial roles in plants under environmental stress [59]. For example, the predicted gene Glyma.13G117700 is the target gene for novel_circ_000232 and encodes WRKY TF 47, whose overexpression was proven to improve P-use efficiency in our subsequent experiments (unpublished). In addition, there are other predicted genes that may be related to P-use efficiency, such as the stress-induced protein SAM22 (Glyma.17G030300), the type IV inositol polyphosphate 5-phosphatase (Glyma.07G107400), and the TF bHLH51 (Glyma.02G258900). Overall, our results not only broaden the understanding of circRNAs in plants but also help elucidate the mechanism of circRNA responses to LP stress in soybean.

Conclusions

Our results revealed that LP stress can significantly alter the genome-wide profiles of circRNAs. The treatment-specific and genotype-specific DE circRNAs identified in this study may be involved in plant responses to LP stress via posttranscriptional regulation of miRNA networks. Moreover, our results provide some clues for further investigation of the functions of circRNAs in plant responses to LP stress.

Supporting information

S1 Fig. Distribution of different length ranges of circRNAs (e.g., 500 represents 0–500).

https://doi.org/10.1371/journal.pone.0227243.s001

(TIF)

S2 Fig. circRNA validation by RT-PCR and Sanger sequencing.

Thirteen circRNAs are shown in the figure. The successfully validated circRNAs included the following: 1, novel_circ_000274; 3, novel_circ_000338; 5, novel_circ_000035; 6, novel_circ_000093; 8, novel_circ_000108; and 11, novel_circ_000237. The unsuccessfully validated circRNAs included the following: 2, 4, 7, 9, 10, 12 and 13 (M, DL500 marker). Detailed information can be found in S3 Table.

https://doi.org/10.1371/journal.pone.0227243.s002

(TIF)

S1 Table. Detailed information of identified circRNAs in soybean roots.

https://doi.org/10.1371/journal.pone.0227243.s003

(XLSX)

S2 Table. Significant DE circRNAs detected among the four comparison groups.

https://doi.org/10.1371/journal.pone.0227243.s004

(XLSX)

S3 Table. Divergent primers used for circRNA validation.

https://doi.org/10.1371/journal.pone.0227243.s005

(XLSX)

S4 Table. GO enrichment analysis of the circRNA-host genes of significantly DE circRNAs.

https://doi.org/10.1371/journal.pone.0227243.s006

(XLSX)

S5 Table. miRNAs sponged by 70 DE circRNAs.

https://doi.org/10.1371/journal.pone.0227243.s007

(XLSX)

S6 Table. miRNAs sponged by DE circRNAs in the HP-BR and LP-BR comparisons.

https://doi.org/10.1371/journal.pone.0227243.s008

(XLSX)

S7 Table. miRNA-targeted mRNAs and their annotations.

https://doi.org/10.1371/journal.pone.0227243.s009

(XLSX)

S8 Table. GO enrichment analysis of the targeted mRNAs of significantly DE circRNAs.

https://doi.org/10.1371/journal.pone.0227243.s010

(XLSX)

S9 Table. KEGG pathway annotation of the predicted target mRNAs of circRNAs.

https://doi.org/10.1371/journal.pone.0227243.s011

(XLSX)

Acknowledgments

We thank Dr. Zhenbin Hu from Kansas State University for his critical reading of the manuscript. Three anonymous reviewers are thanked for their critical and highly valued comments.

References

  1. 1. Scheerer U, Trube N, Netzer F, Rennenberg H, Herschbach C. ATP as phosphorus and nitrogen source for nutrient uptake by fagus sylvatica and populus x canescens roots. Frontiers in Plant Science. 2019;10(378). pmid:31019519
  2. 2. Zhang D, Zhang H, Chu S, Li H, Chi Y, Triebwasser-Freese D, et al. Integrating QTL mapping and transcriptomics identifies candidate genes underlying QTLs associated with soybean tolerance to low-phosphorus stress. Plant Molecular Biology. 2017;93(1):137–50. pmid:27815671
  3. 3. Matsui A, Nguyen AH, Nakaminami K, Seki M. Arabidopsis non-coding RNA regulation in abiotic stress responses. International Journal of Molecular Sciences. 2013;14(11):22642–54. pmid:24252906
  4. 4. Chiou TJ, Aung K, Lin SI, Wu CC, Chiang SF, Su CL. Regulation of phosphate homeostasis by microRNA in Arabidopsis. Plant Cell. 2006;18(2):412–21. pmid:16387831
  5. 5. Li ZX, Zhang XR, Liu XX, Zhao YJ, Wang BM, Zhang JR. miRNA alterations are important mechanism in maize adaptations to low-phosphate environments. Plant Sci. 2016;252:103–17. pmid:27717445
  6. 6. Bao H, Chen H, Chen M, Xu H, Huo X, Xu Q, et al. Transcriptome-wide identification and characterization of microRNAs responsive to phosphate starvation in Populus tomentosa. Functional & integrative genomics. 2019. pmid:31177404
  7. 7. Pant BD, Buhtz A, Kehr J, Scheible W-R. MicroRNA399 is a long-distance signal for the regulation of plant phosphate homeostasis. Plant Journal. 2008;53(5):731–8. pmid:17988220
  8. 8. Tian L, Liu H, Ren L, Ku L, Wu L, Li M, et al. MicroRNA 399 as a potential integrator of photo-response, phosphate homeostasis, and sucrose signaling under long day condition. Bmc Plant Biology. 2018;18. pmid:30463514
  9. 9. Wang Y, Zhang F, Cui WX, Chen KQ, Zhao R, Zhang ZH. The FvPHR1 transcription factor control phosphate homeostasis by transcriptionally regulating miR399a in woodland strawberry. Plant Science. 2019;280:258–68. pmid:30824004
  10. 10. Hsieh L-C, Lin S-I, Kuo H-F, Chiou T-J. Abundance of tRNA-derived small RNAs in phosphate-starved Arabidopsis roots. Plant signaling & behavior. 2010;5(5):537–9. pmid:20404547
  11. 11. Hsieh LC, Lin SI, Shih ACC, Chen JW, Lin WY, Tseng CY, et al. Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing. Plant Physiology. 2009;151(4):2120–32. pmid:19854858
  12. 12. Shah S, Wittmann S, Kilchert C, Vasiljeva L. IncRNA recruits RNAi and the exosome to dynamically regulate pho1 expression in response to phosphate levels in fission yeast. Genes & Development. 2014;28(3):231–44. pmid:24493644
  13. 13. Garg A, Sanchez AM, Shuman S, Schwer B. A long noncoding (lnc)RNA governs expression of the phosphate transporter Pho84 in fission yeast and has cascading effects on the flanking prt lncRNA and pho1 genes. Journal of Biological Chemistry. 2018;293(12):4456–67. pmid:29414789
  14. 14. Zhang Z, Zheng Y, Ham B-K, Zhang S, Fei Z, Lucas WJ. Plant lncRNAs are enriched in and move systemically through the phloem in response to phosphate deficiency. Journal of Integrative Plant Biology. 2019;61(4):492–508. pmid:30171742
  15. 15. Yuan J, Zhang Y, Dong J, Sun Y, Lim BL, Liu D, et al. Systematic characterization of novel lncRNAs responding to phosphate starvation in Arabidopsis thaliana. Bmc Genomics. 2016;17. pmid:27538394
  16. 16. Wang TZ, Zhao MG, Zhang XX, Liu M, Yang CG, Chen YH, et al. Novel phosphate deficiency-responsive long non-coding RNAs in the legume model plant Medicago truncatula. J Exp Bot. 2017;68(21–22):5937–48. pmid:29165588
  17. 17. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8. pmid:23446348
  18. 18. Zhang Y, Zhang XO, Chen T, Xiang JF, Yin QF, Xing YH, et al. Circular Intronic long noncoding RNAs. Mol Cell. 2013;51(6):792–806. pmid:24035497
  19. 19. Zhou R, Zhu Y, Zhao J, Fang Z, Wang S, Yin J, et al. Transcriptome-wide identification and characterization of potato circular RNAs in response to pectobacterium carotovorum subspecies brasiliense infection. International Journal of Molecular Sciences. 2018;19(1). pmid:29280973
  20. 20. Pan T, Sun X, Liu Y, Li H, Deng G, Lin H, et al. Heat stress alters genome-wide profiles of circular RNAs in Arabidopsis. Plant Molecular Biology. 2018;96(3):217–29. pmid:29177640
  21. 21. Xiang L, Cai C, Cheng J, Wang L, Wu C, Shi Y, et al. Identification of circularRNAs and their targets in Gossypium under Verticillium wilt stress based on RNA-seq. Peer J. 2018;6. pmid:29576969
  22. 22. Litholdo CG Jr., da Fonseca GC. Circular RNAs and Plant Stress Responses. In: Xiao J, editor. Circular Rnas: Biogenesis and Functions. Advances in Experimental Medicine and Biology. 10872018. p. 345–53.
  23. 23. Ren Y, Yue H, Li L, Xu Y, Wang Z, Xin Z, et al. Identification and characterization of circRNAs involved in the regulation of low nitrogen-promoted root growth in hexaploid wheat. Biological Research. 2018;51. pmid:30390705
  24. 24. Wang K, Wang C, Guo B, Song K, Shi C, Jiang X, et al. CropCircDB: a comprehensive circular RNA resource for crops in response to abiotic stress. Database: the journal of biological databases and curation. 2019;2019. pmid:31058278
  25. 25. Zhao W, Chu S, Jiao Y. Present Scenario of Circular RNAs (circRNAs) in Plants. Frontiers in Plant Science. 2019;10. pmid:31001302
  26. 26. Wang Z, Liu Y, Li D, Li L, Zhang Q, Wang S, et al. Identification of circular RNAs in kiwifruit and their species-specific response to bacterial canker pathogen invasion. Frontiers in Plant Science. 2017;8(413). pmid:28396678
  27. 27. Zuo J, Wang Q, Zhu B, Luo Y, Gao L. Deciphering the roles of circRNAs on chilling injury in tomato. Biochemical and Biophysical Research Communications. 2016;479(2):132–8. pmid:27402275
  28. 28. Wang YX, Yang M, Wei SM, Qin FJ, Zhao HJ, Suo B. Identification of Circular RNAs and Their Targets in Leaves of Triticum aestivum L. under Dehydration Stress. Frontiers in Plant Science. 2017;7. pmid:28105043
  29. 29. Ye CY, Chen L, Liu C, Zhu QH, Fan LJ. Widespread noncoding circular RNAs in plants. New Phytol. 2015;208(1):88–95. pmid:26204923
  30. 30. Zhao W, Cheng YH, Zhang C, You QB, Shen XJ, Guo W, et al. Genome-wide identification and characterization of circular RNAs by high throughput sequencing in soybean. Sci Rep-Uk. 2017;7. pmid:28717203
  31. 31. Andrews S. FastQC: a quality control tool for high throughput sequence data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc. 2010.
  32. 32. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. pmid:24695404
  33. 33. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36. pmid:23618408
  34. 34. Glazar P, Papavasileiou P, Rajewsky N. circBase: a database for circular RNAs. Rna. 2014;20(11):1666–70. WOS:000344065900002. pmid:25234927
  35. 35. Warnes G, Bolker B, Bonebakker L, Gentleman R, Liaw W, Lumley T, et al. gplots: Various R programming tools for plotting data. R package version. 2009;2(4):1.
  36. 36. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(-Delta Delta C) method. Methods. 2001;25(4):402–8. pmid:11846609
  37. 37. Zhang XO, Dong R, Zhang Y, Zhang JL, Luo Z, Zhang J, et al. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome Res. 2016;26(9):1277–87. pmid:27365365
  38. 38. Gao Y, Wang JF, Zheng Y, Zhang JY, Chen S, Zhao FQ. Comprehensive identification of internal structure and alternative splicing events in circular RNAs. Nat Commun. 2016;7. pmid:27350239
  39. 39. Zhu Y-X, Jia J-H, Yang L, Xia Y-C, Zhang H-L, Jia J-B, et al. Identification of cucumber circular RNAs responsive to salt stress. BMC Plant Biology. 2019;19(1):164. pmid:31029105
  40. 40. Sablok G, Zhao HW, Sun XY. Plant circular RNAs (circRNAs): transcriptional regulation beyond miRNAs in plants. Mol Plant. 2016;9(2):192–4. pmid:26774621
  41. 41. Herridge DF, Peoples MB, Boddey RM. Global inputs of biological nitrogen fixation in agricultural systems. Plant Soil. 2008;311(1–2):1–18.
  42. 42. Zhang D, Song HN, Cheng H, Hao DR, Wang H, Kan GZ, et al. The acid phosphatase-encoding gene GmACP1 contributes to soybean tolerance to low-phosphorus stress. Plos Genet. 2014;10(1). pmid:24391523
  43. 43. Yin J, Liu M, Ma D, Wu J, Li S, Zhu Y, et al. Identification of circular RNAs and their targets during tomato fruit ripening. Postharvest Biology and Technology. 2018;136:90–8.
  44. 44. Chiou TJ, Lin SI. Signaling network in sensing phosphate availability in plants. Annu Rev Plant Biol. 2011;62:185–206. pmid:21370979
  45. 45. Veno MT, Hansen TB, Veno ST, Clausen BH, Grebing M, Finsen B, et al. Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development. Genome biology. 2015;16. pmid:26541409
  46. 46. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, et al. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384–8. pmid:23446346
  47. 47. Du Q, Wang K, Zou C, Xu C, Li W-X. The PILNCR1-miR399 regulatory module is important for low phosphate tolerance in maize. Plant Physiology. 2018;177(4):1743–53. pmid:29967097
  48. 48. Hackenberg M, Shi B-J, Gustafson P, Langridge P. Characterization of phosphorus-regulated miR399 and miR827 and their isomirs in barley under phosphorus-sufficient and phosphorus-deficient conditions. Bmc Plant Biology. 2013;13. pmid:24330740
  49. 49. Liu J-Q, Allan DL, Vance CP. Systemic signaling and local sensing of phosphate in common bean: cross-talk between photosynthate and microRNA399. Mol Plant. 2010;3(2):428–37. pmid:20147371
  50. 50. Aung K, Lin S-I, Wu C-C, Huang Y-T, Su C-L, Chiou T-J. pho2, a phosphate overaccumulator, is caused by a nonsense mutation in a MicroRNA399 target gene. Plant Physiology. 2006;141(3):1000–11. pmid:16679417
  51. 51. Chao LI, Sha AHJCJoOCS. Overexpression of Gm mIR319 in tobacco improving tolerance to phosphorus deficiency. Chinese Journal of Oil Crop Sciences. 2016.
  52. 52. Lei KJ, Lin YM, An GY. miR156 modulates rhizosphere acidification in response to phosphate limitation in Arabidopsis. Journal of Plant Research. 2016;129(2):275–84. pmid:26659856
  53. 53. Li ZY, Xu HY, Li Y, Wan XF, Ma Z, Cao J, et al. Analysis of physiological and miRNA responses to Pi deficiency in alfalfa (Medicago sativa L.). Plant Molecular Biology. 2018;96(4–5):473–92. pmid:29532290
  54. 54. Srivastava S, Srivastava AK, Suprasanna P, D'Souza SFJJoEB. Identification and profiling of arsenic stress-induced microRNAs in Brassica juncea. J Exp Bot. 2013;64(1):303–15. pmid:23162117
  55. 55. Khaldun ABM, Huang W, Lv H, Liao S, Zeng S, Wang Y. Comparative profiling of miRNAs and target gene identification in distant-grafting between tomato and lycium (Goji Berry). Frontiers in Plant Science. 2016;7. pmid:27803702
  56. 56. Zhao WL, Liang DL, Shi M, Bin HU, Xiao R, Wang JW. Effects of phosphate and selenate interactions on uptake of phosphorus and selenium by Pak Choi. Journal of Agro-Environment Science. 2013;32(12):2331–8.
  57. 57. Liu Q, Wang DJ, Jiang XJ, Cao ZH, Health. Effects of the interactions between selenium and phosphorus on the growth and selenium accumulation in rice (Oryza Sativa). J Environmental Geochemistry. 2004;26(2):325–30.
  58. 58. Takken FLW, Albrecht M, Tameling WIL. Resistance proteins: molecular switches of plant defence. Curr Opin Plant Biol. 2006;9(4):383–90. pmid:16713729
  59. 59. Chen WJ, Zhu TJTiPS. Networks of transcription factors with roles in environmental stress response. Trends in Plant Science. 2004;9(12):591–6. pmid:15564126