Combined RNA-seq, small RNA and degradome sequencing approaches insights into salt-stress responses in Beta vulgaris

Background Salinity is one of the most serious threat to agriculture worldwide. Sugar beet is an important sugar-yielding crop and has a certain tolerance to salt. While a genome-wide analysis of genes involved in salt-stress remains largely unknown in B. vulgaris. Results Three high-throughput sequencing approaches, namely, RNA-seq, small RNA and degradome sequencing were used to exploring the molecular basis of salt-resistance in beta vulgaris. A total of 12 230 differentially expressed genes were identified, which were mainly related to transcription factors, protein kinases, and enzymes. The small RNA sequencing resulted in the identification of 476 miRNAs included 219 known and 257 novel miRNAs, of which 130 were differentially expressed under salt-stress. A total of 2534 targets were predicted for 317 miRNAs by bioinformatics and degradome sequencing. Functional analysis of these targets suggested that miRNA-mediated post-transcriptional regulation plays a crucial role in transcriptional reprogramming under stress. The combination analysis of multi-omics revealed 209 negative correlation mRNA-miRNA pairs. All these data were used to construct a hypothetical regulatory network of beta vulgaris in response to salt stress. These components and the preliminary network provides better insights into the molecular mechanism of salt-stress response, and also offers candidate genes for beet improvement.


