Next Article in Journal
Nephrotoxic Potential of Putative 3,5-Dichloroaniline (3,5-DCA) Metabolites and Biotransformation of 3,5-DCA in Isolated Kidney Cells from Fischer 344 Rats
Next Article in Special Issue
Plant Morphological, Physiological and Anatomical Adaption to Flooding Stress and the Underlying Molecular Mechanisms
Previous Article in Journal
Trifostigmanoside I, an Active Compound from Sweet Potato, Restores the Activity of MUC2 and Protects the Tight Junctions through PKCα/β to Maintain Intestinal Barrier Function
Previous Article in Special Issue
Compensation Mechanism of the Photosynthetic Apparatus in Arabidopsis thaliana ch1 Mutants
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 22
Issue 1
10.3390/ijms22010289
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Whole-Transcriptome RNA Sequencing Reveals the Global Molecular Responses and CeRNA Regulatory Network of mRNAs, lncRNAs, miRNAs and circRNAs in Response to Salt Stress in Sugar Beet (Beta vulgaris)

by
Junliang Li
,
Jie Cui
*,
Cuihong Dai
,
Tianjiao Liu
,
Dayou Cheng
* and
Chengfei Luo
School of Chemistry and Chemical Engineering, Harbin Institute of Technology, Harbin 150086, China
*
Authors to whom correspondence should be addressed.
Int. J. Mol. Sci. 2021, 22(1), 289; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22010289
Submission received: 2 December 2020 / Revised: 25 December 2020 / Accepted: 27 December 2020 / Published: 30 December 2020
(This article belongs to the Special Issue Environmental Stress and Plants)
Browse Figures
Review Reports Versions Notes

Abstract

:
Sugar beet is an important sugar-yielding crop with some tolerance to salt, but the mechanistic basis of this tolerance is not known. In the present study, we have used whole-transcriptome RNA-seq and degradome sequencing in response to salt stress to uncover differentially expressed (DE) mRNAs, microRNAs (miRNAs), long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) in both leaves and roots. A competitive endogenous RNA (ceRNA) network was constructed with the predicted DE pairs, which revealed regulatory roles under salt stress. A functional analysis suggests that ceRNAs are implicated in copper redistribution, plasma membrane permeability, glycometabolism and energy metabolism, NAC transcription factor and the phosphoinositol signaling system. Overall, we conducted for the first time a full transcriptomic analysis of sugar beet under salt stress that involves a potential ceRNA network, thus providing a basis to study the potential functions of lncRNAs/circRNAs.

1. Introduction

Crop growth and productivity is limited by high soil salinity, which affects plant physiological and biochemical processes [1,2]. Indeed, soil salinity leads to toxic levels of Na+ and Cl− ions, increasing osmotic stress [3,4]. An excess of salt in the soil also reduces the water potential on the root surface and prevents water absorption [5]. Further, a high cytoplasmic Na+ concentration can disrupt the uptake of important ions such as potassium (K+), with adverse effects on enzyme catalytic activity and various metabolic pathways [6,7,8]. Osmotic and ionic stress can cause secondary stress in plants by accumulation of toxic compounds and reactive oxygen species (ROS) [4]. In response to these, plants have evolved mechanisms that include aquaporins [9], transporters [10] and various enzymes involved in free radical scavenging.
In addition to protein-coding RNAs, gene regulation by non-coding RNAs is also essential for plant adaptation to stressful conditions [11,12]. MicroRNAs (miRNAs) are small, endogenous, non-coding RNAs ranging from 19 to 24 nucleotides (nts) that regulate post-transcriptional events by cleaving their target mRNAs or by preventing translation [13]. Increasing evidence implicates miRNA in the gene regulation during seedling growth [14], floral development [15], male and female sterility [16,17] and responses to environmental stress [18]. Also, exported miRNAs have been shown to inhibit virulent gene expression by a plant fungal pathogen [19].
Long non-coding RNAs (lncRNAs) are another type of noncoding RNA that is longer than 200 nts [20]. In general, lncRNAs act via epigenetic modification, transcriptional regulation and post transcriptional regulation [21]. Recent genome-wide studies have revealed that they may be involved in plant stress; in Populus trichocarpa, a total of 504 lncRNAs were found to be drought responsive [22], and, in Gossypium hirsutum, 44 lncRNAs were responsive to salt stress [23]. Under salt stress, 185 differentially expressed lncRNAs have been reported in duckweed [24], whereas, in cotton, lncRNA973 regulates several genes, some of which are involved in oxygen scavenging [25]. In apples (Malus domestica), lncRNA MSTRG.85814.11 promoted the expression of SAUR32, which activated proton extrusion related to Fe-deficiency response [26].
Circular RNAs (circRNAs) are also non-coding RNAs with covalently closed circular structures, which makes them more stable than linear RNAs [27]. CircRNAs are rich in microRNA binding sites and can absorb cytoplasmic miRNA. Thus, they relieve the inhibition caused by miRNA and therefore increase gene expression [28,29]. High-throughput sequencing and bioinformatics analyses have shown that circRNAs have potent regulatory roles in plants [30]. In wheat, 62 circRNAs related to plant growth and survival were differentially expressed under dehydration [31]. In tomato (solanum lycopersicum L.), 163 circRNAs responded to low temperature stress [32]. In Populus Euphratica Oliv, 18 circRNAs participated in the stress response [33]. In Arabidopsis, overexpression of a circRNA derived from glycerol-3-p acyltransferase (ATS1), Vv-circATS1, improved cold tolerance [34]. Thus, circRNAs are involved in adaptive responses to abiotic stresses, but their role in the regulation of adaptive responses remains to be confirmed.
The competitive endogenous RNA (ceRNA) regulation hypothesis states that mRNA, lncRNA, circRNA and pseudogene transcripts modulate the stability or translational activity of target genes via competitive binding with miRNA to achieve post-transcriptional gene regulation [35]. The ceRNA mechanism invokes miRNA as the vehicle to mutually regulate RNAs, making coding and non-coding genes part of a large and sophisticated regulatory network within the whole transcriptome. From its initial applications in oncology, progress in third generation sequencing technology and bioinformatics analysis methods has expanded ceRNA research to other fields. For example, it is involved in the adaptation to low nitrogen in poplar [36], in the response to copper toxicity in Citrus junos [37] and in the response to heat stress in Cucumis sativus L. [38]. These results greatly enhance our understanding of the mechanisms of abiotic stress response in plants.
Sugar beet (Beta vulgaris ssp. vulgaris) is one of the most important sugar-yielding crops worldwide. As a recently domesticated crop, cultivated beets have inherited salt tolerance traits from their wild ancestor Beta vulgaris ssp. maritima (referred to as B. maritima or ‘sea beet’) [39]. However, research on sugar beet ncRNA is limited. Herein, we have used whole-transcriptome RNA sequencing of cultivar O68, a widely used sugar beet in China, to uncover the global molecular response to salt stress at both protein-coding RNAs (mRNAs) and non-coding levels (lncRNAs, miRNAs, and circRNAs). Degradome sequencing was also performed to validate the ceRNA networks. In summary, our objective was to provide a basis for studying the potential function of ceRNA in response to salt stress in sugar beet by (I) identifying how mRNA and ncRNA are involved in this response and (II) dissecting the underlying ceRNA regulatory networks.

2. Results

2.1. Effects of Salinity on Sugar Beet Physiological Indices

We have shown previously that the relative germination rate of cultivar ‘O68′ under salt treatment was more than 90% and 70% in 200 and 300 mmol·L−1 NaCl, respectively [40], whereas seedling growth in the former case was even stronger than in the control group [40]. In the present study, seedlings exposed to 300 mmol·L−1 NaCl for 12 h showed morphological and physiological changes in both leaves and roots, such as clear wilting of leaves and color darkening of roots. We tested the tolerance to NaCl using physiological indices (Table 1), ‘ck’ represents the control group and ‘st’ represents the treatment group in this study. Salt stress reduced the relative water content and chlorophyll content in leaves. The soluble sugar content increased significantly, and the increase in leaves was greater than that of roots. Malondialdehyde content was slightly increased under salt stress in both leaves and roots. The content of proline in leaves increased significantly but not in roots. Moreover, the activity of peroxidase and superoxide dismutase was significantly enhanced under salt stress in both leaves and roots. The activity of catalase in leaves was significantly higher than that in roots and remained high under salt stress.

2.2. Global Response of mRNA to Salt Stress

A total of twelve transcriptome libraries were constructed from leaves and roots of six-leaf-stage seedlings which were either untreated (control) or treated with 300 mM NaCl. More than 1.5 billion raw reads were generated from RNA-seq with about 0.12 billion reads per sample. After quality control, 93.7% (1.4 billion) of the reads representing valid data were processed for further analysis. These reads were mapped on the beet genome with over 82.3% average mapping ratio. After quantification of gene expression levels as fragments per kilobase of transcript per million mapped reads (FPKM), 39,590 genes were identified as expressed genes. A Violin Plot shows the distribution of gene expression in each sample (Figure 1a). The image data indicated that gene expression levels were generally comparable among different time-point groups and these showed good repeatability. In leaves, 3578 differentially expressed mRNAs (DEmRNAs) (1180 up-regulated and 2398 down-regulated) were identified. A larger number, 4553 (1789 up-regulated, 2764 down-regulated), was identified in roots (Figure 1b and Table S1), which suggests a more robust response to salt stress. Although some DEmRNAs expression was tissue specific, with 2436 and 3363 in leaves and roots, respectively (Figure 1c), 1084 were common to both (Figure 1c), which likely form the basic response to salt stress.
We used gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment to explore DEmRNAs function (Figure 1d–g). GO annotation showed that the oxidation-reduction process was the dominant biological process (BP) in both leaves and roots (Figure 1d,e). However, in leaves, more DEmRNAs were annotated to response to cold and water deprivation, and flavonoid biosynthetic process under BP (Figure 1d). In contrast, in roots, those more abundant were defense response, signal transduction and response to abscisic acid (Figure 1e). For cellular component (CC), maximum DEmRNAs were annotated to chloroplast in leaves but plasma membrane in roots. In addition, a large amount of DEmRNAs belonging to the extracellular region and cell wall was found in both roots and leaves. In leaves, most DEmRNAs under molecular function (MF) were annotated to structural constituent of ribosome, RNA binding and endonuclease activity, while in roots were annotated to sequence-specific DNA binding, hydrolase activity and peroxidase activity.
KEGG enrichment analysis showed enrichment in Ribosome (ko03010), Starch and sucrose metabolism (ko00500) and Isoflavonoid biosynthesis (ko00943) in leaves (Figure 1f), whereas Phenylpropanoid biosynthesis (ko00940), Plant hormone signal transduction (ko04075) and MAPK signaling pathway-plant (ko04016) were significantly enriched in roots (Figure 1g).
Identification of metabolic pathways affected by salt stress, and comparison of different organ responses, was performed with a MapMan analysis (Figure 1h-i and Figure S1). In general, salt stress led to up-regulation of DEmRNAs involved in glycolysis and nucleotide metabolism (degradation) in both leaves and roots (Figure 1h,i). In leaves, up-regulation was observed in DEmRNAs associated with starch and sucrose synthesis, amino acid catabolism, and some secondary metabolites (such as flavonoids and phenols), whereas those related to cell wall, lipid metabolism (FA synthesis), photosynthesis and amino acid synthesis were down-regulated (Figure 1h). In roots, in contrast, those related to cell wall and amino acid synthesis were up-regulated, whereas most related to secondary metabolism were down-regulated (Figure 1i).
Regulation, cellular response, large enzyme families and transcription were used to evaluate metabolism-related DEmRNAs in detail (Figure S1). In both leaves and roots, metabolism for ABA (abscisic acid), ethylene and cytokinin (Figure S1a,e), heat stress and Redox (Figure S1b,f) were the most significantly up-regulated pathways. In leaves, many members of the cytochrome P450 family, peroxidases, UDP glycosyltransferases and the Glutathione-S-transferases family were up-regulated (Figure S1c). In leaves, transcription analysis showed that members of the AP2-EREBP and HSF families were up-regulated in leaves (Figure S1d). In roots, a large number of AP2-EREBP, bZIP, MYB and HB family members were significantly up-regulated (Figure S1h). Thus, these metabolism pathways defined the differential response between leaves and roots to salt stress, even under the same genetic background.

