Genome-Wide Analysis of Long Non-Coding RNAs (lncRNAs) in Two Contrasting Soybean Genotypes Subjected to Phosphate Starvation

DOI: https://doi.org/10.21203/rs.3.rs-67625/v1

Abstract

Background: Phosphorus (P) is essential for plant growth and development, and low-phosphorus (LP) stress is a major factor limiting the growth and yield of soybean. Recently, long noncoding RNAs (lncRNAs) have been reported to be key regulators in response to stress conditions in plants. In soybean, however, how LP stress mediates biogenesis of lncRNAs remains unclear.

Results: In this study, to explore the response mechanisms of lncRNAs to LP stress, we used the roots of two representative soybean genotypes with opposite P deficiency responsiveness, a P-sensitive genotype (Bogao) and a P-tolerant genotype (NN94156), to construct RNA sequencing (RNA-seq) libraries. In total, 4,166 novel lncRNAs including 525 differently expressed (DE) lncRNAs were identified across two genotypes at different P levels. GO and KEGG analyses indicated that numerous DE lncRNAs might be involved in diverse biological processes of phosphate, such as lipid metabolic process, catalytic activity, cell membrane formation, signal transduction, nitrogen fixation. Moreover, lncRNA-mRNA-miRNA and lncRNA-mRNA networks were constructed and several promising lncRNAs were identified, which may have highly valuable for further analysis the mechanism in response to LP stress in soybean.

Conclusions: These results revealed that LP stress can significantly alter the genome-wide profiles of lncRNAs, especially for P sensitive genotype Bogao. Our findings increase the understanding and provides new insights into the function of lncNRAs responses to P stress in soybean.

Background

Long non-coding RNA (lncRNA), in general, refers to transcripts longer than 200 nucleotides and does not encode the open reading frame (ORF) [1]. In eukaryotes, most lncRNAs are transcribed by RNA polymerase II and have a structure similar to that of mRNA, such as 5′ capping, splicing and polyadenylation [2]. A growing body of evidence has shown that lncRNAs play important functional roles in diverse biological processes, such as epigenetic regulation, cell cycle regulation, cellular growth and differentiation, by regulating the level of target genes [3, 4]. LncRNAs are involved in a wide range of regulatory mechanisms impacting on gene expression, including chromatin remodeling, modulation of alternative splicing, fine-tuning of miRNA activity, and the control of mRNA translation or accumulation [5].

Recent advances in biological technologies, such as tiling arrays and RNA deep sequencing (RNA-seq), have made it possible to survey the transcriptomes of many organisms to an unprecedented degree [6]. Currently, lncRNAs have been widely identified in various plants, such as Arabidopsis thaliana [7, 8], rice [9], Zea mays [10], cotton [11]. Emerging researches revealed lncRNAs play important roles in various biological processes, including flowering regulation [12], photomorphogenesis [13], stress responses [14, 15] and other important developmental pathways [16, 17]. For example, the rice-specific lncRNA LDMAR was found to be a key gene in controlling photoperiod-sensitive male sterility [18].

Plants possess an elaborate physiological system that responds to external abiotic stress conditions [19], including phosphorus (P) deficiency. As one of the major mineral macronutrients present in all living things, P is essential for plant growth and development because of its key role in the regulation of energy metabolism and the synthesis of nucleic acids and membranes [20, 21]. Athough P is abundant in soil, it is often limited for direct use for plants due to its low bioavailability. Thus, low phosphorus (LP) stress represents a major limiting factor for plant growth and productivity [22]. P is important for plant growth and agricultural industry, however, it is estimated that P rock reserves will be depleted by 2050 [23]. Therefore, it is critical to understand the molecular mechanism of crops’ responses to LP stress and improving their phosphorus use efficiency. Plants have evolved numerous adaptive developmental and metabolic responses to cope with growth in conditions of limited phosphate, which includes modifying root system architecture (RSA), increasing the acid phosphatase activity (APA), releasing the small molecular organic acids [20]. Many studies showed that many P related genes involved in plant growth and development, such as GmACP1 [22], GmHAD1 [24], PHR1 [25]. Non-coding RNA is one of the key regulators involved in the P starvation response network. Changes in miRNAs constitute an important mechanism in plants to adapt to LP environments, such as miR399 [26] and miR827 [27]. LncRNAs also play key roles in regulating mRNA and/or miRNA levels of a large number of genes associated with P starvation responses [14, 28, 29], suggesting its important functions in regulating plant responses to LP stress. Du et al. found that PILNCR1 (long-noncoding RNA1) can inhibit ZmmiR399-guided cleavage of ZmPHO2, and the interaction between PILNCR1 and miR399 is important for tolerance to LP in maize [28].

Soybean is not only a major crop plant constituting a major worldwide agricultural industry, but also an important seed crop, which is an essential source of protein, oil and micronutrients for human and livestock consumption [30]. As containing higher concentrations of P in soybean seeds than in rice, wheat and corn, soybean requries more P than those crops to maintain growth and development [31]. Previous studies have provided an understanding of the protein-coding genes and miRNAs involved in soybean phosphate starvation [14, 28, 29], but the role of lncRNAs response to LP stress has rarely been reported in soybean.

In this study, two contrast genotypes of soybean, Bogao (a LP-sensitive genotype) and Nannong 94156 (a LP-tolerant genotype), were used to investigate their regulatory mechanism of lncRNAs in P starvation. Through genome-wide high-throughput RNA sequencing (RNA-seq) technology, we identified and characterized 4,166 lncRNAs in total that are responsive to LP stress in the roots of soybean seedlings, validated 14 lncRNAs of them by qPCR, and identified 525 differentially expressed (DE) lncRNAs in regulating soybean tolerance to LP stress. We next performed GO and KEGG analyses and constructed a LP-responsive network to explore putative functions of lncRNAs. The results lay the foundation for further in-depth understanding of the molecular mechanisms of lncRNAs’ role in response to LP stress. This study has enriched the knowledge concerning lncRNAs and provides new insights into the function of lncRNAs in LP stress.