Introduction
Salinity is one of the most severe abiotic threats that affects the growth and development of crops (Pichu Rengasamy 2010 andRana Munns 2015). High concentrations of variety salts, such as chlorides, sulfates and carbonates of sodium, magnesium and calcium, make up most of the saline areas all over the world. Soil-salinization is a growing problem for agriculture that may negatively decrease the quality and yield of crops. The common effect of soil salinity on plants comes from the inhibition of growth by Na + and Clˉ accumulation (Tester 2003). Although the growth of most crops is adversely affected by soil salinity, some cultivars are able to adapt to saline conditions to achieve good harvests. In this context, breeding and propagation of salt-tolerant crops is a promising solution of agronomic and engineering to address these challenges.
Sugar beet (Beta vulgaris ssp. Vulgaris or B. vulgaris) is one of the most important sugar-yielding crops in the world. As recently domesticated crop, cultivated beets inherited certain salt-tolerance traits from its wild ancestor Beta vulgaris ssp. maritima (referred to as B. maritima or 'sea beet') (Monika 2019). In the past few decades, much research has been devoted to the excavation of salttolerant genes in B. vulgaris (Antonio 2003, Wu 2013, Kamila 2014and Kito 2017. In this study, sugar beet 'O68' was selected for studying the molecular response of sugar beet to salt-stress. 'O68' was a high yield and high resistance cultivar breeding from a cross between "Shuang feng 3" (a kind of high yield variety bred and selected from West German K.Lelnwanzlchen-E) and "A Yu 1" (a kind of resistance to Cercospora Leaf Spot and high sugar content varieties bred and selected from Poland K.BusZCYZYnski-CLR). Our previous study showed that compared with the control, the relative germination rate of this cultivar was more than 90% and the seedling can grow well under 200 m mol·L -1 NaCl treatment, while the germination rate and growth state were not significantly affected until the treatment concentration reached 300 m mol·L -1 (Shi 2008). In addition, it has a strong regeneration capability of petiole explants and readily infected by agrobacterium which is ideally suited for use in molecular breeding. At the same time, it is one of the most excellent parent used in traditional crossbreeding. Therefore, O68 is a good choice for studying the mechanism of salt-stress response in B. vulgaris.
Unlike other abiotic stresses, salinity causes both osmotic stress and ion toxicity in plants (Munns 2008). Plant growth may rapidly impaired by osmotic stress in a first phase and then specific ion toxicity primarily from Na + and Clˉ accumulation may cause membrane disorganization, the generation of reactive oxygen species, metabolic toxicity, inhibition of photosynthesis, and the attenuation of nutrient acquisition in a second phase of salt stress (MUNNS 1993 andGarthwaite 2005). In plant, excess ingestion of Na + can affect the absorption of mineral nutrients, such as calcium (Ca 2+ ), magnesium (Mg 2+ ) and potassium (K + ) (Antonio 2003 andParida 2005 andMunns 4 2008). However, absorbs large amounts of Na+ reduced the absorption of K + , the growth and development of B. vulgaris were not significantly affected. This is mainly due to the use of Na + replaces the function of K + , like achieve long-distance transport of anions, stomatal regulation and Osmotic regulation (Subbarao 1999, Subbarao 2000and Pi 2016. Compared to the achievement of physiological level, a genome-wide analysis of genes involved in salt-stress remains largely unknown in B. vulgaris. RNA-seq is an excellent way for the transcript identification and expression quantification which has been widely used to study gene expression changes at the transcriptional level. MicroRNAs (miRNAs) are small, endogenous, noncoding RNAs ranging from 19 to 24 nucleotides that play crucial regulatory roles at the post-transcriptional stage by cleavage or translation repression of their target mRNAs (Bartel 2004). Increasing evidence suggests that miRNA-guided gene regulation is involved in multiple physiological and biochemical processes, such as seedling growth (Hu 2016), floral development (Chen 2004), male and female sterility (Millar 2005 andYan 2015), and responses to environmental stress (Sunkar 2007). Furthermore, a new study has found that plants can export miRNAs to inhibit virulence gene expression in a fungal pathogen (Zhang 2016).Thus, identifying miRNAs and their target mRNAs is essential to elucidating their regulatory mechanisms, and to understand the growth, development, and stress response processes. In plants, most miRNAs have perfect or near-perfect sequences that complement their target mRNAs to function and can be predicted using computational tools (Jones-Rhoades 2006). Moreover, a modified version of 5'-Rapid Amplification of cDNA Ends (RACE) using high-throughput deep sequencing technology named degradome sequencing has been increasingly used to predict targets of miRNAs (Addo-Quaye 2008 and German 2008). Many of the previously validated and predicted targets of miRNAs have been verified using this method, which provides a comprehensive means of analyzing patterns of RNA degradation. This combined with small RNA sequencing has greatly increased the rate of identifying miRNAs and their targets in plants (Gutierrez 2012).
With an objective to get deeper insights into the salt resistance mechanisms in B. vulgaris, the 5 present study employed RNA-seq, small RNA sequencing and degradome sequencing using cultivar 'O68' treated with 300 m mol·L -1 NaCl for different time. A comprehensive and integrated analysis of these different datasets have identified 12 231 differentially expressed genes (DEGs), 130 DE miRNAs and 41 positive correlation DE miRNA-mRNA pairs in B. vulgaris in response to salt-stress. This data will provide some insights into the miRNA-mediated regulatory pathways in salt-stress responses and improve the understanding of the salt-stress response of B. vulgaris.

Transcriptome sequencing
Ten transcriptome libraries was constructed from mixed RNA pools consisting of leaf / root samples from six-leaf-stage seedlings treated with 300 mM NaCl 2 for 5 durations (0 hour, 12 hour, 24 hour, 48 hour and 72 hour). Then, more than 0.2 billion reads were generated from the SE 50 sequencing, with approximately 21.9 million reads for each sample. After quality control, 97.7% (216.3 million) of the reads representing high-quality reads were processed for further analysis. These reads were mapped on the beet genome with over 80% uniquely mapping ratio. The details information of RNA-seq data generated, filtered reads and mapped ratios are given in Table 1. After quantifying gene expression levels as fragments per kilobase of transcript per million mapped reads (FPKM), a total of 24 729 genes were identified as expressed genes in B. vulgaris under salt stress (Supplementary file 1: Table   S1). The number and proportion of genes at each stage indicated the gene expression levels were generally comparable among different samples ( Figure 1a). In each tissue, the gene expression levels at any two given stages were highly correlated, with the average correlation coefficient (r) above 0.8.
Interestingly, the sample dendrogram showed that genes at root_12 and leaf_24 had a relatively far distance compared with other time points (Supplementary file 2: Figure S1a).

Differential gene expression analysis
To investigate the expression patterns of genes in B. vulgaris under salt stress, the genes with very low expression values (FPKM< 1) in all the samples were filtered out. Under this standard, a total of 21 849 genes were found to be expressed in at least one of the samples. After filtering the lowly expressed genes, we identified 12 230 genes as differentially expressed genes (DEGs) using STEM, 6 including7 800 and 7 764 in root and leaf, respectively (Supplementary file 1: Table S2). Under salt stress, the highest number of DEGs (4 874) were identified from root between 12 h and 0 h, and the trend of down-regulation was more intense at 12h (3510) compared to control in root (Figure 1b). A similar situation in leaf was occurred at 24 h (3134) in contrast to control. Further, an overlap between different time points was analyzed in root and leaf, respectively. It was observed that only a small fraction of DEGs was unique for each time point. In root, 12h exhibited highest number of stagespecific stress-responsive DEGs (1921) in contrast to 0h (Supplementary file 2: Figure S1b). While in leaf, 24h exhibited highest number of stage-specific stress-responsive DEGs (838) in contrast to 0h (Supplementary file 2: Figure S1c). The genomic localization and an overall expression of the DEGs at five time points in root and leaf are illustrated in Figure 1 c and d, respectively. Overall, these results indicate tissue-and stage-specific response to salt stress in B. vulgaris.
Out of 12 231 significant DEGs, 11 524 were annotated using blastx. It should be noted that a total of 555 transcription factor (TF) encoding genes belonging to 53 families were identified from DEGs Phenylpropanoid biosynthesis (ko00940) were among the most enriched pathways (Figure 2d).

Expression trends analysis during five time point
A large number of transcriptional processes were involved in salt-response in plant. Therefore, cluster analysis was used to sort out genes with similar expression trends at five time point under salt-stress.
Hierarchical clustering of the DEGs revealed three groups of genes based on their expression trends in root and leaf, respectively. In root, Group I genes exhibited a gradually decreased expression. In contrast, Group III genes had an increased expression. Unlike them, the expression of genes in Group II were predominantly decreased at 12 h stage and increase subsequently ( Figure 3a). While in leaf, Group I genes had gradually increased expression, Group II genes exhibited decreased expression and Group III genes exhibited decrease first and then increase expression (Figure 3c). To explore the possible functions of the DEGs in response to salt-stress, k-means clustering analysis was used to further divide the DEGs into six clusters with clear and distinct expression profiles, simultaneously, GO enrichment analysis and pathway analysis were individually performed for each gene cluster.
In root (Figure 3b and Supplementary file 2: Figure S3 a), cluster 1 contained genes showed repressed expression in contrast to control under salt-stress. The GO enrichment analysis revealed that most of these genes were involved in translation, defense response, plant-type cell wall organization and 8 phosphorelay signal transduction system. A number of ribosome protein genes like 40S ribosomal protein S5 (Bv3_051190_tnuc), 60S ribosomal protein L15 (Bv8_188220_yfhe) along with defense response genes like defensin-like protein (Bv4_083710_eume) and Expansin protein genes like expansin-A12 (Bv_30600_egdm) were significant members of this cluster. The pathway enrichment suggested that these genes were mostly participated in ribosome and metabolism of nucleotide. The expression pattern of genes in cluster 2 revealed significantly up-regulation from 0 h to 24 h and then decreased gradually over time. This cluster mainly contained genes involved in fatty acid biosynthetic process, glucuronoxylan biosynthetic process, one-carbon metabolic process and plant-type secondary cell wall biogenesis. Cellulose synthase A (Bv8_195690_hdar) and 3-ketoacyl-CoA synthase 1 (Bv1_002370_dexr) were included in this cluster. Pathway analysis showed that genes in this cluster were took part in biosynthesis of secondary metabolites and metabolic pathways. Cluster 3 contained genes highly-expressed at 12 h and then down-regulated until stable. These genes were mainly associated with response to heat and water deprivation, glycolytic process and abscisic acid-activated signaling pathway like heat shock 70 kDa protein gene (Bv_29970_uwua), hexokinase-1 (Bv9u_231360_aurc), EID1-like F-box protein (Bv_35170_scdu) and mitogen-activated protein kinase homolog MMK2 (Bv5_094830_onxe). The pathway enrichment results implicated that these genes play an important role in pathogen interaction and MAPK signaling pathway. Genes in cluster 4 were significant down-regulated at 12 h and return to normal level after 24 h. It contained genes involved in microtubule-based movement like kinesin protein NACK1 (Bv_02120_kzid), DNA replication initiation like DNA replication licensing factor MCM2 (Bv8_191930_unzu) and cellulose catabolic process like endoglucanase 17(Bv7_164140_fgmh). Interestingly, a large number of TFs genes like ABI3/VP1, ARF, bHLH, ERF and MYB were also included in this cluster. Pathway enrichment analysis indicated that genes in this cluster were mostly participated in DNA replication and amino sugar and nucleotide sugar metabolism pathways. Cluster 5 and cluster 6 contained genes with gradually increased expression under salt-stress. Further, genes in cluster 5 were induced to express under salt-stress. The genes of this cluster were mainly involved in protein polyubiquitination, malate metabolic process and ATP synthesis coupled electron transport.E3 ubiquitin-protein ligase 9 (Bv6_153190_wfqg), malate dehydrogenase (Bv_33060_gkqg) and TFs encoding for ERF, MYB, ZAT5 and WRKY were included in this cluster. The pathway analysis suggested these genes were primarily included in RNA polymerase and ubiquinone biosynthesis. While cluster 6 contained genes highlyexpressed from 24 h to 72 h. These genes were mainly associated with brassinosteroid biosynthetic process and cytochrome b6f complex assembly like dehydrocholesterol reductase (Bv4u_092080_fdde), thioredoxin (Bv2_034380_eidr). And genes in this cluster were mostly participated in peroxisome and carotenoid biosynthesis.
As well as root, the six clusters in leaf also showed significant differences in gene function ( Figure 3d and Supplementary file 2: Figure S3 b). Cluster 1 contained genes induced to express under saltstress. GO enrichment analysis revealed that these genes were mainly associated with cell wall organization, phosphorelay signal transduction system, plant-type secondary cell wall biogenesis, amine metabolic process, glucuronoxylan biosynthetic process and iron-sulfur cluster assembly.
Genes like cellulose synthase A catalytic subunit 4 (Bv_41430_krdf), two-component response regulator APRR5 (Bv3_051580_qrsy), histidine-containing phosphotransfer protein 4 (Bv9u_231480_afwz), primary amine oxidase (Bv5_093200_xpuu) and SufD (Bv3u_070060_mirg) were included in this cluster. The pathway enrichment results implicated that genes in this cluster were mostly participated in biosynthesis of phenylpropanoid and secondary metabolites. Genes in cluster 2 showed a significantly inhibited expression from 12 h to 48 h and then recovered slightly at 72 h in contrast to control. Most of these genes were involved in metabolic process, trehalose biosynthetic process and ATP synthesis coupled proton transport like trehalose-phosphate phosphatase A (Bv3u_070160_npqs), ATP synthase CF1 (Bv_70050_rwta) along with TFs like ARF, bHLH, ERF, MYB.
The pathway enrichment suggested that these genes were mostly participated in signal transduction and metabolism of starch and sucrose. While the expression of genes in cluster 3 presented repressed compared to control during salt-stress. Most of these genes were associated with translation and nucleosome assembly. A number of ribosome protein genes like 60S ribosomal protein L21 (Bv_02040_uxcz) and 40S ribosomal protein S13 (Bv7_161260_xjzx) were enrich in this cluster.
Consistent with this, genes in this cluster were mostly participated in ribosome and biosynthesis of amino acids pathways. Cluster 4 contained genes significant down-regulated at 12 h and gradually return to normal level after 48 h. It contained genes involved in RNA modification like PPR protein (Bv_51200_cehg) and rRNA processing like rRNA biogenesis protein RRP5 (Bv7u_183220_ogrd).
Pathway enrichment analysis indicated that genes in this cluster were mostly participated in homologous recombination, ribosome biogenesis and nucleotide metabolism. Cluster 5 contained genes had the wave-like increasing tendency under salt-stress. The genes of this cluster were mainly involved in lipid metabolic process and regulation of transcription, DNA-templated. Scarecrow-like protein 15 (Bv_47070_aoao), phosphatidylinositol phosphodiesterase (Bv9_220370_gyof) and TFs encoding for bHLH, bZIP, ERF, EIL, HSF and MYB were included in this cluster. Pathway analysis showed that genes in this cluster were mainly took part in peroxisome and plant hormone signal transduction pathway. The expression pattern of genes in cluster 6 revealed significantly downregulation from 0 h to 24 h and then gradually return to primal level at 72 h. Genes in this cluster were primarily associated with DNA replication, photosystem II repair, cell wall organization or biogenesis, cell surface receptor signaling pathway and response to salicylic acid. G-protein coupled receptor 1 (Bv2u_047890_ddzh), origin recognition complex subunit 3 (Bv_64640_kphq), HtrA serine peptidase 2 (Bv7_168040_mpxt) and UDP-glucose: glycoprotein glucosyltransferase (Bv_23190_dwsd) were included in this cluster. The pathway analysis suggested these genes were abundant included in RNA transport and protein processing in endoplasmic reticulum.

Small RNA sequencing
To identify miRNAs that respond to salt-stress in B. vulgaris, ten small RNA libraries (half leaf and half root) were constructed from six-leaf-stage seedlings treated with NaCl for 0 h, 12 h, 24 h, 48 h, and 72 h and then sequenced. A total of 126.3 million reads with an average of 12.6 million reads for each sample were generated. After filtering the low-quality reads and removing the sequences belonging to mRNA, repeat sequence, and Rfam RNA tags, approximately 66 million reads with length range from 18 nt to 25 nt were selected as valid reads for further analysis (Table 1). The length distribution of the unique small RNA reads indicated that 24 nt (36.6%) sRNAs were the most abundant class, followed by 23 nt (18.6%), 22 nt (12.7%) and 21 nt (11.3%) sRNAs (Supplementary file 2: Figure S4 a). From the broader perspective of high-throughput sequencing of small RNAs from B. vulgaris, we observed that sRNAs of 24 nt dominated the library of unique species, in concurrence with previous reported for other plant species, such as Arabidopsis thaliana (Fahlgren 2007), Medicago truncatula (Wang 2011), and tea plants (Zhang 2014), that may implies the abundant representation of endogenous siRNAs in all the samples analyzed. In each sample, the miRNA expression levels at any two given stages were highly correlated, with the average correlation coefficient (r) above 0.97 (Supplementary file 2: Figure   S4 b). For instance, a high correlation was observed between root and leaf, despite their distinct difference on transcription level, and the sample dendrogram showed that miRNAs in control had a relatively far distance compared with treatment groups in either root or leaf.

Identification of known and novel miRNAs
There are different strategies for identifying miRNA genes from sequences depending on whether the target species has known miRNA genes. Since there are no B. vulgaris miRNAs reported in miRBase, all viridiplantae miRNA data were used as a reference to identify known miRNAs. The unannotated sRNA tags were mapped onto B. vulgaris genome with no mismatch to predict novel miRNAs, and all Mappable reads should be carefully evaluated for their precursor structure, including length, hairpinstructure, mismatch, (G + C) content, MFE, and miRNA* (Henderson 2006 andBerezikov 2006).
Identification of anti-sense miRNAs (miRNAs*) from the candidate miRNAs provides credible evidence that these miRNAs are authentic, and many of our candidate miRNAs had miRNAs* in their sequencing data. Finally, a total of 476 miRNA candidates including 219 known miRNAs (root: 148 and leaf: 196) and 257 novel miRNAs (root: 223 and leaf: 198) were identified in our sequencing results.
Once in Dhom's B. vulgaris genome sequencing research, 522 miRNAs has been predicted based on bioinformatics among which 258 were covered by their small RNA reads (Juliane 2014). Compared with their results, 158 miRNAs in our data were consistent with their sequencing results, other 8 miRNAs in our data covered their bioinformatics prediction and the remaining 310 miRNAs were newly identified by our sequencing.
The length distribution of identified miRNAs ranged from 18 to 25 nt with 21 nt (28.36%) and 24 nt (30.88%) being the most abundant (Additional file 2: Figure S4 c). The length between 21 and 24 nt account for the majority of candidates which are agreement with the typical characteristic attribute of DICER-LIKE (DCL) cleavage in plant (Zhixin Xie 2014). The nucleotide composition of the mature miRNAs was evaluated. Our analysis indicated that each high-abundance miRNA cluster (21 nt or 24 nt) has its distinct sequence bias, with a 5′-uridine predominating in the 21 nt cluster and a 5′adenosine in the 24 nt miRNA cluster (Additional file 2: Figure S4 d). Similar phenomena has been observed in many plant species and also accord with the specificity of DCL restriction sites. The average GC content of B. vulgaris miRNAs was about 45.9% which was similar with Medicago (44%), Arabidopsis (45%) and soybean (46%)(Additional file 1: Table S3). According to miRBase, the 219 conserved miRNAs were belonging to 40 families, among which MIR159 was the largest family with 25 members followed by MIR156 with 21 members, MIR5268, MIR6300 and MIR6478 tied for third with 16 members. The 257 novel miRNAs were grouped into 238 novel families based on their precursors.

Differentially expressed miRNAs under salt-stress
Global normalization was carried out to correct copy numbers among different samples to eliminate the influence sequencing discrepancy on the calculation of small RNA expression. Therefore the calculated gene expression can be directly used for comparing the difference of gene expression among samples. Our analysis indicated that most of identified miRNAs were expressed in more than one sample (Additional file 1: Table S3). The expression patterns of miRNAs under salt-stress (five time-points) in each tissues were evaluated by STEM to distinguish potential miRNAs response to saltstress. In the present study, 130 miRNAs including 60 in root and 82 in leaf showed significant changes in response to salt-stress. Cluster analysis was conducted to further elucidate the expression patterns of DE miRNAs in leaf and root, respectively (Additional file 2: Figure S5a and b). The majority of differentially expressed miRNAs (DE miRNAs) in root showed a trend of up-regulation under saltstress, especially at 12 h. While in leaf, half of the DE miRNAs presented a trend of down-regulation with salt-stress, and the other half showed a trend of increase from 0 h to 24 h and then decrease. A significant number of DE miRNAs were differentially expressed in tissues-specific manner. Under stress, around 37% of the differentially expressed miRNAs were tissues-specific in root and nearly 13 54% DE miRNAs in leaf. Interestingly, despite the expression pattern of miRNA is tissue-specific, however in different tissues, different miRNAs from the same family are often present. For instance, miR159k and miR166c-5p were differential expression under salt-stress in root, while miR159e and miR166a-3p were differential expression under salt-stress in leaf. In addition, most of the novel miRNAs showed low expression levels and some were only observed at one stage during exposure to salt-stress, but nov-miR6-3p, nov-miR51, nov-miR64 and nov-miR70 were expressed during all five stages in both root and leaf (Additional file 1: Table S3).

Target prediction of miRNAs via bioinformatics and degradome sequencing
Unlike animals, plant miRNAs usually have perfect or near-perfect complementarity with their targets to induce target gene splicing and thus regulate gene expression, and this feature can be used to predict the targets of plant miRNA. Using psRNATarget server, we could identify 2534 targets for a total of 317 (187 known and 130 novel) miRNAs (Additional file 1: Table S4). Accurate identification of target genes is essential to reveal the regulatory function of miRNAs. Since the prediction method cannot distinguish the authenticity of the targets, and in some cases like a high mismatch between miRNA and their targets, the prediction method may lose many real target genes. Thus, degradome sequencing were further used to identify the target genes of miRNA detected in this study. Two degradome libraries were constructed and sequenced from seedling of six-leaf-stage B. vulgaris for leaf and root, respectively. Ultimately, a total of 51.4 million raw reads were generated (Table 2) Based on the abundance of degradome tags at the target sites, these cleaved targets were classified into five categories (0 to 4). In brief, 439 and 317 non-redundant targets were identified in root and leaf, respectively. Interestingly, same cleavage site sometimes may belong to different category between two samples, for instance, miR159p were identified to cleave Bv4_074580_ygxg at the site of 1533 with "category 1" in root, while "category 0" in leaf. P-value in this study is not a screening parameter but for reference only as the reliability of target gene cleavage prediction. As can be seen 14 from the results, most P-values of "category 0" and "category 1" were less than 0.05, however, most P-values of "category 2", "category 3" and "category 4" were greater than 0.05 (Additional file 1: Table S4). Finally, 352 target genes for a total of 167 (167 known and 28 novel) miRNAs were identified through degradome sequencing. These cleavage sites were represented in the form of target plots (T-plots). Moreover, the difference in the target genes between the two libraries suggested that the cleavage of targets by miRNAs had organization specificity.
Through in silico and degradome analysis, a total of 2534 targets were identified for 317 miRNAs (Additional file 1: Table S4). It was seen that the maximum targets were obtained for members of MIR159 family (353), followed by miR156 (300) and MIR5368 (270). Functional annotation revealed that TFs and kinases represented a considerable fraction of miRNA targets. In total, 92 TFs belonging to 30 families were identified as targets. The maximum targets were from AP2 family (15.2%), followed by MYB (12.0%) and bHLH (8.7%). Besides, varieties of PPR protein genes, channel protein genes, oxidase genes, reductase genes, stress-responsive genes, disease resistance genes and phytohormones also occupied a large proportion among the targets. GO classification analysis was further performed to elucidate the potential role of miRNA targets in B. vulgaris. These targets were assigned to 460 biological processes, 514 cellular components and 334 molecular functions (Additional file 2: Figure S6). In addition, KEGG analysis was conducted to further elucidate the biological pathways of the miRNA target genes, which were classified into five groups including: "organismal system", metabolism", "genetic information processing", "environment information processes" and "cellular processes" (Additional file 2: Figure S7). The categories of "carbohydrate metabolism", "Folding, sorting and degradation", "Transcription" and "Environmental adaptation" were the most enriched pathways. Moreover, all the targets could be used for KOG classifications which were classified into 22 groups (Additional file 2: Figure S8). The most dominative group was the "Function unknown", followed by the "Transcription" and "Signal transduction mechanisms". This result is consistent with the functional annotation that transcription factors and kinases occupy a large part of targets. These results suggested that miRNAs may regulate transcription factors and signal transduction processes to take part in salt-stress response in B. vulgaris.
Further, to investigate the association of salt-stress responsive miRNAs with their targets, network analysis was performed using Cytoscape platform. For network construction, DE miRNAs which were considered as potential miRNAs response to salt stress and their targets involved in salt stress response were included. These stress-responsive genes included genes encoding for many Kinases like G-type lectin S-receptor-like serine/threonine-protein kinase (B120) and probable membraneassociated kinase regulator 2 (MAKR2), heat shock proteins (HSP), ABC transporter family members, detoxification enzymes like catalase and peroxidase, aquaporin and phytohormones-associated proteins like indole-3-acetic acid-amido synthetase and gibberellin 2-beta-dioxygenase. The network also included a lot of TFs like NAC, WRKY, ARF and so on (Figure 4).

Correlation analysis of miRNAs expression profiles and their target genes
To find out the mediatory role of miRNAs and their targets under salt stress in B. vulgaris, the expression pattern of both differential expression miRNAs (from small RNA-seq) and their differentially expressed target genes (from RNA-seq) were analyzed comprehensively. A total of 242 and 314 miRNA-mRNA interaction pairs were identified under stress conditions in root (Additional file 1: Figure S9a) and leaf (Additional file 1: Figure S9b), respectively. Among these interaction pairs, 209 positive correlation pairs which showed an opposite expression pattern between miRNA and mRNA were screened out for detailed analyzed since they were more consistent with the negative regulation of miRNA on mRNA. The rest negative correlation pairs were not focused here, nevertheless, it does not mean that there is no negative relationship between them, while the expression pattern of these target genes may be affected by complicated multi-factors. In root, 106 positive correlation pairs including 91 genes and 24 miRNAs were found (Figure 5a). Similarly, 103 coherent pairs consisting of 80 genes and 29 miRNAs were found in leaf (Figure 5b). It was seen that a higher number of miRNAs were up-regulated and their corresponding targets were down-regulated in root when compared with leaf. In root, miR396h was down-regulated and its targets nuclear transcription factor Y subunit A-3 (Bv3_063280_fahh) and catalase-like gene (Bv6u_155110_qtds) were up-regulated. Other miRNAs like miR159j, miR6300p, and a novel miRNA nov-miR46 were up-regulated and their targets ABC transporter B family member 11-like (Bv_20580_kwui), serine/threonine protein kinase PBS1-like (Bv2_028420_pmtd) and transcription factor TCP 21-like (Bv9_206510_huxy), respectively, were down-regulated. In leaf, about half miRNAs like miR156i, miR164e and a novel miRNA nov-miR217 were down-regulated and their targets two-component response regulator ARR2 (Bv1u_022840_ktmf), cyclic dof factor 1-like (Bv6_131470_guat), and eukaryotic translation initiation factor 3 subunit I-like (Bv2_046750_tipa), respectively, were up-regulated. The remaining half such as miR159n, miR2938a and miR2938b were up-regulated and their target genes heat shock 70 kDa protein 18-like (Bv1_010750_yzwe), pentatricopeptide repeat-containing protein (Bv1_007030_uxgf) and NB-ARC domain disease resistance protein (Bv8_184400_xkic), respectively, were down-regulated.
According to NR annotation and STRING protein interaction database (Christian 2005), PPI (proteinprotein interaction) analysis was carried out for DEGs, ultimately, a total of 213 potential salt stress responsive-associated genes with 18 confirmed miRNA-mRNA pairs were selected to construct a regulatory network of B. vulgaris in respond to salt-stress (Figure 7). Among these DEGs, nine KEGG pathway were significant enrichment including: Plant hormone signal transduction (bvg04075: 46 genes) and MAPK signaling pathway-plant (bvg04016: 20 genes) were marked in purple in figure 7, Porphyrin and chlorophyll metabolism (bvg00860: 17 genes), Carbon fixation in photosynthetic organisms (bvg00710: 10 genes) and Photosynthesis (bvg00195: 5 genes) were marked in green, Glutathione metabolism (bvg00480: 25 genes) and Ascorbate and aldarate metabolism (bvg00053: 5 genes) were marked in orange, Spliceosome (bvg03040: 10 genes) were marked in blue, Peroxisome (bvg04146: 9 genes) were marked in yellow. Based on KEGG classification, DEGs can be classified into four categories according to functions. The first category is genes involved in signal transduction like Auxin-responsive protein IAA29 (Bv8_198330_iozs) and Abscisic acid receptor PYL9 (Bv1_003450_qfuz), the second category is genes involved in transcriptional regulation like miR167d and its target ethylene-responsive transcription factor 014 (Bv1_000100_pnxw), the third category is genes that maintain the intensity of photosynthesis like Ribulose bisphosphate carboxylase small chain (Bv2_027420_aizd) and Photosystem II core complex proteins psbY (Bv6_146250_qqmp), and the last category is genes related to reactive oxygen species (ROS) scavenging system like Superoxide dismutase [Cu-Zn] 2 (Bv_01630_nkex) and Glutathione S-transferase U10 (Bv5_113320_zudp).

Quantitative RT-PCR validation of differential gene and miRNA expression
The expression patterns of genes and miRNAs obtained from RNA and small RNA sequencing were validated by the quantitative reverse transcription-polymerase chain reaction (qRT-PCR). The expression of fifteen randomly selected miRNAs (including 8 from root and 7 from leaf ) and twelve of their target genes were validated via qRT-PCR (Figure 8 and Additional file 2: Figure S9). Among of them, six novel miRNAs were also validated. Similar expression trends of miRNAs and target genes were observed between qRT-PCR and high-throughput sequencing. These results indicate the high reliability of high-throughput sequencing.

Discussion
High-throughput sequencing and the combination analysis of multi-omics provide an effective solution to explore the molecular mechanism of stress response in plant. In the present study, three highthroughput approaches including RNA-seq, small RNA sequencing and degradome sequencing were used for illuminating the genetic and molecular basis of salinity resistance in B. vulgaris. For plants, expose to salt-stress usually results in a lag phase of growth, which is followed by growth recovery after adapting to stress, albeit the observed growth vigour is often lower than control conditions (Magdalena 2015). When treated directly with 300 m mol NaCl, the seedlings of cultivated 'O68' wilted within 1 hour due to osmotic pressure changes (Additional file 2: Figure S10 a). This phenomenon recovered in about four hours, suggesting the plant has preliminary adapted to the osmotic pressure. After 12 h of treatment, some root tips were observed to turn black, indicating the toxic effect of salt on the root tips (Additional file 2: Figure S10 b). While salt will be transported from root to the above-ground tissue for storage. When the treatment time reaches 72h, the old leaves were observed to turn yellow (Additional file 2: Figure S10 a). Under this background, one control node (0 h), and four treatment nodes (12h, 24h, 48h and 72h) were set up to explore the salttolerance genes in B. vulgaris. As salt tolerance is a multigene trait, the identification of salt-related genes facilitates an understanding of the physiological processes and molecular mechanisms of plant responses to salt-stress. A total of 12 231 genes showed differential expression pattern under saltstress, which can be loosely classified into two broad categories: In the present study, nearly 500 protein kinases genes showed significant differential expression under salt-stress, including nine CDPK genes, sixteen CaM genes, twelve CBL genes, and eight annexin genes. PPI analysis of proteins encoded by these genes indicated that CDPK18 (Bv6_127310_prrx) interaction with calcium/calmodulin-regulated receptor-like kinase 1 (Bv_11660_smrf), and this calmodulin interaction with CBL-interacting serine/threonine-protein kinase 4 (Bv7_161470_jxag) and CBL-interacting serine/threonine-protein kinase 14 (Bv7_156600_urmk) together with another calmodulin-binding protein 60 A isoform X1 (Bv4_076600_wwmd) (Figure 9).

Possible function of transcription factors in beet salt-stress response.
Unlike animals, plants have to adapt to environmental changes continuously, they need to adjust their growth and processes of life timely in response to such alterations. Transcription factors are able to 20 alter the cell transcriptome through a very sophisticated network and play central roles in those complex adaptation processes. A total of 555 differentially expressed TF coding genes were identified by RNA-seq (Supplementary file 2: Figure S2), including some families that have been reported to participate in plant abiotic stress response like AP2/ERF, bZIP, WRKY, NAC, etc. Salt-responsive genes are often categorized into abscisic acid (ABA)-independent or ABA-dependent groups. ABAindependent genes usually contain dehydration responsive elements/C-repeats (DRE/CRTs) in their promoters and are regulated by AP2/ERF transcription factors. Our results showed that ethyleneresponsive transcription factor ERF073 gene (Bv_50650_ryod) was the most abundant among these DE transcription factor coding genes (Figure 7), with a fifteen-flod increase compared with control, implying that it may play an important role in regulating the expression of downstream salt-tolerance genes. In addition, two significantly up-regulated genes encoding for DREB1 (Bv3u_069010_ignp) and DREB2 (Bv2u_047070_drzk) were also identified under salt stress in B. vulgaris (Figure 7). Another category genes responsive to salt-stress require abscisic acid (ABA) mediation, these ABA-inducible genes usually contain a conserved ABA-responsive cis-acting element (ABRE) in their promoter regions. ABREs can be recognized by transcription factors containing a basic leucine zipper structure (bZIP). We identified a gene encoding bZIP transcription factor (Bv1_015900_smoy) that was significantly induced by salt stress in B. vulgaris (Figure 7). WRKY proteins are also important components of plant hormone signaling network in response to abiotic stimuli. Previous studies have shown that overexpression this family genes (eg. GmWRKY13, AtWRKY57, MusaWRKY71) confers tolerance to salt stress in transgenic plant by regulating plant growth, osmotic balance, Na+/K+ homeostasis, and the antioxidant system (Zhou 2008, Shekhawat 2011and Jiang 2012. Three homologous genes encoding for WRKY13 (Bv6_130210_rsgt), WRKY57 (Bv6_133740_qpam), WRKY71 (Bv1_004110_zfhn) were also significantly up-regulated in B. vulgaris during salt-stress (Figure 7). The NAC genes are specific and abundant to plants, which seem to play various roles not only in growth and development but also in the recognition of environmental stimuli. Our RNA-seq data revealed that three NAC domain-containing protein coding genes (Bv4_078440_hzcx, Bv3_051270_usim and Bv6_149120_hdfh) were highly expressed under salt stress (Figure 7).

Possible action of posttranscriptional control on beet salt-stress response.
A series of modifications to precursor mRNAs begin almost immediately with the initiation of transcription, that are essential for their nuclear export, maturation, and subsequent translation in eukaryotes (Tom 2018). Among these modifications, alternative splicing plays a key role in gene expression by tight contacting with transcription, translation, downstream mRNA metabolic events and even miRNA processing (Kang Yan 2012). Pre-mRNA splicing is performed by the spliceosome, a large dynamic protein complex composed of five small nuclear riboproteins (snRNPs) and a number of non-snRNP proteins, such as serine/arginine-rich (SR) proteins and glycine-rich RNA-binding proteins (GRPs). In plants, this post-transcriptional mechanism may save the time required for changes in transcriptional activation and pre-mRNA accumulation which allowing rapid adaptation to adverse environmental conditions. Recent studies have shown that salt stress changes the alternative splicing pattern of over 6000 arabidopsis genes, and a splicing factor, ski-interacting protein (SKIP), is report to confer osmotic tolerance under salt stress by controlling the alternative splicing of salt-tolerance genes (Jinlin F 2015). Our results identified significant differences in the expression of nineteen snRNP genes, thirteen SKIP genes, seven GRP genes and seven SR genes, which may function in alternative splicing of beet in response to salt stress (Figure 9).
The understanding of post-transcriptional regulation has been greatly broadened and improved since miRNAs were discovered over the past 20 years. Many stress responsive miRNAs have been identified by microarray and deep sequencing in various crop plants. In beet also, a few studies reporting the genome-wide discovery of miRNAs (Juliane C 2014). In this study, a total of 476 miRNAs were identified by our sequencing, of which 130 demonstrated differential expression patterns under saltstress. The function of miRNA can be inferred by their targets, and a very high number of targets (2534) were predicted for the identified miRNAs by using both computation and degradome sequencing. These targets included transcription factors (NAC, ERF, ARF, etc.), protein kinases (calcium-dependent protein kinase 1-like, G-type lectin S-receptor-like serine/threonine-protein kinase, wall-associated receptor kinase, etc.), ubiquitinated proteins (E3 ubiquitin-protein ligase ATL6like), hormone synthesis (jasmonic acid-amido synthetase JAR1-like), channel proteins (K + efflux 22 antiporter 3 and aquaporin PIP-type), and detoxifying enzymes (thioredoxin F-type, catalase and peroxidase), which were involved in various aspects of growth, development, signal transduction and stress response in B. vulgaris. The integrated analysis of DE miRNAs and mRNAs revealed several salt responsive miRNAs. For instance, a number of miR164 family members were seen to target the NACs (Figure 7), which is one of the largest families of transcriptional regulators in plants ( Xin Lu 2017).
Similarly, miR160d and miR396b-5p was found to target ARF and GRF (Figure 7), respectively, both of which may play an important role in balance of survival and growth timely under salt stress (Shimizu-Sato 2001).The up-regulation of various miRNAs like miR159n, miR159f, miR2592a, and nov-miR46 and the down-regulation of their targets, heat shock 70 kDa protein, E3 ubiquitin-protein ligase, Gtype lectin S-receptor-like serine/threonine-protein kinas and TCP21, respectively, suggests that miRNA regulated resistance might be nurtured during salt-stress. The degradome sequencing also confirmed the cleavage of some genes by their respective miRNAs. Recent advances in highthroughput sequencing methods have revolutionized the identification of low-abundance, novel miRNAs in various species. In this study, many novel miRNAs with low abundance were discovered from various stress stages. The expression profiles of these novel potential miRNAs during the different stages might provide valuable information on the response to salt-stress in B. vulgaris.
Dicer-like (DCL) proteins are double-stranded RNA (dsRNA)-specific endoribonuclease, play a key role in either double-stranded RNAs cleavage or single-stranded RNAs which bearing stem-loop structures (such as miRNA precursor) cleavage (Akihito 2017). Seven DE DCL genes were found by RNA-seq in this study, in addition, miR162-3p was found to target one of them, a Dicer homolog 1 gene (Bv1_001950_imjw). Argonaute proteins are direct binding to small RNAs and essential players in all small RNAs guided gene-silencing processes (G Meister 2013). Six differential expression Argonaute genes were existed in our results, interestingly, miR168b was found to target one of them, an argonaute 1 gene (Bv8_192000_ctnt). Moreover, recent research suggests that Argonaute proteins are also related to chromatin modification and alternative splicing [ ! ] , the above results indicate that there are interactions among several posttranscriptional regulatory pathways, which constitute complex posttranscriptional mechanism under salt-stress.
Possible effect of functional genes in beet salt-stress response.
The first step in adapting to salt stress is to solve the osmotic stress caused by high salinity in the environment, glycine betaine plays a crucial role in this process as intracellular osmoprotectants.

Plant materials and RNA extraction
The cultivar sugar beet "O68" were from our own laboratory (Heilongjiang, China). Let the seeds soaked in water for 10h, then sterilized in 0.1% (v/v) HgCl 2 for 10 min, washed repeatedly with distilled water, and germinated on wet filter paper in germination box at 26 °C for 2 days. After 25 germination, seedlings were selected, transferred to plastic pots (43.5 cm × 20 cm × 14 cm, 10plants per pot) filled with quarter-strength Hoagland solution. The germinating seeds were cultivated under 14/10 light photoperiod at 24 °C (day)/18 °C (night) in a phytotron (Friocell 707, Germany). Six-leafstage seedlings were treated with 300mM NaCl for 12, 24, 48 and 72h and untreated were used as controls. The experiment included 5 groups with 20 plants for per group. Finally, 10 uniform plants were selected from the 20 plants in each group (Additional file 2: Figure S10). The samples were collected from control and salt treatment leaves and roots as shown in Figure S10, respectively. And the samples were immediately frozen in liquid nitrogen and stored at -80 °C until further use.
Total RNA was extracted from mixed leaf tissues and root tissues by TRIZOL (Invitrogen, USA) according to the manufacturer's instructions. A total of 0.5 g leaves or roots from 10 plants were mixed for each sample. The qualitative and quantitative assessment of these total RNA samples were conducted using Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit). The total RNA from each tissue was used for library construction for RNA-seq, small RNA and degradome sequencing. The library construction and sequencing were conducted by LC Sciences, Hangzhou, China.

RNA-seq and data analysis
The high-quality total RNA (RIN 8) was used for RNA-seq library construction using TruSeq RNA finally sequenced by Illumina HiSeq2500 platforms. 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 B. vulgaris miRNA was unrecorded in miRBase 21.0 (http://www.mirbase.org), the filtered unique reads with length 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 (http://bvseq.molgen.mpg.de) to determine their genomic locations using Bowtie [ ! ] . 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 predicated from the flank 120 nt sequences using the mfold software (Zuker 2003).
A global normalization method was used to correct copy numbers among different samples followed the procedures as described in a previous study (Xiang Li 2016).This method is able to eliminate the influence sequencing discrepancy on the calculation of small RNA expression. Therefore the calculated gene expression can be directly used for comparing the difference of gene expression among samples. The DE miRNAs were identified followed the same procedures as mentioned in "RNAseq data analysis" based on normalized counts. Target Finder was used to identify putative targets for all identified miRNA candidates.

Degradome sequencing and analysis
For degradome library construction, Approximately 20 μg of total RNA from all five time points were pooled together to generate a mix-library for leaf and root, respectively. The library construction was performed as previously described by Zhaorong Ma (2010)  remove low-quality reads, reads with N s and any reads with adaptor and primer contamination using ACGT101-DEG (Lc-bio Sciences, Houston, Texas, USA). The clean reads were used to identify the degraded fragments of mRNA after removed the reads which can be annotated into different ncRNA database. The potential miRNA editing sites were identified use the small RNA sequencing data by CleaveLand 3.0 (Addo-Quaye 2009).The identified sites were categorized into five categories (0, 1, 2, 3 and 4) based on the read abundance at that position as previous study (Xu, X 2013). Based on the signatures and abundances along the B. vulgaris RNA-seq data, t-plots were built for the high efficiency analysis of the potential miRNA targets. Finally, all of the identified potential target genes were annotated using the information from NCBI. The functional classification and pathways based on the three database: KOG (Clusters of Orthologous Groups of proteins), KO (KEGG Ortholog database) and GO (Gene Ontology).

Real-time quantitative Reverse transcription PCR (RT-qPCR) analysis
For validation of the gene and miRNA expression obtained from high-throughput sequencing, qRT-PCR of randomly selected genes and miRNAs was performed. The qRT-PCR reactions were performed using the High-Capacity cDNA Reverse Transcription Kits (Thermo Fisher SCIENTIFIC, CA) and iTaq Universal SYBR® Green Supermix (BIO-RAD, Hercules, CA) on the CFX Real-time PCR system (BIO-RAD, CA). All the primers were listed in Additional file 1: Table S5. 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 (Kramer 2011). The U6 gene was used as an endogenous control.
For mRNA, PP2A+ UBQ5 and PP2A+25S were used as endogenous controls for root and leaf, respectively, which were selected from 15 candidate genes by evaluated using Bestkeeper (Pfaffl 2004), NormFinder (Andersen 2004) and GeNorm (Vandesompele 2002). To avoid non-specific amplification, melting curve was carried out for each PCR product. The expression level of the miRNAs in different samples were calculated by comparative 2 −△△ CT method.
Declarations Figure 2 Annotation and analysis of differentially expressed genes (DEGs).
38 Figure 3 Expression patterns of differentially expressed genes (DEGs) in B. vulgaris under salt stress.
39 Figure 4 A network representing the relationships between miRNAs and their targets associated with salt-stress.
40 Figure 5 Expression profiles of important miRNA-mRNA interaction pairs in root and leaf.
42 Figure 7 Overview of important genes in response to salt stress in B. vulgaris.
43 Figure 8 Comparison between qRT-PCR and deep sequencing miRNA expressions during the five stages of salt-stress.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download. Table1.xlsx