2.3. Global Response of lncRNA to Salt Stress

Apart from mRNAs, CNCI and CPC analysis identified 9076 potential lncRNAs from the RNA-seq data (Table S2). Since no sugar beet lncRNA was included, all lncRNAs identified here are novel lncRNAs. When mRNA and lncRNA were compared, more than 60% of the mRNA was more than 1000 bp long, whereas more than 60% of the lncRNA was less than 500 bp long. Also, one to three exons were observed for a higher percentage of lncRNAs, which also had a shorter ORF length and lower FPKM value (Figure 2a–d). Finally, the number of DElncRNAs in leaves was 66 (32 up-regulated, 34 down-regulated), and this number was much higher in roots, 453 (218 up-regulated, 235 down-regulated) (Figure 2e–g), suggesting a more robust response to salt stress at the lncRNA level (Figure 2e). The general expression profiles of DElncRNAs (Figure 2f,g) show a similar proportion of up- and down-regulated.
To explore the potential functions of these DElncRNAs, their targeted mRNAs were predicted by cis-regulated analysis, 100 K upstream and downstream. In leaves, no suitable targeted mRNAs were found, but 55 were identified for 61 DElncRNAs in roots (Table S3). The enriched GO terms of these targets are shown in Figure 2h, whereas enriched pathways indicated by KEGG enrichment analysis (Figure 2i) highlight targets such as Sulfur metabolism (ko00920), Pyrimidine metabolism (ko00240), Purine metabolism (ko00230) and Base excision repair (ko03410).

2.4. Global Responses of circRNA to Salt Stress

In the leaf and root, 2625 potential circRNAs were identified (Table S4), which showed no chromosome distribution preference (Figure 3a). Among these, 4.38, 10.52, and 85.10% belong to the intergenic region type, intron type and exon type, respectively (Figure 3b). In leaves, 13 DEcircRNAs were identified (eight up-regulated, five down-regulated), whereas 30 were found in roots (17 up-regulated, 13 down-regulated) (Figure 3d and Table S4) and one was common to leaves and roots (Figure 3c). More DEcircRNAs were highly expressed in the treatment group, compared to leaves (Figure 4e) and roots (Figure 4f), suggesting that circRNAs are also involved in the response to salt stress.

2.5. Global Responses of microRNA to Salt Stress

For miRNAs sequencing, we constructed twelve small RNA libraries from six-leaf-stage seedlings from four groups (treatment and control in roots and leaves) with three replicates. A total of 150 million raw reads with an average of 12.5 million per sample were generated. After filtering low-quality reads and removing sequences belonging to mRNA, repeat sequences and Rfam RNA tags, about 60 million reads with lengths 18 to 25 nt were selected as valid reads for further analysis. Finally, 494 mature miRNA candidates, including 165 known miRNAs (leaf: 153 and root: 130) and 329 novel miRNAs (leaf: 312 and root: 216) were identified in our small RNA sequencing data (Table S5). These mature miRNAs can be mapped to 487 pre-miRNAs. In Dhom’s genome sequencing study for B. vulgaris, 522 ncRNA sequences were classified as pre-miRNAs [41], whereas 161 mature miRNAs from our small RNA sequencing can be mapped to these sequences (Table S5). Those 21 nt long (37%) and 24 nt long (48%) account for the majority of candidate miRNAs, in agreement with the pattern of DICER-LIKE (DCL) cleavage in plant [42] (Figure 4a). Nucleotide composition of mature miRNAs indicated that both high-abundance miRNA clusters (21 nt or 24 nt) have a distinct sequence bias; 5′-uridine predominates in the first, and 5′-adenosine in the second (Figure 4b). This also has been observed in other plant species and is in agreement with the specificity of DCL restriction sites [42]. The distribution of identified miRNAs among different samples is shown in Figure 4c.
A total of 122 miRNAs showed significant changes in response to salt stress (Table S6), 73 in leaves (20 up-regulated, 53 down-regulated) and 64 in roots (20 up-regulated, 44 down-regulated). However, although the DEmiRNAs expression seems tissue specific, some of these DEmiRNAs may belong to a same family. For instance, three MIR166 family members, mtr-miR166e-5p, mtr-miR166e-5p_2ss3AT10GT and mtr-miR166e-5p_L-2R+2, were differentially expressed in leaves, whereas mtr-miR166a and mtr-miR166e-5p did so in roots. The general expression profiles of these DEmiRNAs (Figure 4d,e) show obvious differences between ck and salt-treated samples and between root and leaf samples. About one-third were up-regulated in both leaves and roots under salt stress, whereas two-thirds were down-regulated.
Unlike in animals, plant miRNAs have usually a perfect or near-perfect complementarity with their targets to induce target gene splicing and thus regulate gene expression; this feature can be used to predict plant miRNA targets. A total of 4150 target transcripts were identified for 432 miRNAs using the psRNATarget server and TargetFinder (Table S7). About 87% of miRNAs obtained had more than one target (Table S7). For example, a novel miRNA, PC-5P-3682_979, had 156 targets (Table S7). This high number of target genes supports the conclusion that a single miRNA has the ability to cleave multiple targets. The GO annotation of DEmiRNAs shows the main targets in leaves (Figure 4f) and in roots (Figure 4g). The significantly enriched GO terms of these targets in leaves an in roots are shown in Figure S2a,b, respectively. KEGG enrichment analysis showed that the significantly enriched pathways of these targets in leaves were Tyrosine metabolism (ko00350), Isoquinoline alkaloid biosynthesis (ko00950) and Sphingolipid metabolism (ko00600) (Figure S2c), and Glyoxylate and dicarboxylate metabolism (ko00630), Carotenoid biosynthesis (ko00906) and Cysteine and methionine metabolism (ko00270) in roots (Figure S2d).

2.6. CeRNA Regulatory Network in Response to Salt Stress

To reveal the global regulatory network of mRNAs and ncRNAs under salt stress, a ceRNA network was constructed using DEmiRNAs, DEmRNAs, DElncRNAs, and DEcircRNAs based on ceRNA theory. After interaction analysis and correlation filtered, 15 DEmRNAs, 23 DElncRNAs and three DEcircRNAs were predicted to be targets of nine DEmiRNAs. These filtered genes were used to construct the salt stress response ceRNA network in leaves and roots (Figure 5). P threshold can be temporarily relaxed in some specific cases, such as the one-sided significance level was set at 0.2 in the study of Erlotinib alone or with Bevacizumab as first-line therapy in patients with lung cancer [43]. For leaves, since no ceRNAs had been obtained in leaves, the criteria of DEmRNA and DEncRNA were temporarily relaxed from p < 0.05 to p < 0.1 to explore potential ceRNAs, then qPCR was used to verify the expression correlation of all ceRNAs in subsequent experiments. Four novel miRNAs, PC-3p-125293_44, PC-5p-109182_51, PC-5p-7955_523 and PC-5p-3682_979, participated in the network, with the latter having an important role in both leaf and root under salt stress. The coding DEmRNAs of the ceRNA network in leaves included basic blue protein (Bv1_023200_jmkt.t1), uclacyanin-2 (Bv6_144730_qgsa.t1), stress protein DDR48 (Bv1_011770_noyh.t1), G-type lectin S-receptor-like serine/threonine protein kinase LECRK3 (Bv4_071790_jfup.t1) and tocopherol cyclase (Bv2_032060_rrem.t1). While the coding DEmRNAs of ceRNA network in roots included tubulin beta-1 chain (Bv5_100510_ttjr.t1), acid beta-fructofuranosidase (Bv5_097930_juac.t1), hypothetical protein LOC104907984 (Bv5_124210_skcr.t1), NAC domain-containing protein 21/22 (Bv_000470_dyzk.t1), NAC domain-containing protein 100 (Bv5_114390_pjnp.t1), uncharacterized protein CENP-C (Bv9_220500_yszy.t2), inositol hexakisphosphate and diphosphoinositol-pentakisphosphate kinase 2 (Bv2_040500_uums.t1), hypothetical protein BVRB_5g117790 (Bv5_117790_mxei.t1), malate synthase (Bv5_110660_jxri.t1), and hypothetical protein BVRB_9g213960 (Bv9_213960_qddk.t1).

2.7. Verification of the Cleavage of miRNA to ceRNAs by Degradome Sequencing

Degradome sequencing was used to verify the cleavage of targets by miRNA. For leaves and roots, a respective degradome library was constructed and sequenced. A total of 51.4 million raw reads were generated. The clean tags were mapped to sugar beet reference genome and transcriptome data for identification of cleavage sites after consecutive steps of filtering. Based on the abundance of degradome tags at the target sites, cleaved targets were classified into five categories (0 to 4). The p-value here is not a screening parameter but for reference only; most p-values of Category 0 and Category 1 were less than 0.05 (Table S8). A total of 752 miRNA–mRNA pairs were identified for 141 miRNAs and 444 transcripts (including 406 mRNAs and 38 lncRNAs) and six pairs of miRNA–mRNA, which were used to build the ceRNA network. Some representative cleavage sites are shown as target plots (T-plots) (Figure 6). Interestingly, degradome sequencing confirmed that a single miRNA can act on multiple targets, e.g., mtr-miR408-3p_L-1R+1 can cleave both basic blue protein (bv1_023200_jmkt) and uclacyanin-2 (bv6_144730_qgsa) (Figure 6b,c). In addition, different miRNAs from the same family can work together to regulate the same genes, e.g., basic blue protein (bv1_023200_jmkt) was cleaved by gma-miR408a-3p_L-1R+5, mtr-miR408-3p_L-1R+1 and mtr-MIR408-p3_2ss18GT19GA (Figure 6a,c,d). Combining the expression patterns of DEmRNAs and DEmiRNAs, 16 and 17 negative correlations in expression were identified in miRNA–mRNA pairs in leaves and roots, respectively, where 11 pairs were confirmed by degradome sequencing (Table 2).