Results

Identification and characterization of lncRNAs across two soybean genotypes under different P levels

To identify lncRNAs LP responsive expressed in soybean root under, we constructed 12 cDNA libraries from soybean root samples supplied with high/normal phosphorus ( HP, 500 µM, control ) and low phosphorus (LP, 5 µM) conditions, with two genotypes Bogao (BG, a LP-sensitive genotype) and Nannong 94156 (NN94156, a LP-tolerant genotype) with contrast responsiveness to LP stress. Three biological replicates for each condition were used to minimize individual variation. The libraries were sequenced using the Illumina HiSeq 4000 platform and 125 bp paired-end reads were generated. Approximately 1,087 million of raw sequencing reads were generated from all 12 libraries, and each sample contained reads ranging from 75.5 to 100.7 million. After discarding adaptor sequences and low-quality reads (Q-value ≤ 20), more than 90% of total reads were retained [32]. We mapped those clean reads to the soybean reference genome sequence (Wm82.a2.v1). In total, 4,166 novel lncRNAs were predicted using coding-non-coding index (CNCI) [33] and coding potential calculator (CPC) [34] in all testing conditions (Table S1).

The classification of these lncRNAs showed that, of the 4.166 lncRNAs, majority of which (2,865, 68.77%) were located in the intergenic regions, and the remaining 1,301 (31.23%) resided within the genic regions, including 454 bidirectional lncRNAs, 498 antisense lncRNAs, 121 sense lncRNAs, and 228 others that were not classified into these types (Fig. 1a). The type of lncRNA may be related to their functions, for example, overexpress LAIR (a lncRNA transcribed from the antisense of neighbouring gene LRK cluster) regulated the expression of several LRK genes and increases grain yield in rice[35]. We next analyzed the chromosomal location of all the lncRNAs on the soybean genome. The distribution of lncRNAs was uneven, with chr13 and chr18 containing greater than 250 lncRNAs, and chr05, chr11, chr16 having approximately 150 (Fig. 1b). In addition, we also analyzed the number of exons and introns in each lncRNA transcript. Most of the lncRNAs contained one exon and without intron (3,597), and the number of exon and intron can be as many as seven and six, respectively (Fig. 1c). The GC content of the lncRNAs varied greatly ranging from 20.68–64.1% with an average of 35.88%. The majority of lncRNAs have the GC percent of 30–45% (Fig. 1d). A majority (94.43%) of the lncRNAs were shorter than 2,000 nucleotides (Fig. 1e).

 

Differentially expressed (DE) lncRNAs under different P levels in two soybean genotypes

To identify the lncRNAs that were responsive to LP stress, we identified differentially expressed (DE) transcripts of the lncRNAs in the pairwise comparison between two soybean genotypes under HP and LP conditions. The FPKM (Fragments Per Kilobaseof transcript per Million mapped reads) value was used to evaluate the transcript abundance of lncRNAs. Differently expressed lncRNAs (refer as DE lncRNAs hereafter) were defined as those lncRNAs with Log2FC > 1 and FDR < 0.05. In total, 525 DE lncRNAs were identified among two different genotypes in HP and LP conditions, including 116 DE lncRNAs between different P levels in the same genotype and 456 DE lncRNAs between different genotypes at the same P level, and 47 shared DE lncRNAs (Table S2). To identify the effect of LP stress on lncRNAs, we compared the DE lncRNAs of different genotypes under the same P condition, and in same genotype at different P levels, respectively (Fig. 2). Generally, it can be seen from the volcano plot that upon LP treatment, Bogao and NN94156 have more down-regulated DE lncRNAs than up-regulated DE lncRNAs, and down-regulated DE lncRNAs has a more dramatic change in differential expression than that of up-regulated (Fig. 2a, Fig. 2b). The number and expression fold change of up-regulated and down-regulated DE lncRNAs are relatively consistent in Bogao and NN94156 under same P level (Fig. 2c, Fig. 2d, Fig. S1).

 

Because two genotypes differed greatly in response to LP stress, we performed a Venn diagram analysis to elucidate the DE lncRNAs between two genotypes under LP condition. The number of common and unique DE lncRNAs between two genotypes was indicated in the Venn diagram (Fig. 3a). NN94156 and Bogao shared 21 common DE lncRNAs in HP/LP comparisons, and Bogao contained more genotype specific DE lncRNAs (72) than that of NN94156 (23) (Fig. 3b), which is consistent with the result shown in the volcano plot (Fig. 2a, Fig. 2b). We found that the 21 common DE lncRNAs in Bogao were all down regulated in LP condition, while most of them (20) were also down regulated in NN94156, except for TCONS_00029009 (Fig. 3d). To determine whether the effect of LP stress on lncRNA is related to genotypes, we compared the changes of DE lncRNAs between Bogao and NN94156 under LP or HP levels. There were 133 and 139 DE lncRNAs unique to LP and HP condition, respectively (Fig. 3c). And the 184 common DE lncRNAs showed same trend in the up or down regulation, with the 123 down regulated and 61 up regulated (Fig. 3e).

 

Validation and quantification of lncRNAs

To validate the expression of these LP responsive lncRNAs, 14 lncRNAs were randomly selected and analyzed by quantitative PCR (qPCR). As shown in Fig. 4a, the expression patterns of the LP/HP lncRNAs as investigated by RNA-seq and qPCR were relatively consistent with similar trends. Both qPCR and RNA-seq assays presented a positive correlation in the expression fold-change with an R2 of 0.7878 (Fig. 4b), indicating the robustness of our analysis and that the lncRNA expression patterns as identified in the current study are reliable. These findings confirm that these lncRNAs are responsive to LP stress in soybean roots.

 

Functions and expression patterns of DE lncRNAs and their target genes

To reveal potential functions of the differentially expressed lncRNAs under LP stress in two contrast genotypes, we predicted the candidate targets of cis-, trans- and antisense acting DE lncRNA. In total, 785 targets for 374 DE lncRNAs were identified, and 960 pairs for one lncRNA may have several targets and/or one mRNA targets may be target of several lncRNAs (Table S3). To explore the putative functions of DE lncRNAs, we analyzed Gene Ontology (GO) terms (Table S4) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for their putative target genes (Table S5).

For DE lncRNAs in one genotype at different P levels, GO analysis revealed that 403 GO terms (195 under the biological process category, 146 under the molecular function category, and 62 under the cellular component category) were significantly enriched (P < 0.05) (Fig. 5a). For DE lncRNAs in Bogao or NN94156 supplied with the same P level, 1,086 GO terms (497 under the biological process category, 362 under the molecular function category, and 227 under the cellular component category) were significantly enriched (P < 0.05) (Fig. 5b). Although the number of GO terms in two types were different, the trend of them were relatively similar. In brief, the most significant GO terms for biological process were metabolic process, single-organism process, cellular process. As molecular functions are concerned, catalytic activity and binding were the important significantly enriched GO terms. The cell, cell part, membrane and organelle were the most important significant term for cellular components. Taken together, these results show that these lncRNAs may play roles in a variety of biological processes that are responsive to LP stress.

 

We next analyzed the enrichment of the predicted target genes of DE lncRNAs in KEGG pathway (Table S5). The targets of DE lncRNAs, same genotype at different P levels, were enriched in 42 KEGG pathways. Of these, several KEGG pathways related to carbohydrate metabolism, lipid metabolism, Amino acid metabolism were enriched (Fig. 6a). For example, propanoate metabolism, glycolysis/ gluconeogenesis, pyruvate metabolism, starch and sucrose metabolism, pentose phosphate pathway were all belonging to carbohydrate metabolism. Besides, fatty acid biosynthesis, alpha-linolenic acid metabolism and fatty acid degradation were lipid metabolism. For targets of DE lncRNAs between Bogao and NN94156 in same condition (HP and LP), 74 KEGG terms were enriched, including environmental adaptation, carbohydrate metabolism, biosynthesis of other secondary metabolites, lipid metabolism, signal transduction. Of the top 20 enriched pathways (Fig. 6b), circadian rhythm-plant and plant-pathogen interaction belonged to environmental adaptation, and were statistically significant enriched (q-values < 0.05). Terms related to three secondary metabolites pathways were enriched, including flavonoid biosynthesis, isoflavonoid biosynthesis, phenylpropanoid biosynthesis. These findings suggest that DE lncRNAs may regulate genes involved in many biological processes, including molecule metabolism, energy synthesis and signal transduction, in response to LP stress.

 

Putative P related lncRNAs based on miRNAs

MiRNAs are endogenous noncoding RNAs 20 to 24 nucleotides in size that are generated from a single-stranded RNA precursor with a hairpin secondary structure. LncRNAs can be spliced by miRNA into multiple small RNAs by which function of lncRNA can be regulated by miRNA via post-transcriptional regulation. For example miR399 is the first miRNA demonstrated to be upregulated specifically by P deficiency and rapidly decreased upon P re-addition, and certain miRNA families are commonly responsive to P deficiency among species [36].

As our research mainly focused on LP stress, thus the targets including lncRNAs and mRNAs of P related miRNAs, such as miR399, miR827, miR395, miR319, miR156, miR159, miR166, miR169, miR398 and miR447 [27, 36, 37], were selected. The lncRNA, TCONS_00090111, was identified as target of five miRNAs, including gma-miR156aa, gma-miR156z, gma-miR159b-3p, gma-miR159c, gma-miR159f-3p. Similarly, TCONS_00015352 was predicted as target of miR447-y (Table 1). Then the P related miRNA targets were predicted and we found that there were 9 mRNAs were targets of lncRNAs as well (Table S6). As shown in Fig. 7, Glyma.19G121000, Glyma.02G109500 and TCONS_00068024 (a novel identified mRNA) were targets of lncRNA and several miRNAs. Glyma.06G290000 and Glyma.12G117000 were both annotated as ethylene-responsive transcription factor 9-like mRNA and were targets of two lncRNAs, TCONS_00030280 and TCONS_00068008, respectively. Both genes were also targets of gma-miR169l-3p, which is a P related miRNA [38]. Another example came from Glyma.19G193900 that was predicted as purple acid phosphatase 22-like, and it is the target of TCONS_00105416 and miR398-x, which belonged to miR398 with demonstrated role in coping with P starvation stress [36].

Table 1

The identified lncRNAs as targets of P related miRNAs

MiRNA

Target lncRNA

gma-miR156aa

TCONS_00090111

gma-miR156z

TCONS_00090111

gma-miR159b-3p

TCONS_00090111

gma-miR159c

TCONS_00090111

gma-miR159f-3p

TCONS_00090111

 

Transcription factors, P-related and planthormone associated lncRNA-mRNA network construction