2.8. qRT-PCR Validated Expression Correlation between miRNAs and ceRNAs under Salt Stress

To confirm the results of RNA-seq and to validate the expression correlation of miRNAs and their targets, all components of the ceRNA network were selected for qRT-PCR detection (Figure 7). The qRT-PCR results agree well with RNA-seq data regarding expression trends. However, qPCR fold-change of some genes was different from the sequencing data, especially for ncRNA, possibly because the different principles underlying these two detection methods. The expression of all miRNAs was negatively correlated with ceRNAs, whereas those of ceRNAs were positively correlated. For example, gma-miR6300_L-1R+1 was up-regulated and all of its targets were down-regulated; PC-5p-109182_51 was down-regulated and all of its predicted targets were up-regulated. This result not only demonstrates that the RNA-seq data in this study is reliable, but also validates the negatively correlated expression between miRNAs and ceRNAs. No expression was detected for two down-regulated lncRNAs in the ST group (MSTRG.30182.1 and MSTRG.28802.1), possibly because of low abundance.

3. Discussion

Salt stress threatens the normal growth and development of plants by directly inducing osmotic imbalance, ion injury and reactive oxygen species (ROS). Thus, response to salt stress involves a variety of physiological and biochemical processes, such as osmotic regulation, ion balance or reactive oxygen scavenging. These are regulated by the expression of a large number of genes. Although miRNAs, lncRNAs and circRNAs have been identified in plants, there are no studies of ncRNAs in sugar beet. In this study, DEmRNAs, DElncRNAs, DEcircRNAs and DEmiRNAs were identified in leaves and roots under salt stress. These results support the ceRNA hypothesis which has progressively wider acceptance since it was first proposed [35]. These RNAs compete for binding to common MREs and thus are mutually regulated, with implications for posttranscriptional gene regulation in physiological and pathophysiological processes [44]. From this perspective, we used our sequencing data in sugar beet to construct a ceRNA regulatory network to elucidate the mechanisms underlying tolerance to salt. The fact that the number of salt stress-response genes (including mRNA and ncRNA) in roots was higher than that in leaves, and the functional divergence observed in these two organs, suggests that roots have a more robust and complex response to salt.

3.1. Analysis of Salt Stress Response in Sugar Beet Leaves

Functional analysis of differentially expressed genes (DEGs) in leaves showed that leaf cells appear to regulate gene expression with the maintenance of photosynthesis as the core. Some genes directly involved in photosynthesis were significantly up-regulated under salt stress. For example, NAD(P)H-quinone oxidoreductase chain 2 (Bv7_171100_eawo.t1) involved in NADH dehydrogenase (ubiquinone) activity, proton gradient regulation protein 5 (Bv3_048340_odem.t1) involved in photosynthetic electron transport in photosystem I or STN7 protein kinase (Bv6_155590_ugjw.t1) involved in the regulation of light reaction. In addition, some genes located in the chloroplast that play a protective role in photosynthesis were also significantly up-regulated under salt stress. For example, choline monooxygenase (Bv6_146100_wdro.t1) and betaine aldehyde dehydrogenase (Bv5_116230_ntjn.t1), which are involved in the synthesis of betaine [45]; glutathione peroxidase (Bv2_044150_kdwp.t1), ferrochelatase 2 (Bv2_033330_nrat.t1) and soul heme-binding family protein (Bv4_094190_wifn.t1) involved in heme biosynthesis in non-photosynthetic tissues and induced by oxidative stress in photosynthetic tissues to supply heme for defensive hemoproteins. Studies of rhodiola crenulata have shown that glutathione peroxidase can physically interact with specific proteins in multiple pathways which participate in photosynthesis, respiration and stress tolerance [46]. Further, alpha amylase (Bv1_021820_ccxs.t1) was significantly up-regulated under salt stress, which suggests faster starch degradation by mesophyll cells. This may increase soluble sugar content and provide energy to regulate osmotic pressure and respond to salt stress. The increase in soluble sugar content in leaves also supports this hypothesis (Table 1).
Copper ion plays an irreplaceable role in chlorophyll synthesis, stabilization and photosynthesis; thus, regulating metal homeostasis in chloroplasts is extremely important [47]. Plants may have to redistribute a limited copper ion population in leaves, since salt stress results in reduced transpiration and hindered transport of copper ions from the roots. Indeed, we found that basic blue protein (Bv1_023200_jmkt.t1) and uclacyanin-2 (Bv6_144730_qgsa.t1) were down-regulated by the mtr-MIR408-p3_2ss18GT19GA-mediated ceRNA network (Figure 8), which may be related to copper distribution. This is in agreement with our previous studies on the proteome of sugar beet under salt stress revealing up-regulation of plastocyanin in leaves [48], consistent with copper ions being preferentially supplied to photosynthesis. Other studies have shown that the expression of miR408 was up-regulated under abiotic stress in Arabidopsis thaliana [49] and that the overexpression of miR408 enhanced drought tolerance in chickpea (Cicer arietinum L.) [50].
Under salt stress, increased membrane permeability and disruption of Na+/K+ homeostasis is enabled by α-tocopherol (α-T) by increased unsaturation of the membrane [51], where tocopherol cyclase may play an important role. Our results show that PC-5p-3682_979-mediated ceRNAs up-regulated the expression of tocopherol cyclase (Bv2_032060_rrem.t1) (Figure 8). Under salt stress, the content of tocopherol was significantly decreased in beet leaves, while the content of MDA did not change significantly (Table 1). We hypothesized that sugar beets protect their membrane lipids from peroxidation by consuming tocopherol. The same Figure shows that the ceRNA network revealed up-regulation of ABA induced stress protein DDR48 (Bv1_011770_noyh.t1) and receptor-like protein kinase 1 (RPK1) (Bv4_071790_jfup.t1) located in the endomembrane system; these may be involved in the response to salt stress through ABA signaling and phosphorylation pathways. In Arabidopsis thaliana, RPK1 is involved in regulating the shutdown of guard cells that is key for plant stress response [52].
Moreover, a number of genes associated with ABA and ethylene metabolism and signal transduction were up-regulated in leaves: protein phosphatase 2C (Bv8_197440_jnqo.t1), zeaxanthin epoxidase (Bv_003210_fkxn.t1), carotenoid cleavage dioxygenase (Bv5_107230_qtix.t1) and ERF3 (Bv2_031130_cghd.t1), 1-aminocyclopropane-1-carboxylate oxidase (Bv8_185300_wndh.t1). Indeed, Nicotinamide adenine dinucleotide (NAD) responds to salt stress by affecting the accumulation of ABA and proline [53] and up-regulation of NAD(P)H-quinone oxidoreductase and ABA-related genes and the increase in proline content support this conclusion. 1-aminocyclopropane-1-carboxylate oxidase (ACO) is a key enzyme in the ethylene biosynthesis pathway. The more than 20-fold change of carotenoid cleavage dioxygenase has received the most attention, suggesting that it may play an important role in the response to salt stress. These results suggest that leaves are regulated to maintain photosynthesis under salt stress (Figure 8).

3.2. Analysis of Salt Stress Response in Sugar Beet Roots

In contrast to what we observed in leaves, roots showed a broader mobilization under salt stress. Functional analysis of DEGs in roots indicated that root cells shift resources from growth and development to survival by adjusting several processes, e.g., glycometabolism and metabolism of fatty acids and energy and secondary metabolism. Sucrose is the end product of photosynthesis, yielding hexoses (Hexes) necessary to generate energy and synthesize diverse biomolecules. A significant up-regulation of sucrose synthases (Bv7_163460_jmqz.t1 and Bv8_190960_nnjy.t2) was detected in root cells under salt stress. Sucrose metabolism has been associated with high sensitivity of development to abiotic stress in plants, since a reduction in hexose triggers a downstream response to stress [54]. Also, pressure can affect the expression of acid beta-fructofuranosidase, involved in the metabolic signaling of primary metabolism and defense response [55,56]. Our results show that gma-miR6300_L-1R+1-mediated ceRNAs down-regulate the expression of Bv5_100510_ttjr.t1 and Bv5_097930_juac.t1. Bv5_100510_ttjr.t1 (Figure 8). These genes encode tubulin beta-1 chain protein mostly expressed in endodermal and phloem cells of primary roots. The role of microtubules as signal transducers for plant stress response has been reported [57]. Bv5_097930_juac.t1 encodes acid beta-fructofuranosidase which participates in the degradation of sucrose and is associated to hexose signals. Interestingly, a series of glycolytic-related genes were up-regulated under salt stress in roots, such as hexokinase (Bv9_224680_aote.t1), phosphofructokinase (bv9_212010_mhum.t1) and pyruvate kinase (Bv7_160680_tjji.t1 and Bv2_034650_iitw.t1), which are the rate-limiting enzymes in glycolysis (Figure 8). Pyruvate dehydrogenase (Bv_004020_oxcs.t1) and NADP-dependent malate dehydrogenase (Bv8_197290_sfck.t1 and Bv_006710_gkqg.t1) were also up-regulated in roots. Further, PC-5p-109182_51-mediated ceRNAs up-regulate the expression of Bv5_110660_jxri.t1 and Bv9_213960_qddk.t1. Bv5_110660_jxri.t1 encodes malate synthase, involved in the conversion of fatty acids to carbohydrates and important in maintaining gluconeogenesis [58]. Expression of phosphoenolpyruvate carboxykinase (Bv6_152970_ygjo.t1) was reduced by 20-fold under salt stress. Thus, salt stress may weaken the gluconeogenic pathway, and acetyl-coa and oxaloacetic acid is probably provided to the tricarboxylic acid cycle by other pathways. These results suggest that salt stress increases energy metabolism in root cells.
NAC domain-containing proteins are one of the largest plant-specific transcription factor families. Proteins in this family are characterized by a highly-conserved N-terminus, known as the DNA-binding NAC domain (NACs), and a highly divergent C-terminal region containing the transcription regulatory region. Besides their critical functions in plant growth and development [59], NACs have also been implicated in the response to abiotic and biotic stresses [60]. We show here that two NAC family members, NAC domain-containing protein 21/22 (Bv_000470_dyzk.t1) and NAC domain-containing protein 100 (Bv5_114390_pjnp.t1), are down-regulated by mtr-miR164d-mediated ceRNAs under salt stress (Figure 8). We confirmed cleavage of these two genes by mtr-miR164d by degradation sequencing. This down-regulation may activate multiple downstream genes, which is in agreement with a Nicotiana benthamiana study infected by beet necrotic yellow vein virus [61]. PC-5P-3682_979 regulates the expression of tocopherol cyclase in leaves, where it may perform a different function. In roots, a PC-5P-3682_979-mediated ceRNA seems to be involved in the phosphoinositol signaling system since it down-regulated inositol hexakisphosphate and diphosphoinositol-pentakisphosphate kinase 2 (Bv2_040500_uums.t1) and F-box/kelch-repeat protein At3g06240-like (Bv5_117790_mxei.t1) (Figure 8).
In addition, the ceRNA network also revealed the down-regulation of two functional unknown proteins uncharacterized protein CENP-C (Bv9_220500_yszy.t2) and hypothetical protein LOC104907984 (Bv5_124210_skcr.t1), which may play roles in response to salt stress by unknown pathway in sugar beet (Figure 8). Additionally, more remarkable, many AP2-EREBP, bZIP and MYB family transcription factors were significantly up-regulated under salt stress, such as ERF73 (Bv7_163350_ryod.t1), DREBP1 (Bv3_066590_ignp.t1), TGA3 (Bv1_021380_gmre.t1), ABSCISIC ACID-INSENSITIVE 5 (Bv7_159570_afnu.t1), MYB39 (Bv8_181400_zguf.t1), etc. It has been reported that the DREB transcription factor is involved in the salt stress response of tomatoes and arabidopsis thaliana [62] and that the overexpression of DREB transcription factors can increase the salt tolerance of potatoes (solanum tuberosum) [63]. These transcription factors may play a crucial regulation role in the response to salt stress in roots. Similar to leaves, a large number of genes related to ABA and ethylene metabolism were also significantly upregulated in roots under salt stress (Figure 8), for example, PP2C 8 (Bv3_052290_rcph.t1) and PP2C 24 (Bv7_173230_chqm.t1), 1-aminocyclopropane-1-carboxylate synthase (Bv1_004040_gcto.t1), and so on. Further, peroxidase 5-like (Bv4_079730_uxri.t1), glutathione S-transferase (Bv4_095800_rgfa.t1), and some cell wall-related genes such as pectinesterase 11 (Bv1_006810_mswa.t1) and pectin acetylesterase 8 (Bv5_098360_owsf.t1) were also significantly up-regulated under salt stress (Figure 8). The above results support the idea that an extensive metabolic regulation was carried out in roots in favor of cell survival under salt stress (Figure 8).
Overall, these results indicate that miRNAs–lncRNAs-circRNAs–mRNAs networks are involved in regulating gene expression to modify growth, improve photosynthesis, adjust plasma membrane permeability, promote glycometabolism and energy metabolism, regulate transcription factor, and participate in the phosphoinositol signaling system during adaptation to salt stress.