TFs regulate a diverse group of genes during stress responses and are important components of gene regulatory networks, and many TFs belonging to some families have been proved to play an important role in the maintenance of P homeostasis, such as PHR (phosphate starvation response), bHLH, WRKY, ZAT, MYB and so on [39]. P regulation of root system architecture is driven by the local perception of PO4− at the root tip and involves changes in multiple plant hormones, such as auxin, gibberellins and ethylene, the hormonal changes coordinate together the root developmental responses to P availability [37]. And P related genes such as PHO2, PHR1 play important roles in P starvation response. To further study the function of lncRNAs in the process of soybean roots in response to LP stress, we constructed a lncRNA-mRNA network of interested mRNAs (including transcription factors, P-related and planthormone targets) and corresponding lncRNAs, according to GO, KEGG and functional annotation of target genes (Table S7). As shown in Fig. 8, the lncRNA-mRNA network consisted 52 lncRNAs and 109 targets in total. 23 lncRNAs might be involved in the regulation of gene transcription, because their target genes have transcription factor activity, of them 3 lncRNAs each have two transcription factor targets, and 20 lncRNAs with only one. And the number of lncRNAs targets varied from one to five, Glyma.02G226800 and Glyma.02G226700 were targets of 5 and 4 lncRNAs. And 26 transcription factors belonged to diverse families such as MYB, bHLH, NAC, AP2. Of them, MYB and bHLH have been reported to play roles in the maintenance of P homeostasis [39]. Interestingly, we found 8 lncRNAs might be involved in ethylene regulation, as their targets were annotated as ethylene responsive transcription factor. Among the P-related genes, Glyma.17G172700 and Glyma.19G193900 were annotated as purple acid phosphatase. And Glyma.20G021600 were predicted as phosphate transporter, which were known as PHTs and involved in LP stress.

 

Discussion

LncRNAs play important roles in a wide range of biological processes, especially in plant reproductive development and response to stresses [40]. However, little is known about their roles in LP stress, which is a major limiting factor for plant growth and agricultural industry. Here, we undertook a genome wide ananlysis of lncRNAs in two contrasting soybean genotypes subjected to phosphate starvation.

The number of lncRNAs varies greatly across plant species. For example, 48,345 lncRNAs were identified in maize transcriptome [15], 1,212 novel lncRNAs were found in Arabidopsis seedlings grown under P-sufficient and P-deficient conditions [14]. In this study, 4,166 lncRNAs were identified. The lncRNAs identified here share most of the common features of lncRNAs reported with those in other plants, such as short length, single exons, low GC percent, which may be responsible for the common and ancient evolutionary origin. In addition, we found that the sequence length of several (70) lncRNAs were longer than 3,000 bp in soybean roots, which is similar to previous research in which 285 lncRNAs length were longer than 3,000 bp and the length of 28 lncRNAs were more than 10,000 bp [41], indicating that a small number of long lncRNAs exist in plants. The type of lncRNAs were also highly variable in plants. We identified more intergenic lncRNAs (2,865, 68.77%) than other types including antisense lncRNAs (498, 11.95%), in contrast, the number of antisense lncRNAs were greater than other types in Arabidopsis [14]. We are curious whether the number of exons is related to the length of the gene in lncRNA. After inspection, we found that the length of lncRNAs with a large number of exons did not increase significantly.

To identify the lncRNAs that were responsive to P stress, we identified differentially expressed (DE) transcripts of the lncRNAs in the pairwise comparison between two soybean genotypes under HP and LP conditions, and 525 DE lncRNAs were identified among two different genotypes in HP and LP conditions in total. As shown in Fig. 2a, Fig. 2b and Fig. S1, the number of DE lncRNAs under LP stress in Bogao is more than that in NN94156, indicating that the Bogao is more sensitive to LP treatment, which is consistent with our previous research in which Bogao was a P-sensitive genotype, while NN94156 was a P-tolerant genotype [30, 32]. Furthermore, previous and our results reveals that all the mRNA, circular RNA and lncRNA, Bogao is more sensitive to LP stress than that of NN94156 [30, 32]. The two genotypes differed greatly in response to LP stress, NN94156 and Bogao shared 21 common DE lncRNAs in HP/LP comparisons, and Bogao contained more genotype specific DE lncRNAs (72) than that of NN94156 (23) (Fig. 3b). Most of the common DE lncRNAs were consistitutively down regulated in two genotypes, indicating conserved biological mechanism for lncRNAs that involved in basal responsiveness to LP stress in both soybean genotypes. Although the similar trend of response to LP stress for these shared lncRNAs, the expression change degree of response is significantly different (Fig. 3d, Fig. 3e). To determine whether the effect of LP stress on lncRNA is related to genotypes, we compared the changes of DE lncRNAs between Bogao and NN94156 under LP or HP levels. We identified 317 DE lncRNAs (181 down-regulated, 136 up-regulated) in LP condition, suggesting that they were constitutively but differentially expressed between the two genotypes under LP condition. The 133 DE lncRNAs unique to LP condition may be play a role in LP tolerance. In contrast, 139 lncRNAs were differentially expressed in two genotypes only in HP condition, suggesting that they are differentially expressed specific to LP stress.