4. Materials and Methods

4.1. Plant Cultivation and Treatments

Sugar beet ‘O68′ was used in this study. Seeds from our laboratory (Heilongjiang, China) were soaked in water for ten hours, sterilized in 0.1% (v/v) HgCl2 for 10 min, washed repeatedly with distilled water and germinated on wet filter paper in a germination box at 26 °C for 2 days. After germination, seedlings were transferred to plastic pots (43.5 cm × 20 cm × 14 cm, 10 plants per pot) filled with quarter-strength Hoagland solution. The germinating seeds were cultivated under 16/8 light photoperiod at 24 °C (day)/18 °C (night) in a phytotron (Friocell 707, Germany). Six-leaf-stage seedlings (28 days) were treated with 300 mmol NaCl for 12 h and untreated plants served as controls.

4.2. Measurement of Physiologic Indicators and Harvest

Fresh samples were used to measure physiologic indicators. All of the physiological indicators listed in Table 1 were carried out in accordance with the Experimental Guidelines for Plant Physiology [64]. Samples were collected from leaves and roots (about 6 cm from the root tip) as shown in Figure S3 for sequencing. Variability was minimized by constructing the library with samples taken from six different plants. All collected samples were immediately frozen in liquid nitrogen for half an hour and stored at −80 °C until further use.

4.3. RNA Extraction, Library Preparation, and RNA Sequencing

A total of 0.5 g of leaves, or roots, was taken from six plants and mixed to obtain each sample. Total RNA was extracted from the roots and leaves of ck and st samples by TRIZOL (Invitrogen, CA, USA) according to the manufacturer’s instructions. Quantity and purity of total RNA were analyzed using Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number > 7.0. High-quality total RNA from each sample was used to construct libraries for RNA-seq and small RNA sequencing. Total RNA from leaves and roots was mixed for degradome sequencing. Library construction and sequencing were conducted by LC Sciences (Hangzhou, China).
For mRNA, LncRNA and circRNA sequencing, about 5 µg of total RNA from each sample was used to deplete ribosomal RNA according to the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, USA). After removing ribosomal RNAs, RNAs were fragmented by divalent cations under a high temperature. Cleaved fragments were reverse-transcribed to create the cDNA, which was used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I, RNase H and dUTP. There were three biological replicates per group and a total of 12 libraries (leaf_ck, leaf_st, root_ck, root_st). The average insert size in the final cDNA library was 300 bp (±50 bp). Finally, we performed paired-end sequencing on an Illumina Novaseq™ 6000 (LC Bio, China) following the vendor’s recommended protocol. For small RNA sequencing, 12 small RNA libraries (leaf_ck, leaf_st, root_ck, root_st) were constructed according to the TruSeq Small RNA Sample Preparation (Illumina, San Diego, CA). The small RNA fragments were finally sequenced by Illumina HiSeq2500 platforms. The raw data of RNA-seq (BioProject ID: PRJNA666384) and small RNA sequencing (BioProject ID: PRJNA666142) were deposited in the NCBI SRA database.

4.4. Read Mapping and Transcriptome Assembly