Studies have shown that lncRNA can directly bind to mRNA by which affecting translation, shearing, and degradation of mRNA, and can also indirectly influence the expression of mRNA [17]. Thus far, the mechanism of interaction between lncRNA and mRNA has not yet become clear. To reveal potential functions of the differentially expressed lncRNAs under LP stress in two contrast genotypes, we predicted the candidate targets of DE lncRNA, then analyzed GO terms and KEGG pathways for their putative target genes (Table S5). For DE lncRNAs in one genotype at different P levels and DE lncRNAs in Bogao or NN94156 supplied with the same P level, 403 and 1,086 GO terms with 42 and 74 KEGG pathways were significantly enriched (P < 0.05), respectively (Fig. 5). Our enrichment results showed that LP stress is a complex regulatory network involved in diverse biological processes such as lipid metabolic process, catalytic activity, cell membrane formation, signal transduction, nitrogen fixation (Fig. 5, Fig. 6), which was supported by previous studies focusing LP in soybean [32, 42]. Previous research showed that NtMYB12 acts as a phosphorus starvation response enhancement factor and regulates NtCHS and NtPT2 expression, which results in increased flavonol and P accumulation and enhances tolerance to LP stress [43]. In this study, targets of lncRNAs were enriched as KEGG pathway including flavonoid, isoflavonoid and phenylpropanoid biosynthesis, indicating that those lncRNAs may be involved in secondary metabolites to regulate P responsive genes, which needs further research.

LncRNAs can be spliced by miRNA into multiple small RNAs by which function of lncRNA can be regulated by miRNA via post-transcriptional regulation [36]. Targets of ten P related miRNAs (miR399, miR827, miR395, miR319, miR156, miR159, miR166, miR169, miR398 and miR447) were predicted, and nine of them were targets of lncRNAs as well (Table S6). Two lncRNAs, TCONS_00030280 and TCONS_00068008, possessed shared mRNA targets (Glyma.06G290000 and Glyma.12G117000, annotated as ethylene-responsive transcription factor 9-like) with the P-related miR169l-3p (Fig. 7). This finding was further supported by our recent study where NN94156 has ability to tolerate LP stress by ethylene regulator-mediated enhanced P uptake and use efficiency in roots [44]. Therefore, lncRNAs might, in part, involved in ethylene-mediated LP stress tolerance and both genes were candidate genes that merit further investigation to gain further understanding of how lncRNA involved the LP stress tolerance. Glyma.19G193900, which was predicted as purple acid phosphatase 22-like, is the target of TCONS_00105416 and miR398-x, which belonged to miR398 with demonstrated role in coping with P starvation stress [36]. Purple acid phosphatases is widely recognized as an adaption of plants to phosphorus deficiency and the secretion of PAPs (purple acid phosphatases) play important roles in P acquisition [45]. Identification of other genes with our preliminary scenario suggested that lncRNAs involve in LP stress responsive through the manipulation of genes with a variety of functionality and many of them that might also be the co-target of P associated miRNA. Detailed investigation on those genes may gain increased understanding.

Enhanced root hair production, which increases the root surface area for nutrient uptake, is a typical adaptive response of plants to phosphate starvation [46]. And ethylene plays an important role in root hair development induced by P starvation via controlling root hair elongation [37]. According to GO, KEGG and functional annotation of targets genes of DE lncRNAs, transcription factors, P-related and planthormone targets were selected to construct lncRNA-mRNA network (Fig. 8). Interestingly, we found 8 lncRNAs might be involved in ethylene regulation, as their targets were annotated as ethylene responsive transcription factor. Ethylene-responsive transcription factor belongs to APETALA2 (AP2)/ Ethylene Response Factor (ERF) family, which existed widely in plants. In Arabidopsis, the AP2/ERF TF superfamily comprises 147 members, and AP2/ERF proteins are known to regulate the responses of plants to various biotic and abiotic stresses and developmental processes. RNA interference and overexpressing of AtERF070 (AT1G71130), an ethylene response factor, resulted in altered morphophysiological traits of roots, and express change of a number of P starvation responsive genes, suggesting a potential role of this TF in maintaining P homeostasis [39]. Besides ethylene, auxin, GA and salicylic acid may be involved in response to LP stress in soybean. Induction and secretion of acid phosphatases (APases) is considered to be an important strategy for improving plant growth under conditions of low inorganic phosphate. PAPs, are an important class of plant APases that could be secreted into the rhizosphere to utilize organic phosphorus for plant growth and development [47]. Among the P-related genes, Glyma.17G172700 and Glyma.19G193900 were annotated as purple acid phosphatase. And Glyma.20G021600 were predicted as phosphate transporter, which were known as PHTs, and have important roles in P acquisition, allocation, and signal transduction. We speculate their corresponding lncRNAs may have similar functions and be involved in LP stress.

Conclusion

The main aim of this research was to identify the potential lncRNAs response to LP stress in soybean, and the differences in lncRNAs response between different P efficiency genotypes. Our results showed that a total of 4,166 lncRNAs, including 525 DE lncRNAs were identified by using roots of two representative genotypes upon HP and LP conditions. LP stress can alter the genome-wide expression levels of lncRNAs, especially for P-sensitive genotypes Bogao. These finding may provide a first look at the landscape of lncRNAs of soybean in response to LP stress. Moreover, we have identified several promising lncRNAs, which may have potential value for further analysis the mechanism in response to LP stress in soybean. Overall, this study has enriched the knowledge concerning lncRNAs and provides some clues to explore the function of lncRNAs in LP stress in soybean.

Methods

Plant materials and growth conditions

The LP-tolerant genotype Nannong 94–156 (NN94156) and LP-sensitive genotype Bogao were grown hydroponically as described previously [42]. In Brief, the seeds were surface-sterilized and germinated and grown in an artificial intelligence climate chamber at 28/20˚C and with a photoperiod of 10/14 h of light/dark. 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, KH2PO4, sufficient P, HP). Three days later, half of the seedlings were transferred to a low supply of P (5 µM P, LP) Hoagland’s nutrient solution, the other half the seeds remained in the P-deficient conditions. The soybean plants were placed in the hydroponics box using a completely randomized block design. The solution was replenished every 3 d, and at 10 d after the plants were transferred to the P-deficient conditions, three independent biological replicates of the roots of seedlings were collected (12 samples in total) and stored at -80˚C for extraction of total RNA.

RNA extraction, library construction, and Illumina sequencing

The process of RNA extraction and purity, library construction, and Illumina sequencing were according to Lv et al. [32]. Briefly, after total RNA was extracted, rRNAs were removed using an Epicentre Ribo-zero rRNA Kit (Epicentre, USA, cat: MRZSR116). Then, 1 µg rRNA-depleted RNA per sample were used to generate sequencing libraries according to the manual provided by Gene Denovo Biotechnology Co. (Guangzhou, China). The qualified libraries were then constructed and sequenced on an Illumina HiSeq 4000 platform. The 12 gene expression libraries were named as HP-NR-1, HP-NR-2, HP-NR-3; LP-NR-1, LP-NR-2, LP-NR-3; HP-BR-1, HP-BR-2, HP-BR-3; LP-BR-1, LP-BR-2, and LP-BR-3.

Identification of lncRNAs

Raw data were preprocessed to filter out adapters, reads containing more than 10% of unknown nucleotides and reads containing more than 50% of low quality (Q-value ≤ 20) bases were removed. Then, the obtained clean reads were aligned to soybean reference genome, Williams 82 Wm82.a2.v1, using the splice read aligner TopHat2 [48]. The length of transcripts longer than 200 bp and the exon number more than one were selected. Two softwares CNCI (coding-non-coding index) [33] and CPC (coding potential calculator) [34] were used to assess the protein-coding potential of transcripts by default parameters. The intersection of both non protein-coding potential results were chosen as long non-coding RNAs (lncRNAs). LncRNAs differential expression analysis was performed by DESeq [49] software between two different groups. The transcripts with the parameter of false discovery rate (FDR) below 0.05 were considered differentially expressed genes. Venn diagrams were generated using Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

Quantitative PCR (qPCR) validation of lncRNAs

To validate the expression data from RNA-seq, 14 lncRNAs were selected randomly for quantitative PCR (qPCR) analysis. And qPCR experiments were performed on an ABI 7500 system (Applied Biosystems, Foster City, CA, USA). Each PCR reaction contained 10 µL qPCR SYBR MIX (Toyobo, USA), 50 ng cDNA and 0.5 µL of 10 µmol L− 1 gene-specific primers. The PCR amplification procedure was 95˚C for 5 min, followed by 40 cycles of 95˚C for 15 s, then 60˚C for 60 s. The Tubulin (GenBank accession: AY907703) gene in soybean was used as internal control, and the cDNA template replaced by ddH2O was used as negative control. This experiment was performed with three technical replicates and three biological replicates, and the relative expression of lncRNA were analyzed using the 2−ΔΔCT method [50]. The genes and their primers are listed in Table S8.

Target gene prediction and functional analysis

We searched for coding genes 10 kb upstream and downstream of the identified lncRNAs and then predicted them as cis-acting lncRNAs target neighboring genes. Some antisense lncRNAs may regulate gene silencing, transcription and mRNA stability. The software RNAplex [51] (https://www.tbi.univie.ac.at/RNA/RNAplex.1.html) was used to predict the complementary correlation of antisense lncRNA and mRNA, and mRNAs were predicted as antisense genes of lncRNAs. Another function of lncRNAs is trans-regulation of co-expressed genes not adjacent to lncRNAs. The correlation of expression between lncRNAs and protein-coding genes were analyzed to identify trans-acting genes of lncRNAs. These lncRNA target genes were functionally annotated using the GO (http://geneontology.org/) and KEGG (http://www.genome.jp/kegg/) databases.

Analyses of lncRNAs and/or mRNAs with miRNAs

To find potential miRNA precursor, lncRNAs were aligned to miRBase and those with more than 90% were selected. In addition, the software miRPare, which is based on SVM method was also used to predict miRNA precursors. The software patmatch (v1.2) was used to predict target genes of miRNA. The lncRNA-mRNA-miRNA network was constructed by Cytoscape 3.8.0 [52] software.

Abbreviations

lncRNAs: long non-coding RNAs; ORF: Open reading frame; APA: acid phosphatase activity; FDR: False Discovery Rate; DE: differently expressed; RNA-seq: RNA sequencing; FPKM: Fragments Per Kilobase of transcript per Million mapped reads; RSA: root system architecture; P: Phosphorus; LP: low-phosphorus; CNCI: coding-non-coding index; CPC: coding potential calculator; qPCR: quantitative PCR; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; PAPs: purple acid phosphatases; AP2/ERF: Apetala2/Ethylene Response Factor ; TFs: Transcription factors

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

All data generated or analyzed during this study are included in this article and its additional files.

Competing interests

The authors declare that they have no competing interest.

Funding

This research was funded by the Ministry of Science and Technology of China (2016YFD0100500), the key scientific and technological project of Henan Province (192102110023), the Henan agricultural university science and technology innovation fund (KJCX2019C02) and the Key Scientific Research Projects of Higher Education Institutions in Henan Province (15B210007, 20A210017, 18B210004).

Authors’ Contributions

DZ and ZH conceived and designed the experiments. XZ, YY and HX conducted the experiment. JZ and DZ performed data analysis. JZ and XZ wrote the manuscript. All authors read and approved the manuscript.

Acknowledgments

The authors thank to lab members for assistance.

Authors' information

1 Collaborative Innovation Center of Henan Grain Crops, College of Agronomy, Henan Agricultural University, Zhengzhou, China. 2 School of Life Science and Technology, Henan Institute of Science and Technology/Henan Collaborative Innovation Center of Modern Biological Breeding, Xinxiang 453003, China.

References

  1. Djebali S, Davis CA, Merkel A, Dobin A, Lassmann T, Mortazavi A, Tanzer A, Lagarde J, Lin W, Schlesinger F, et al. Landscape of transcription in human cells. Nature. 2012;489(7414):101–8.
  2. Marchese FP, Raimondi I, Huarte M. The multidimensional mechanisms of long noncoding RNA function. Genome Biology. 2017; 18.
  3. Rinn JL, Chang HY. Genome Regulation by Long Noncoding RNAs. Annu Rev Biochem. 2012;81:81:145–66.
  4. Mishra A, Bohra A. Non-coding RNAs and plant male sterility: current knowledge and future prospects. Plant Cell Rep. 2018;37(2):177–91.
  5. Ariel F, Romero-Barrios N, Jegu T, Benhamed M, Crespi M. Battles and hijacks: noncoding transcription in plants. Trends Plant Sci. 2015;20(6):362–71.
  6. Moran VA, Perera RJ, Khalil AM. Emerging functional and mechanistic paradigms of mammalian long non-coding RNAs. Nucleic Acids Res. 2012;40(14):6391–400.
  7. Liu ZW, Zhao N, Su YN, Chen SS, He XJ. Exogenously overexpressed intronic long noncoding RNAs activate host gene expression by affecting histone modification in Arabidopsis. Sci Rep. 2020;10(1):3094.
  8. Wang T, Xing J, Liu Z, Zheng M, Yao Y, Hu Z, Peng H, Xin M, Zhou D, Ni Z. Histone acetyltransferase GCN5-mediated regulation of long non-coding RNA At4 contributes to phosphate starvation response in Arabidopsis. J Exp Bot. 2019;70(21):6337–48.
  9. Chen L, Shi S, Jiang N, Khanzada H, Wassan GM, Zhu C, Peng X, Xu J, Chen Y, Yu Q, et al. Genome-wide analysis of long non-coding RNAs affecting roots development at an early stage in the rice response to cadmium stress. BMC Genomics. 2018;19(1):460.
  10. Li L, Eichten SR, Shimizu R, Petsch K, Yeh CT, Wu W, Chettoor AM, Givan SA, Cole RA, Fowler JE, et al. Genome-wide discovery and characterization of maize long non-coding RNAs. Genome Biol. 2014;15(2):R40.
  11. Zhao T, Tao X, Feng S, Wang L, Hong H, Ma W, Shang G, Guo S, He Y, Zhou B, et al. LncRNAs in polyploid cotton interspecific hybrids are derived from transposon neofunctionalization. Genome Biol. 2018;19(1):195.
  12. Zhao X, Li J, Lian B, Gu H, Li Y, Qi Y. Global identification of Arabidopsis lncRNAs reveals the regulation of MAF4 by a natural antisense RNA. Nat Commun. 2018;9(1):5056.
  13. Wang Y, Fan X, Lin F, He G, Terzaghi W, Zhu D, Deng XW. Arabidopsis noncoding RNA mediates control of photomorphogenesis by red light. Proc Natl Acad Sci U S A. 2014;111(28):10359–64.
  14. Yuan J, Zhang Y, Dong J, Sun Y, Lim BL, Liu D, Lu ZJ. Systematic characterization of novel lncRNAs responding to phosphate starvation in Arabidopsis thaliana. BMC Genomics. 2016;17:655.
  15. Huanca-Mamani W, Arias-Carrasco R, Cardenas-Ninasivincha S, Rojas-Herrera M, Sepulveda-Hermosilla G, Caris-Maldonado JC, Bastias E, Maracaja-Coutinho V. Long Non-Coding RNAs Responsive to Salt and Boron Stress in the Hyper-Arid Lluteno Maize from Atacama Desert. Genes (Basel). 2018; 9(3).
  16. Wang HV, Chekanova JA. Long Noncoding RNAs in Plants. Adv Exp Med Biol. 2017;1008:133–54.
  17. Chekanova JA. Long non-coding RNAs and their functions in plants. Curr Opin Plant Biol. 2015;27:207–16.
  18. Ding J, Lu Q, Ouyang Y, Mao H, Zhang P, Yao J, Xu C, Li X, Xiao J, Zhang Q. A long noncoding RNA regulates photoperiod-sensitive male sterility, an essential component of hybrid rice. Proc Natl Acad Sci U S A. 2012;109(7):2654–9.
  19. Hirayama T, Shinozaki K. Research on plant abiotic stress responses in the post-genome era: past, present and future. Plant J. 2010;61(6):1041–52.
  20. Puga MI, Rojas-Triana M, de Lorenzo L, Leyva A, Rubio V, Paz-Ares J. Novel signals in the regulation of Pi starvation responses in plants: facts and promises. Curr Opin Plant Biol. 2017;39:40–9.
  21. Wang F, Deng M, Xu J, Zhu X, Mao C. Molecular mechanisms of phosphate transport and signaling in higher plants. Semin Cell Dev Biol. 2018;74:114–22.
  22. Zhang D, Song H, Cheng H, Hao D, Wang H, Kan G, Jin H, Yu D. The acid phosphatase-encoding gene GmACP1 contributes to soybean tolerance to low-phosphorus stress. PLoS Genet. 2014;10(1):e1004061.
  23. Cordell D, Drangert J-O, White S. The story of phosphorus: Global food security and food for thought. Glob Environ Change. 2009;19(2):292–305.
  24. Cai Z, Cheng Y, Xian P, Ma Q, Wen K, Xia Q, Zhang G, Nian H. Acid phosphatase gene GmHAD1 linked to low phosphorus tolerance in soybean, through fine mapping. Theor Appl Genet. 2018;131(8):1715–28.
  25. Rubio V, Linhares F, Solano R, Martin AC, Iglesias J, Leyva A, Paz-Ares J. A conserved MYB transcription factor involved in phosphate starvation signaling both in vascular plants and in unicellular algae. Genes Dev. 2001;15(16):2122–33.
  26. Fujii H, Chiou TJ, Lin SI, Aung K, Zhu JK. A miRNA involved in phosphate-starvation response in Arabidopsis. Curr Biol. 2005;15(22):2038–43.
  27. Li Z, Zhang X, Liu X, Zhao Y, Wang B, Zhang J. miRNA alterations are important mechanism in maize adaptations to low-phosphate environments. Plant Sci. 2016;252:103–17.
  28. Du Q, Wang K, Zou C, Xu C, Li WX. The PILNCR1-miR399 Regulatory Module Is Important for Low Phosphate Tolerance in Maize. Plant Physiol. 2018;177(4):1743–53.
  29. Hu B, Zhu C, Li F, Tang J, Wang Y, Lin A, Liu L, Che R, Chu C. LEAF TIP NECROSIS1 plays a pivotal role in the regulation of multiple phosphate starvation responses in rice. Plant Physiol. 2011;156(3):1101–15.
  30. Li H, Yang Y, Zhang H, Chu S, Zhang X, Yin D, Yu D, Zhang D. A Genetic Relationship between Phosphorus Efficiency and Photosynthetic Traits in Soybean As Revealed by QTL Analysis Using a High-Density Genetic Map. Front Plant Sci. 2016;7:924.
  31. XiHuan L, WenSuo C, CaiYing Z. Advances of soybean (Glycine max L.) phosphorus nutrition and high P-eicient germplasms screening in China Soybean Science. 2011; 30(02):322–327.
  32. Lv L, Yu K, Lü H, Zhang X, Liu X, Sun C, Xu H, Zhang J, He X, Zhang D. Transcriptome-wide identification of novel circular RNAs in soybean in response to low-phosphorus stress. PLoS One. 2020;15(1):e0227243.
  33. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.
  34. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007; 35(Web Server issue):W345-349.
  35. Wang Y, Luo X, Sun F, Hu J, Zha X, Su W, Yang J. Overexpressing lncRNA LAIR increases grain yield and regulates neighbouring gene cluster expression in rice. Nature Communications. 2018;9(1):3516.
  36. Kuo HF, Chiou TJ. The role of microRNAs in phosphorus deficiency signaling. Plant Physiol. 2011;156(3):1016–24.
  37. Oldroyd GED, Leyser O. A plant’s diet, surviving in a variable nutrient environment. Science. 2020; 368(6486).
  38. Ning L-H, Du W-k, Song H-N, Shao H-B, Qi W-C, Sheteiwy MSA, Yu D-y. Identification of responsive miRNAs involved in combination stresses of phosphate starvation and salt stress in soybean root. Environ Exp Bot. 2019;167:103823.
  39. Ramaiah M, Jain A, Raghothama KG. ETHYLENE RESPONSE FACTOR070 Regulates Root Development and Phosphate Starvation-Mediated Responses. Plant Physiol. 2014;164(3):1484.
  40. Zhang YC, Chen YQ. Long noncoding RNAs: new regulators in plant development. Biochem Biophys Res Commun. 2013;436(2):111–4.
  41. Xu W, Yang T, Wang B, Han B, Zhou H, Wang Y, Li DZ, Liu A. Differential expression networks and inheritance patterns of long non-coding RNAs in castor bean seeds. Plant J. 2018;95(2):324–40.
  42. Zhang D, Zhang H, Chu S, Li H, Chi Y, Triebwasser-Freese D, Lv H, Yu D. Integrating QTL mapping and transcriptomics identifies candidate genes underlying QTLs associated with soybean tolerance to low-phosphorus stress. Plant Mol Biol. 2016.
  43. Song Z, Luo Y, Wang W, Fan N, Wang D, Yang C, Jia H. NtMYB12 Positively Regulates Flavonol Biosynthesis and Enhances Tolerance to Low Pi Stress in Nicotiana tabacum. Front Plant Sci. 2020;10:1683–3.
  44. Zhang H, Yang Y, Sun C, Liu X, Lv L, Hu Z, Yu D, Zhang D. Up-regulating GmETO1 improves phosphorus uptake and use efficiency by promoting root growth in soybean. Plant Cell Environ. 2020.
  45. Li C, Li C, Zhang H, Liao H, Wang X. The purple acid phosphatase GmPAP21 enhances internal phosphorus utilization and possibly plays a role in symbiosis with rhizobia in soybean. Physiol Plant. 2017;159(2):215–27.
  46. Song L, Yu H, Dong J, Che X, Jiao Y, Liu D. The Molecular Mechanism of Ethylene-Mediated Root Hair Development Induced by Phosphate Starvation. PLOS Genetics. 2016;12(7):e1006194.
  47. Kong Y, Li X, Wang B, Li W, Du H, Zhang C. The Soybean Purple Acid Phosphatase GmPAP14 Predominantly Enhances External Phytate Utilization in Plants. Front Plant Sci. 2018; 9(292).
  48. 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.
  49. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
  50. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–8.
  51. Tafer H, Hofacker IL. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics. 2008;24(22):2657–63.
  52. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.