For RNA-Seq data, the raw reads obtained from sequencing were processed to remove primer/adaptor contamination, low-quality bases and undetermined bases using Cutadapt [65]. FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to verify sequence quality. Bowtie2 [66] and Hisat2 [67] were used to map reads to the reference genome RefBeet-1.2 (http://bvseq.boku.ac.at/Genome/Download/RefBeet-1.2/). The mapped reads of each sample were assembled using StringTie [68].

4.5. Differentially Expressed mRNA and Bioinformatics Analysis

All transcriptomes from all samples were merged to reconstruct a comprehensive transcriptome using Perl scripts, after which StringTie [68] and edgeR [69] estimated the expression levels of all transcripts. StringTie was used to perform the expression level of mRNAs by calculating FPKM (fragments per kilobase of transcript per million mapped reads). The differential expression analysis was performed using edgeR R-packages. The p-value < 0.05 and |log2 (fold change)| > 1 were defined as the cutoff criteria for DEGs. UpSetR was used to visualize differences in gene expression patterns between groups [70]. By default, 0.000001 was added to all FPKM values to avoid computing log (0).
All transcripts were annotated by conducting Blastx similarity searches against NCBI non-redundant protein databases (Nr) with an E-value threshold of 10−5. All transcripts were blasted to the Gene Ontology database and calculated gene numbers for each term. A hypergeometric test was used to find significantly enriched GO terms in the input gene list [71]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) [72] was used to perform pathway enrichment analysis. Mapman software was used to analyze DEGs (http://mapman.gabip.org, version: 3.6.0).

4.6. Identification and Analysis of lncRNAs

Transcripts that overlapped with known mRNAs, and transcripts shorter than 200 bp, were discarded during the identification of novel lncRNAs. We used CPC [6] and CNCI [7] to predict transcripts with coding potential. All transcripts with CPC score < 0.5 and CNCI score < 0 were removed and those remaining were considered to be lncRNAs. StringTie was used to perform expression level for mRNAs and lncRNAs by calculating FPKM. The differentially expressed lncRNAs (DE lncRNAs) were selected with |log2 (fold change)| > 1 and p-value < 0.05 by edgeR. To explore the function of lncRNAs, we predicted their cis-target genes. LncRNAs may play a cis role acting on neighboring target genes. In this study, coding genes in 100,000 upstream and downstream were selected by Python script. Then, we showed functional analysis of the target genes for lncRNAs by using the BLAST2GO [73]. Significance was expressed as a p-value < 0.05.

4.7. Identification and Analysis of CircRNAs

In this study, Bowtie2 [66] and Hisat2 [67] were used to map reads to the genome of sugar beet. The remaining reads (unmapped reads) were still mapped to genome using tophat-fusion [74]. CIRCExplorer2 [75,76] and CIRI [77] were used to de novo assemble the mapped reads to circular RNAs at first; then, back splicing reads were identified in unmapped reads by tophat-fusion [74]. All samples were generated unique circular RNAs. Differentially expressed circRNAs (DE circRNAs) were extracted with |log2 (fold change)| >1 and with statistical significance (p-value < 0.05) by R package–edgeR.

4.8. Identification and Analysis of miRNAs

The raw reads obtained from sequencing were subjected to ACGT101-miR (LC Sciences, Houston, Texas, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats. As no sugar beet miRNA was recorded in miRBase 21.0, the filtered unique reads with lengths between 18 and 25 nt were mapped to viridiplantae mature miRNAs and precursors by BLAST search to identify conserved and previously reported miRNAs, and the mapped pre-miRNAs were further aligned against the B. vulgaris genomes to determine their genomic locations using Bowtie [78]. Only miRNAs matched to known miRNAs with no more than three mismatches in the miRBase database, and whose precursors could fold into stem-loop structures, were considered to be known miRNAs. Subsequently, the unmapped sequences were blasted against the B. vulgaris genomes to identify potential novel miRNAs, and the hairpin structures-containing sequences were predicted from the flanking 200 nt sequences using the mfold software [79]. A global normalization method was used to correct copy numbers among different samples [80]. The latter method can eliminate the effect of sequencing discrepancy on the calculation of small RNA expression. Therefore, the calculated gene expression can be directly used to compare the difference of gene expression among samples. Only miRNAs with Fold Change 1.5 and p-value < 0.05 were classified as DE miRNAs.

4.9. Degradome Sequencing Validation of the Cleavage of miRNA to Target Genes

For degradome library construction, about 20 μg of total RNA from both treatment and control groups were pooled together to generate a mixed-library for leaf and root. The library construction was performed as described with some modifications [81]. Single-end sequencing (36 bp) of the purified cDNA libraries was performed on an Illumina Hiseq2500 at the LC-BIO (Hangzhou, China) following the vendor’s recommended protocol. Raw sequencing reads were obtained using Illumina’s Pipeline v1.5 software following sequencing image analysis by Pipeline Firecrest Module and base-calling by Pipeline Bustard Module. The raw reads were processed to remove low-quality reads, reads with ‘N’ and any reads with adaptor and primer contamination using ACGT101-DEG (Lc-bio Sciences, Houston, USA). The clean reads were used to identify the degraded fragments of mRNA after removing the reads that can be annotated into a different ncRNA database. Potential miRNA editing sites were identified using the small RNA sequencing data by CleaveLand 3.0 [82]. The identified sites were classified into five categories (0, 1, 2, 3 and 4) based on read abundance at that position as reported previously [83]. Based on the signatures and abundance in the sugar beet RNA-seq data, T-plots were built for high efficiency analysis of potential miRNA targets.

4.10. Construction and Analysis of ceRNAs Regulatory Network

In plants, there are two patterns of complementarity between miRNA and their target gene mRNA; one leads to complete target degradation, whereas, in the other, miRNA is sequestered and therefore inactivated by miRNA target mimics [84]. Based on the ceRNA hypothesis, a ceRNA network can be built by predicting miRNA-binding RNA. In this study, miRNA–mRNA and miRNA–ncRNA pairs were predicted by TargetFinder (https://github.com/carringtonlab/TargetFinder) and Ssearch36 (36.3.6) (https://fasta.bioch.virginia.edu/fasta_www2/fasta_down.shtml), respectively. In brief, (I) the bulge must be on the ncRNA and located in the middle of the miRNA, (II) except for the position in the middle of miRNA, a maximum of 4 mismatches are allowed and there cannot be more than two continuous mismatches and (III) no bulge is allowed in a non-middle position. A ceRNA network was visualized with Cytoscape software (version 3.6.1) [85].

4.11. qRT-PCR Validation of DEmRNAs, DElncRNAs, and DEmiRNAs

qRT-PCR was performed to validate the expression level of DEmRNAs, DE lncRNAs, DE circRNAs and DE miRNAs. Total RNA and small RNA was extracted from leaves and roots samples using an GeneJET Plant RNA Purification Mini Kit (Thermo Scientific, Cat#K0801, Lithuania) and mirVana™ miRNA Isolation Kit (Invitrogen, Cat#K157001, USA), respectively. Firststrand cDNA was synthesized from 1 μg total RNA with a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Cat#4368814, USA) for qPCR of mRNA and miRNA, and with the Verso cDNA Synthesis Kit (Thermo Scientific, Cat#AB1453A, Lithuania) for qPCR of lncRNA and circRNA. The qPCR reactions were performed on CFX Real-time PCR system (BIO-RAD, USA) using iTaq Universal SYBR® Green Supermix (BIO-RAD, Cat#1725125, USA) following the manufacturer’s instructions. For miRNAs, the primers including miRNA-specific stem-loop RT, forward primers and universal reverse primer for the selected miRNAs were designed according to Kramer [86]. The 5S RNA + U6 gene was used as endogenous controls. For mRNA, lncRNA and circRNA, PP2A + UBQ5 and PP2A + 25S RNA genes were used as endogenous controls for root and leaf, which were selected from 15 candidate genes evaluated using Bestkeeper [87], NormFinder [88] and GeNorm [89]. To avoid non-specific amplification, a melting curve was carried out for each PCR product. The expression level of the miRNAs in different samples was calculated by the comparative 2−△△CT method. For the data of the physiological parameters and qPCR analysis, the mean and SD were calculated from three repeats per treatment, and differences were analyzed by Duncan’s multiple range test (p < 0.05) and an independent samples t-test (p < 0.05). All primers are shown in Table S9.

5. Conclusions

We have used whole-transcriptome RNA-seq to reveal the divergent response to salt stress in leaves and roots, identifying 3578 and 4553 DEmRNAs, 66 and 453 DElncRNAs, 73 and 64 DEmiRNAs and 13 and 30 DEcircRNAs, respectively. Our results suggest that specific lncRNAs and circRNAs function as ceRNAs in response to salt stress. GO term, KEGG pathway and Mapman analyses showed that competitive lncRNA/circRNA/mRNA-miRNA regulatory networks are implicated in copper redistribution, plasma membrane permeability, glycometabolism and energy metabolism, NAC transcription factor and the phosphoinositol signaling system. In a nutshell, the metabolic regulation of leaves revolves around ensuring photosynthesis to obtain the energy necessary to cope with environmental pressure whereas roots obtain energy by enhancing glycometabolism and fatty acid metabolism to cope with these environmental changes. Many transcription factors such as MYB, DREBP, NAC and plant hormones such as ABA and ethylene are also involved in the salt stress response. Overall, these results lay the foundation for studying lncRNAs/circRNAs response to salt stress in different organs of sugar beet.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/1422-0067/22/1/289/s1, Table S1: Details information of identified differentially expressed mRNAs, Table S2: Details information of identified lncRNAs, Table S3: Target Prediction of DElncRNA in roots, Table S4: Details information of identified circRNAs, Table S5: Details information of identified miRNAs (known and novel), Table S6: Details information of DEmiRNAs, Table S7: Target Prediction of miRNA by bioinformatics, Table S8: Target identified of miRNAs by degradome sequencing, Table S9: List of primers used for qRT-PCR experiments, Figure S1: Mapman analysis of DEmRNAs in leaves (a–d) and roots (e–h) under salt stress, Figure S2: GO and KEGG enrichment analysis of DEmRNAs targets, Figure S3: Sampling location.

Author Contributions

J.L.: Conceptualization, Data curation, Formal analysis, Methodology, Writing–original draft; J.C.: Conceptualization, Methodology, Project administration; C.D.: Supervision, Validation, Writing–review & editing; T.L.: Supervision, Validation; D.C.: Conceptualization, Project administration; C.L.: Investigation, Resources. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (31571731 and 31771864).

Acknowledgments

We thank the laboratory of exploration and utilization of functional genes in sugar beet for their participation in the study.

Conflicts of Interest

The authors declare no conflict of interest.

Ethics Approval and Consent to Participate

Not applicable.

Availability of Data and Materials

All sequencing data was deposited in the NCBI Short Read Archive (SRA) database under the BioProject ID: PRJNA666384 (RNA-seq) and BioProject ID: PRJNA666142 (miRNA and degradome sequencing). Relevant supporting data can be found within the article and additional files.

Abbreviations

DEmRNADifferentially expressed mRNA
DElncRNADifferentially expressed lncRNA
DEcircRNADifferentially expressed circRNA
DEmiRNADifferentially expressed miRNA
DEGDifferentially expressed gene
MDAmalonaldehyde
PODperoxidase
SODsuperoxide dismutase
CATcatalase
ABAabscisic acid
ROSreactive oxygen species
FPKMfragments per kilo-base per million reads
GOgene ontology
KEGGKyoto encyclopedia of genes and genomes
qRT-PCRquantitative real-time polymerase chain reaction

References

  1. Rengasamy, P. Soil processes affecting crop production in salt-affected soils. Funct. Plant Biol. 2010, 37, 613–620. [Google Scholar] [CrossRef]
  2. Munns, R.; Gilliham, M. Salinity tolerance of crops-what is the cost? New Phytol. 2015, 208, 668–673. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Tester, M.; Davenport, R. Na+ Tolerance and Na+ Transport in Higher Plants. Ann. Bot. 2003, 91, 503–527. [Google Scholar] [CrossRef] [PubMed]
  4. Munns, R.; Tester, M. Mechanisms of Salinity Tolerance. Annu. Rev. Plant Biol. 2008, 59, 651–681. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Flowers, T.J.; Colmer, T.D. Plant salt tolerance: Adaptations in halophytes. Ann. Bot. 2015, 115, 327–331. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Garthwaite, A.J.; Bothmer, R.V.; Colmer, T.D. Salt tolerance in wild Hordeum species is associated with restricted entry of Na+ and Cl- into the shoots. J. Exp. Bot. 2005, 56, 2365–2378. [Google Scholar] [CrossRef]
  7. Deinlein, U.; Stephan, A.B.; Horie, T.; Luo, W.; Xu, G.; Schroeder, J.I. Plant salt-tolerance mechanisms. Trends Plant Sci. 2014, 19, 371–379. [Google Scholar] [CrossRef] [Green Version]
  8. Feng, K.; Yu, J.; Cheng, Y.; Ruan, M.; Wang, R.; Ye, Q.; Zhou, G.; Li, Z.; Yao, Z.; Yang, Y.; et al. The SOD Gene Family in Tomato: Identification, Phylogenetic Relationships, and Expression Patterns. Front. Plant Sci. 2016, 7, 1279. [Google Scholar] [CrossRef] [Green Version]
  9. Kapilan, R.; Vaziri, M.; Zwiazek, J.J. Regulation of aquaporins in plants under stress. Biol. Res. 2018, 51, 4. [Google Scholar] [CrossRef]
  10. Liu, X.S.; Feng, S.J.; Zhang, B.Q.; Wang, M.Q.; Cao, H.W.; Rono, J.K.; Chen, X.; Yang, Z.M. OsZIP1 functions as a metal efflux transporter limiting excess zinc, copper and cadmium accumulation in rice. BMC Plant Biol. 2019, 19, 283. [Google Scholar] [CrossRef]
  11. Chekanova, J.A. Long non-coding RNAs and their functions in plants. Curr. Opin. Plant Biol. 2015, 27, 207–216. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Yu, Y.; Zhang, Y.; Chen, X.; Chen, Y. Plant Noncoding RNAs: Hidden Players in Development and Stress Responses. Annu. Rev. Cell Dev. Biol. 2019, 35, 407–431. [Google Scholar] [CrossRef] [PubMed]
  13. Bartel, D.P. MicroRNAs, Genomics, Biogenesis, Mechanism, and Function. Cell 2004, 116, 281–297. [Google Scholar] [CrossRef] [Green Version]
  14. Hu, J.; Jin, J.; Qian, Q.; Huang, K.; Ding, Y. Small RNA and degradome profiling reveals miRNA regulation in the seed germination of ancient eudicot Nelumbo. Nucifera. BMC Genom. 2016, 17, 684. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Chen, X. A MicroRNA as a translational repressor of APETALA2 in Arabidopsis flower development. Sci. N. Y. 2004, 303, 2022–2025. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Millar, A.A.; Gubler, F. The Arabidopsis GAMYB-like Genes, MYB33 and MYB65, are microRNA-regulated genes that redundantly facilitate anther development. Plant Cell 2005, 17, 705–721. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Yan, J.; Zhang, H.; Zheng, Y.; Ding, Y. Comparative expression profiling of miRNAs between the cytoplasmic male sterile line MeixiangA and its maintainer line MeixiangB during rice anther development. Planta 2015, 241, 109–123. [Google Scholar] [CrossRef] [PubMed]
  18. Sunkar, R.; Chinnusamy, V.; Zhu, J.; Zhu, J.-K. Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 2007, 12, 301–309. [Google Scholar] [CrossRef] [PubMed]
  19. Zhang, T.; Zhao, Y.-L.; Zhao, J.-H.; Wang1, S.; Jin1, Y.; Chen, Z.-Q.; Fang, Y.-Y.; Hua, C.-L.; Ding, S.-W.; Guo, H.-S. Cotton plants export microRNAs to inhibit virulence gene expression in a fungal pathogen. Nat. Plants 2016, 2, 16153. [Google Scholar] [CrossRef] [PubMed]
  20. Liu, J.; Wang, H.; Chua, N.-H. Long noncoding RNA transcriptome of plants. Plant Biotechnol. J. 2015, 13, 319–328. [Google Scholar] [CrossRef] [PubMed]
  21. Chen, L.-L. Linking Long Noncoding RNA Localization and Function. Trends Biochem. Sci. 2016, 41, 761–772. [Google Scholar] [CrossRef] [PubMed]
  22. Shuai, P.; Liang, D.; Tang, S.; Zhang, Z.; Ye, C.-Y.; Su, Y.; Xia, X.; Yin, W. Genome-wide identification and functional prediction of novel and drought-responsive lincRNAs in Populus Trichocarpa. J. Exp. Bot. 2014, 65, 4975–4983. [Google Scholar] [CrossRef] [PubMed]
  23. Deng, F.; Zhang, X.; Wang, W.; Yuan, R.; Shen, F. Identification of Gossypium hirsutum long non-coding RNAs (lncRNAs) under salt stress. BMC Plant Biol. 2018, 18, 23. [Google Scholar] [CrossRef] [PubMed]
  24. Fu, L.; Ding, Z.; Tan, D.; Han, B.; Sun, X.; Zhang, J. Genome-wide discovery and functional prediction of salt-responsive lncRNAs in duckweed. BMC Genom. 2020, 21, 212. [Google Scholar] [CrossRef]
  25. Zhang, X.; Dong, J.; Deng, F.; Wang, W.; Cheng, Y.; Song, L.; Hu, M.; Shen, J.; Xu, Q.; Shen, F. The long non-coding RNA lncRNA973 is involved in cotton response to salt stress. BMC Plant Biol. 2019, 19, 459. [Google Scholar] [CrossRef]
  26. Sun, Y.; Hao, P.; Lv, X.; Tian, J.; Wang, Y.; Zhang, X.; Xu, X.; Han, Z.; Wu, T. A long non-coding apple RNA, MSTRG.85814.11, acts as a transcriptional enhancer of SAUR32 and contributes to the Fe-deficiency response. Plant J. 2020, 103, 53–67. [Google Scholar] [CrossRef]
  27. Chen, G.; Cui, J.; Wang, L.; Zhu, Y.; Lu, Z.; Jin, B. Genome-Wide Identification of Circular RNAs in Arabidopsis thaliana. Front. Plant Sci. 2017, 8, 1678. [Google Scholar] [CrossRef]
  28. Hansen, T.B.; Jensen, T.I.; Clausen, B.H.; Bramsen, J.B.; Finsen, B.; Damgaard, C.K.; Kjems, J. Natural RNA circles function as efficient microRNA sponges. Nature 2013, 495, 384–388. [Google Scholar] [CrossRef]
  29. Conn, S.J.; Pillman, K.A.; Toubia, J.; Conn, V.M.; Salmanidis, M.; Phillips, C.A.; Roslan, S.; Schreiber, A.W.; Gregory, P.A.; Goodall, G.J. The RNA binding protein quaking regulates formation of circRNAs. Cell 2015, 160, 1125–1134. [Google Scholar] [CrossRef] [Green Version]
  30. Ye, C.; Chen, L.; Liu, C.; Zhu, Q.; Fan, L. Widespread noncoding circular RNAs in plants. New Phytol. 2015, 208, 88–95. [Google Scholar] [CrossRef]
  31. Wang, Y.; Yang, M.; Wei, S.; Qin, F.; Zhao, H.; Suo, B. Identification of Circular RNAs and Their Targets in Leaves of Triticum aestivum L. under Dehydration Stress. Front. Plant Sci. 2016, 7, 2024. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Zuo, J.; Wang, Q.; Zhu, B.; Luo, Y.; Gao, L. Deciphering the roles of circRNAs on chilling injury in tomato. Biochem. Biophys Res. Commun. 2016, 479, 132–138. [Google Scholar] [CrossRef] [PubMed]
  33. Li, C.; Qin, S.; Bao, L.; Guo, Z.; Zhao, L. Identification and functional prediction of circRNAs in Populus Euphratica Oliv. heteromorphic leaves. Genomics 2020, 112, 92–98. [Google Scholar] [CrossRef] [PubMed]
  34. Gao, Z.; Li, J.; Luo, M.; Li, H.; Chen, Q.; Wang, L.; Song, S.; Zhao, L.; Xu, W.; Zhang, C.; et al. Characterization and Cloning of Grape Circular RNAs Identified the Cold Resistance-Related Vv-circATS1. Plant Physiol. 2019, 180, 966–985. [Google Scholar] [CrossRef]
  35. Salmena, L.; Poliseno, L.; Tay, Y.; Kats, L.; Pandolfi, P.P. A ceRNA hypothesis: The Rosetta Stone of a hidden RNA language? Cell 2011, 146, 353–358. [Google Scholar] [CrossRef] [Green Version]
  36. Lu, Y.; Deng, S.; Li, Z.; Wu, J.; Liu, Q.; Liu, W.; Yu, W.; Zhang, Y.; Shi, W.; Zhou, J.; et al. Competing Endogenous RNA Networks Underlying Anatomical and Physiological Characteristics of Poplar Wood in Acclimation to Low Nitrogen Availability. Plant Cell Physiol. 2019, 60, 2478–2495. [Google Scholar] [CrossRef]
  37. Fu, X.; Zhang, X.; Qiu, J.; Zhou, X.; Yuan, M.; He, Y.; Chun, C.; Cao, L.; Ling, L.; Peng, L. Whole-transcriptome RNA sequencing reveals the global molecular responses and ceRNA regulatory network of mRNAs, lncRNAs, miRNAs and circRNAs in response to copper toxicity in Ziyang Xiangcheng (Citrus junos Sieb. Ex Tanaka). BMC Plant Biol. 2019, 19, 509. [Google Scholar] [CrossRef] [Green Version]
  38. He, X.; Guo, S.; Wang, Y.; Wang, L.; Shu, S.; Sun, J. Systematic identification and analysis of heat-stress-responsive lncRNAs, circRNAs and miRNAs with associated co-expression and ceRNA networks in cucumber (Cucumis sativus L.). Physiol. Plant 2020, 168, 736–754. [Google Scholar] [CrossRef]
  39. Skorupa, M.; Gołębiewski, M.; Kurnik, K.; Niedojadło, J.; Kęsy, J.; Klamkowski, K.; Wójcik, K.; Treder, W.; Tretyn, A.; Tyburski, J. Salt stress vs. salt shock-the case of sugar beet and its halophytic ancestor. BMC Plant Biol. 2019, 19, 57. [Google Scholar] [CrossRef] [Green Version]
  40. Shi, S.; Cui, J.; Lu, Z.; Cheng, D.; Luo, C. Screening of tolerance to NaCl in sugar beet germplasms. China Beet Sugar 2008, 4, 7–9. [Google Scholar]
  41. Dohm, J.C.; Minoche, A.E.; Holtgra¨we, D.; Capella-Gutie’rrez, S.; Zakrzewski, F.; Tafer, H.; Rupp, O.; RosleffSo¨rensen, T.; Stracke, R.; Reinhardt, R.; et al. The genome of the recently domesticated crop plant sugar beet (Beta vulgaris). Nature 2014, 505, 546–549. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Xie, Z.; Johansen, L.K.; Gustafson, A.M.; Kasschau, K.D.; Lellis, A.D.; Zilberman, D.; Jacobsen, S.E.; Carrington, J.C. Genetic and Functional Diversification of Small RNA Pathways in Plants. PLoS Biol. 2004, 2, e104. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Seto, T.; Kato, T.; Nishio, M.; Goto, K.; Atagi, S.; Hosomi, Y.; Yamamoto, N.; Hida, T.; Maemondo, M.; Nakagawa, K.; et al. Erlotinib alone or with bevacizumab as first-line therapy in patients with advanced non-squamous non-small-cell lung cancer harbouring EGFR mutations (JO25567): An open-label, randomised, multicentre, phase 2 study. Lancet Oncol. 2014, 15, 1236–1244. [Google Scholar] [CrossRef]
  44. Ala, U.; Karreth, F.A.; Bosia, C.; Pagnani, A.; Taulli, R.; Léopold, V.; Tay, Y.; Provero, P.; Zecchina, R.; Pandolfi, P.P. Integrated transcriptional and competitive endogenous RNA networks are cross-regulated in permissive molecular environments. Proc. Natl. Acad. Sci. USA 2013, 110, 7154–7159. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Joseph, S.; Murphy, D.; Bhave, M. Glycine betaine biosynthesis in saltbushes (Atriplex spp.) under salinity stress. Biologia 2013, 68, 879–895. [Google Scholar] [CrossRef]
  46. Zhang, L.; Wu, M.; Yu, D.; Teng, Y.; Wei, T.; Chen, C.; Song, W. Identification of Glutathione Peroxidase (GPX) Gene Family in Rhodiola crenulata and Gene Expression Analysis under Stress Conditions. Int. J. Mol. Sci. 2018, 19, 3329. [Google Scholar] [CrossRef] [Green Version]
  47. Bashir, K.; Rasheed, S.; Kobayashi, T.; Seki, M.; Nishizawa, N.K. Regulating Subcellular Metal Homeostasis: The Key to Crop Improvement. Front. Plant Sci. 2016, 7, 1192. [Google Scholar] [CrossRef] [Green Version]
  48. Li, J.; Cui, J.; Cheng, D.; Dai, C.; Liu, T.; Wang, C.; Luo, C. iTRAQ protein profile analysis of sugar beet under salt stress: Different coping mechanisms in leaves and roots. BMC Plant Biol. 2020, 20, 347. [Google Scholar] [CrossRef]
  49. Ma, C.; Burd, S.; Lers, A. miR408 is involved in abiotic stress responses in Arabidopsis. Plant J. 2015, 84, 169–187. [Google Scholar] [CrossRef]
  50. Hajyzadeh, M.; Turktas, M.; Khawar, K.M.; Unver, T. miR408 overexpression causes increased drought tolerance in chickpea. Gene 2015, 555, 186–193. [Google Scholar] [CrossRef]
  51. Chen, X.; Zhang, L.; Miao, X.; Hu, X.; Nan, S.; Wang, J.; Fu, H. Effect of salt stress on fatty acid and alpha-tocopherol metabolism in two desert shrub species. Planta 2018, 247, 499–511. [Google Scholar] [CrossRef] [PubMed]
  52. Shang, Y.; Yang, D.; Ha, Y.; Shin, H.-Y.; Nam, K.H. Receptor-like protein kinases RPK1 and BAK1 sequentially form complexes with the cytoplasmic kinase OST1 to regulate ABA-induced stomatal closure. J. Exp. Bot. 2020, 71, 1491–1502. [Google Scholar] [PubMed]
  53. Wei, M.; Zhuang, Y.; Li, H.; Li, P.; Huo, H.; Shu, D.; Huang, W.; Wang, S. The cloning and characterization of hypersensitive to salt stress mutant, affected in quinolinate synthase, highlights the involvement of NAD in stress-induced accumulation of ABA and proline. Plant J. 2020, 102, 85–98. [Google Scholar] [CrossRef]
  54. Ruan, Y.-L. Sucrose metabolism: Gateway to diverse carbon use and sugar signaling. Annu. Rev. Plant Biol. 2014, 65, 33–67. [Google Scholar] [CrossRef] [PubMed]
  55. Mclaughlin, J.E.; Boyer, J.S. Sugar-responsive gene expression, invertase activity, and senescence in aborting maize ovaries at low water potentials. Ann. Bot. 2004, 94, 675–689. [Google Scholar] [CrossRef]
  56. Roitsch, T.; Balibrea, M.E.; Hofmann, M.; Proels, R.; Sinha, A.K. Extracellular invertase: Key metabolic enzyme and PR protein. J. Exp. Bot. 2003, 54, 513–524. [Google Scholar] [CrossRef] [Green Version]
  57. Ma, H.; Liu, M. The microtubule cytoskeleton acts as a sensor for stress response signaling in plants. Mol. Biol. Rep. 2019, 46, 5603–5608. [Google Scholar] [CrossRef]
  58. Pua, E.-C.; Chandramouli, S.; Han, P.; Liu, P. Malate synthase gene expression during fruit ripening of Cavendish banana (Musa acuminata cv. Williams). J. Exp. Bot. 2003, 54, 309–316. [Google Scholar] [CrossRef]
  59. Mao, C.; Lu, S.; Lv, B.; Zhang, B.; Shen, J.; He, J.; Luo, L.; Xi, D.; Chen, X.; Ming, F. A Rice NAC Transcription Factor Promotes Leaf Senescence via ABA Biosynthesis. Plant Physiol. 2017, 174, 1747–1763. [Google Scholar] [CrossRef] [Green Version]
  60. Thirumalaikumar, V.P.; Devkar, V.; Mehterov, N.; Ali, S.; Ozgur, R.; Turkan, I.; Mueller-Roeber, B.; Balazadeh, S. NAC transcription factor JUNGBRUNNEN1 enhances drought tolerance in tomato. Plant Biotechnol. J. 2018, 16, 354–366. [Google Scholar] [CrossRef]
  61. Liu, J.; Fan, H.; Wang, Y.; Han, C.; Wang, X.; Yu, J.; Li, D.; Zhang, Y. Genome-Wide microRNA Profiling Using Oligonucleotide Microarray Reveals Regulatory Networks of microRNAs in Nicotiana benthamiana During Beet Necrotic Yellow Vein Virus Infection. Viruses 2020, 12, 310. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Hichri, I.; Muhovski, Y.; Clippe, A.; Žižková, E.; Dobrev, P.I.; Motyka, V.; Lutts, S. SlDREB2, a tomato dehydration-responsive element-binding 2 transcription factor, mediates salt stress tolerance in tomato and Arabidopsis. Plant Cell Env. 2016, 39, 62–79. [Google Scholar] [CrossRef] [PubMed]
  63. Bouaziz, D.; Pirrello, J.; Amor, H.B.; Hammami, A.; Charfeddine, M.; Dhieb, A.; Bouzayen, M.; Gargouri-Bouzid, R. Ectopic expression of dehydration responsive element binding proteins (StDREB2) confers higher tolerance to salt stress in potato. Plant Physiol. Biochem. 2012, 60, 98–108. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Gao, J. Plant Physiology Experiment Guide; Higher Education Press: Beijing, China, 2006. [Google Scholar]
  65. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet 2011, 17, 10–12. [Google Scholar] [CrossRef]
  66. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef] [Green Version]
  67. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef] [Green Version]
  68. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.-C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290–295. [Google Scholar] [CrossRef] [Green Version]
  69. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2009, 26, 139–140. [Google Scholar] [CrossRef] [Green Version]
  70. Conway, J.R.; Lex, A.; Gehlenborg, N. UpSetR: An R package for the visualization of intersecting sets and their properties. Bioinformatics 2017, 33, 2938–2940. [Google Scholar] [CrossRef] [Green Version]
  71. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T.; et al. Gene ontology: Tool for the unification of biology. Nat. Genet. 2000, 25, 25–29. [Google Scholar] [CrossRef] [Green Version]
  72. Kanehisa, M.; Araki, M.; Goto, S.; Hattori, M.; Hirakawa, M.; Itoh, M.; Katayama, T.; Kawashima, S.; Okuda, S.; Tokimatsu, T.; et al. KEGG for linking genomes to life and the environment. Nucleic. Acids Res. 2007, 36, D480–D484. [Google Scholar] [CrossRef] [PubMed]
  73. Conesa, A.; Götz, S.; García-Gómez, J.M.; Terol, J.; Talón, M.; Robles, M. Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 2005, 21, 3674–3676. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  74. Kim, D.; Salzberg, S.L. TopHat-Fusion: An algorithm for discovery of novel fusion transcripts. Genome. Biol. 2011, 12, R72. [Google Scholar] [CrossRef] [Green Version]
  75. Zhang, X.; Wang, H.; Zhang, Y.; Lu, X.; Chen, L.L.; Yang, L. Complementary sequence-mediated exon circularization. Cell 2014, 159, 134–147. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  76. Zhang, X.; Dong, R.; Zhang, Y.; Zhang, J.; Luo, Z.; Zhang, J.; Chen, L.-L.; Yang, L. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome. Res. 2016, 26, 1277–1287. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  77. Gao, Y.; Wang, J.; Zhao, F. CIRI: An efficient and unbiased algorithm for de novo circular RNA identification. Genome. Biol. 2015, 16, 4. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  78. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome. Biol. 2009, 10, R25. [Google Scholar] [CrossRef] [Green Version]
  79. Zuker, M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic. Acids Res. 2003, 31, 3406–3415. [Google Scholar] [CrossRef]
  80. Li, X.; Shahid, M.Q.; Wu, J.; Wang, L.; Liu, X.; Lu, Y. Comparative Small RNA Analysis of Pollen Development in Autotetraploid and Diploid Rice. Int. J. Mol. Sci. 2016, 17, 499. [Google Scholar] [CrossRef]
  81. Ma, Z.; Coruh, C.; Axtell, M.J. Arabidopsis lyrata small RNAs: Transient MIRNA and small interfering RNA loci within the Arabidopsis genus. Plant Cell 2010, 22, 1090–1103. [Google Scholar] [CrossRef] [Green Version]
  82. Addo-Quaye, C.; Miller, W.; Axtell, M.J. CleaveLand: A pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics 2009, 25, 130–131. [Google Scholar] [CrossRef] [PubMed]
  83. Xu, X.; Yin, L.; Ying, Q.; Song, H.; Xue, D.; Lai, T.; Xu, M.; Shen, B.; Wang, H.; Shi, X. High-Throughput Sequencing and Degradome Analysis Identify miRNAs and Their Targets Involved in Fruit Senescence of Fragaria ananassa. PLoS ONE 2013, 8, e70959. [Google Scholar] [CrossRef] [PubMed]
  84. Meng, Y.; Shao, C.; Wang, H.; Jin, Y. Target mimics_ an embedded layer of microRNA-involved gene regulatory networks in plants. Bmc Genom. 2012, 13, 197. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  85. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideke, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res 2003, 13, 2498–2504. [Google Scholar] [CrossRef] [PubMed]
  86. Kramer, M.F. Stem-Loop RT-qPCR for miRNAs. Curr. Protoc. Mol. Biol. 2011, 95, 15.10.1–15.10.15. [Google Scholar] [CrossRef] [PubMed]
  87. Pfaffl, M.W.; Tichopad, A.; Prgomet, C.; Neuvians, T.P. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper-Excel-based tool using pair-wise correlations. Biotechnol. Lett. 2004, 26, 509–515. [Google Scholar] [CrossRef]
  88. Andersen, C.L.; Jensen, J.L.; Ørntoft, T.F. Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004, 64, 5245–5250. [Google Scholar] [CrossRef] [Green Version]
  89. Vandesompele, J.; Preter, K.D.; Pattyn, F.; Poppe, B.; Roy, N.V.; Paepe, A.D.; Speleman, F. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome. Biol. 2002, 3, research0034.1. [Google Scholar] [CrossRef] [Green Version]
Ijms 22 00289 g001 550
Figure 1. Identification and analysis of differentially expressed mRNAs (DEmRNAs) under salt stress. (a) Violin Plot of gene expression patterns for each sample, with origin representing the median; ck represents the control group and ST represents the treatment group. (b) Volcano Plot of log2 FC(ST/CK) of leaves and roots mRNAs. (c) Visualization of DEmRNAs in each sample, where the intersection size represents those common to all samples, intersection size stands for the common DEmRNAs among different samples. (d–e) gene ontology (GO) annotation in leaves (d) and roots (e). (f,g) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment in leaves (f) and roots (g). (h,i) Mapman analysis of metabolism in leaves (h) and in roots (i), where each square represents mapping to the metabolic pathway and color represents up-regulation (red) or down-regulation (blue).
Figure 1. Identification and analysis of differentially expressed mRNAs (DEmRNAs) under salt stress. (a) Violin Plot of gene expression patterns for each sample, with origin representing the median; ck represents the control group and ST represents the treatment group. (b) Volcano Plot of log2 FC(ST/CK) of leaves and roots mRNAs. (c) Visualization of DEmRNAs in each sample, where the intersection size represents those common to all samples, intersection size stands for the common DEmRNAs among different samples. (d–e) gene ontology (GO) annotation in leaves (d) and roots (e). (f,g) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment in leaves (f) and roots (g). (h,i) Mapman analysis of metabolism in leaves (h) and in roots (i), where each square represents mapping to the metabolic pathway and color represents up-regulation (red) or down-regulation (blue).
Ijms 22 00289 g001
Ijms 22 00289 g002 550
Figure 2. Identification and analysis of DElncRNAs under salt stress. (a–d) Comparison of transcript length, exon number, fragments per kilobase of transcript per million mapped reads (FPKM) value, numbers and ORF length in lncRNA versus mRNA. (e) Venn diagram of DElncRNAs in leaves and roots. (f,g) Heat map of DElncRNAs in leaves (g) and in roots. (h) GO enrichment of targets of DElncRNAs in roots; (i) KEGG enrichment of targets of roots DElncRNAs.
Figure 2. Identification and analysis of DElncRNAs under salt stress. (a–d) Comparison of transcript length, exon number, fragments per kilobase of transcript per million mapped reads (FPKM) value, numbers and ORF length in lncRNA versus mRNA. (e) Venn diagram of DElncRNAs in leaves and roots. (f,g) Heat map of DElncRNAs in leaves (g) and in roots. (h) GO enrichment of targets of DElncRNAs in roots; (i) KEGG enrichment of targets of roots DElncRNAs.
Ijms 22 00289 g002
Ijms 22 00289 g003 550
Figure 3. Identification and analysis of DEcircRNAs under salt stress. (a,b) All identified circRNAs chromosome distribution (a) and type distribution (b). (c) Venn diagram of DEcircRNAs in leaves and roots. (d) Number of DEcircRNAs in leaf and root. (e,f) Heat map of DEcircRNAs in leaves (e) and in roots (f).
Figure 3. Identification and analysis of DEcircRNAs under salt stress. (a,b) All identified circRNAs chromosome distribution (a) and type distribution (b). (c) Venn diagram of DEcircRNAs in leaves and roots. (d) Number of DEcircRNAs in leaf and root. (e,f) Heat map of DEcircRNAs in leaves (e) and in roots (f).
Ijms 22 00289 g003
Ijms 22 00289 g004 550
Figure 4. Identification and analysis of DEmiRNAs under salt stress. (a,b) Identified miRNAs length distribution (a) and 5′ nucleotide composition (b). (c) Venn diagram of identified miRNAs in each sample. (d,e) Heat map of DEmiRNAs in leaves (d) and in roots (e). (f,g) GO annotation of targets of DEmiRNAs in leaves (f) and in roots (g).
Figure 4. Identification and analysis of DEmiRNAs under salt stress. (a,b) Identified miRNAs length distribution (a) and 5′ nucleotide composition (b). (c) Venn diagram of identified miRNAs in each sample. (d,e) Heat map of DEmiRNAs in leaves (d) and in roots (e). (f,g) GO annotation of targets of DEmiRNAs in leaves (f) and in roots (g).
Ijms 22 00289 g004
Ijms 22 00289 g005 550
Figure 5. Salt stress response ceRNA network constructed with DEmRNAs, DElncRNAs and DEmiRNAs in leaves (a) and roots (b). The color represents up-regulation (red) or down-regulation (green). Red solid lines represent regulation verified by degradome sequencing.
Figure 5. Salt stress response ceRNA network constructed with DEmRNAs, DElncRNAs and DEmiRNAs in leaves (a) and roots (b). The color represents up-regulation (red) or down-regulation (green). Red solid lines represent regulation verified by degradome sequencing.
Ijms 22 00289 g005
Ijms 22 00289 g006 550
Figure 6. T-plots and alignments of miRNA–mRNA pairs validated by degradome sequencing. (a) bv1_023200_jmkt.t1 cleaved by gma-miR408a-3p_L-1R+5; (b) bv6_144730_qgsa.t1 cleaved by mtr-miR408-3p_L-1R+1; (c) bv1_023200_jmkt.t1 cleaved by mtr-miR408-3p_L-1R+1; (d) bv1_023200_jmkt.t1 cleaved by mtr-MIR408-p3_2ss18GT19GA; (e) bv6_144730_qgsa.t1 cleaved by mtr-MIR408-p3_2ss18GT19GA; (f) bv5_114390_pjnp.t1 cleaved by mtr-miR164d; (g) bv6_136680_juie.t1 cleaved by PC-3p-154_19269; (h) bv6_136690_fscn.t1 cleaved by PC-3p-154_19269; (i) bv_009820_mrec.t1 cleaved by mtr-miR395k_1ss1TC. The ‘lightning’ shape indicates the position of the cleavage site.
Figure 6. T-plots and alignments of miRNA–mRNA pairs validated by degradome sequencing. (a) bv1_023200_jmkt.t1 cleaved by gma-miR408a-3p_L-1R+5; (b) bv6_144730_qgsa.t1 cleaved by mtr-miR408-3p_L-1R+1; (c) bv1_023200_jmkt.t1 cleaved by mtr-miR408-3p_L-1R+1; (d) bv1_023200_jmkt.t1 cleaved by mtr-MIR408-p3_2ss18GT19GA; (e) bv6_144730_qgsa.t1 cleaved by mtr-MIR408-p3_2ss18GT19GA; (f) bv5_114390_pjnp.t1 cleaved by mtr-miR164d; (g) bv6_136680_juie.t1 cleaved by PC-3p-154_19269; (h) bv6_136690_fscn.t1 cleaved by PC-3p-154_19269; (i) bv_009820_mrec.t1 cleaved by mtr-miR395k_1ss1TC. The ‘lightning’ shape indicates the position of the cleavage site.
Ijms 22 00289 g006
Ijms 22 00289 g007 550
Figure 7. qRT-PCR analysis of all 50 components of the ceRNA network under salt stress in (a) leaves and (b) roots. FC(ST/CK) represents fold changes of relative expression levels between ST and CK, log2FC(ST/CK) is from the mean of three replicates, bars stand for ±SD, asterisk (*) indicates that the gene is not detected in ST and its value is artificially set to −4 instead of calculating log2(0).
Figure 7. qRT-PCR analysis of all 50 components of the ceRNA network under salt stress in (a) leaves and (b) roots. FC(ST/CK) represents fold changes of relative expression levels between ST and CK, log2FC(ST/CK) is from the mean of three replicates, bars stand for ±SD, asterisk (*) indicates that the gene is not detected in ST and its value is artificially set to −4 instead of calculating log2(0).
Ijms 22 00289 g007
Ijms 22 00289 g008 550
Figure 8. Proposed model of ceRNA network that regulates response to salt stress in sugar beet. SS: sucrose synthase, PFK: fructose phosphate kinase, HK: hexokinase, PK: pyruvate kinase, PDH: pyruvate dehydrogenase, MDH: malic dehydrogenase. Red filling represents up-regulation, whereas green filling represents down-regulation.
Figure 8. Proposed model of ceRNA network that regulates response to salt stress in sugar beet. SS: sucrose synthase, PFK: fructose phosphate kinase, HK: hexokinase, PK: pyruvate kinase, PDH: pyruvate dehydrogenase, MDH: malic dehydrogenase. Red filling represents up-regulation, whereas green filling represents down-regulation.
Ijms 22 00289 g008
Table
Table 1. Physiological parameters under salt stress.
Table 1. Physiological parameters under salt stress.
Physiological Indicesck (Leaf)st (Leaf)ck (Root)st (Root)
Relative water content (%)90.17 ± 0.06 *87.73 ± 0.15 *--
Chlorophyll (mg·g−1 FW)2.104 ± 0.069 *1.759 ± 0.096 *--
Tocopherol (μg·g−1 FW)223.22 ± 14.70 *164.50 ± 10.87 *--
Soluble sugar (mg·g−1 FW)0.3664 ± 0.0186 *0.5835 ± 0.0121 *0.2792 ± 0.0067 *0.3309 ± 0.0151 *
MDA (nmol·g−1 FW)8.3427 ± 0.14649.4161 ± 0.84506.6839 ± 0.08458.1476 ± 0.5915
Proline (μg·g−1 FW)1.4131 ± 0.1915 *4.0466 ± 0.0335 *0.5334 ± 0.02130.5897 ± 0.0563
POD activity (U·g−1 FW)48.14 ± 1.2644.43 ± 1.29121.85 ± 13.85 *206.85 ± 8.66 *
SOD activity (U·g−1 FW)24.57 ± 0.54 *41.07 ± 6.28 *24.62 ± 1.33 *64.09 ± 2.49 *
CAT activity (U·g−1 FW)568.78 ± 2.71569.36 ± 2.8936.45 ± 1.28 *21.07 ± 1.15 *
Data are means ± SD (n = 3). * indicate p < 0.05 by Tukey test. FW: fresh weight; MDA: malonaldehyde; POD: peroxidase; SOD: superoxide dismutase; CAT: catalase.
Table
Table 2. Negative regulation DEmiRNA–mRNA pairs in leaves and roots under salt stress.
Table 2. Negative regulation DEmiRNA–mRNA pairs in leaves and roots under salt stress.
miR_nameup/downTranscriptup/downTissueDegradome Detection
mtr-miR408-3p_L-1R+1upBv6_144730_qgsa.t1downleafY
mtr-miR408-3p_L-1R+1upBv1_023200_jmkt.t1downleafY
gma-miR408a-3p_L-1R+5upBv1_023200_jmkt.t1downleafY
PC-3p-154_19269upBv6_136680_juie.t1downleafY
PC-3p-154_19269upBv6_136690_fscn.t1downleafY
mtr-miR395a_L-1downBv_009820_mrec.t1upleafY
mtr-miR395k_1ss1TCdownBv_009820_mrec.t1upleafY
gma-miR403a_1ss20TCdownBv4_094840_csri.t1upleafY
mtr-miR395a_L-1downBv_010680_tsww.t1upleafN
lus-MIR396e-p5downBv1_015660_nghm.t1upleafN
PC-3p-726_3835downMSTRG.2715.1upleafN
PC-5p-109182_51downBv4_078360_ftyu.t1upleafN
PC-5p-36464_148downBv1_011830_oexd.t1upleafN
PC-5p-3682_979downBv5_114160_hsdp.t1upleafN
PC-5p-3682_979downBv8_193090_kary.t1upleafN
PC-5p-49652_112downMSTRG.29155.1upleafN
csi-miR156a-5p_R+1_1ss9GTupBv6_136190_cygi.t1downrootY
mtr-miR164dupBv5_114390_pjnp.t1downrootY
PC-5p-7955_523upBv5_124210_skcr.t1downrootY
gma-miR6300_L-1R+1upBv5_097930_juac.t1downrootN
gma-miR6300_L-1R+1upBv5_100510_ttjr.t1downrootN
gma-miR6300_L-1R+1upBv9_208900_gijo.t1downrootN
PC-5p-3682_979upBv1_012390_utfq.t1downrootN
PC-5p-3682_979upBv1_017030_kxfa.t1downrootN
PC-5p-3682_979upBv5_117790_mxei.t1downrootN
PC-5p-3682_979upBv6_132150_dsqx.t1downrootN
PC-5p-3682_979upBv6_135930_aphq.t1downrootN
PC-5p-3682_979upBv8_193090_kary.t1downrootN
ptc-miR399edownBv5_121080_kswa.t1uprootN
PC-5p-109182_51downBv5_110660_jxri.t1uprootN
PC-5p-109182_51downBv9_213960_qddk.t1uprootN
PC-5p-130970_42downMSTRG.29564.3uprootN
PC-5p-36464_148downMSTRG.25036.1uprootN
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, J.; Cui, J.; Dai, C.; Liu, T.; Cheng, D.; Luo, C. Whole-Transcriptome RNA Sequencing Reveals the Global Molecular Responses and CeRNA Regulatory Network of mRNAs, lncRNAs, miRNAs and circRNAs in Response to Salt Stress in Sugar Beet (Beta vulgaris). Int. J. Mol. Sci. 2021, 22, 289. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22010289

AMA Style

Li J, Cui J, Dai C, Liu T, Cheng D, Luo C. Whole-Transcriptome RNA Sequencing Reveals the Global Molecular Responses and CeRNA Regulatory Network of mRNAs, lncRNAs, miRNAs and circRNAs in Response to Salt Stress in Sugar Beet (Beta vulgaris). International Journal of Molecular Sciences. 2021; 22(1):289. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22010289

Chicago/Turabian Style

Li, Junliang, Jie Cui, Cuihong Dai, Tianjiao Liu, Dayou Cheng, and Chengfei Luo. 2021. "Whole-Transcriptome RNA Sequencing Reveals the Global Molecular Responses and CeRNA Regulatory Network of mRNAs, lncRNAs, miRNAs and circRNAs in Response to Salt Stress in Sugar Beet (Beta vulgaris)" International Journal of Molecular Sciences 22, no. 1: 289. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22010289

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop