Integration of mRNA and miRNA Analysis Reveals the Molecular Mechanism Underlying Salt and Alkali Stress Tolerance in Tobacco

Salinity is one of the most severe forms of abiotic stress and affects crop yields worldwide. Plants respond to salinity stress via a sophisticated mechanism at the physiological, transcriptional and metabolic levels. However, the molecular regulatory networks involved in salt and alkali tolerance have not yet been elucidated. We developed an RNA-seq technique to perform mRNA and small RNA (sRNA) sequencing of plants under salt (NaCl) and alkali (NaHCO3) stress in tobacco. Overall, 8064 differentially expressed genes (DEGs) and 33 differentially expressed microRNAs (DE miRNAs) were identified in response to salt and alkali stress. A total of 1578 overlapping DEGs, which exhibit the same expression patterns and are involved in ion channel, aquaporin (AQP) and antioxidant activities, were identified. Furthermore, genes involved in several biological processes, such as “photosynthesis” and “starch and sucrose metabolism,” were specifically enriched under NaHCO3 treatment. We also identified 15 and 22 miRNAs that were differentially expressed in response to NaCl and NaHCO3, respectively. Analysis of inverse correlations between miRNAs and target mRNAs revealed 26 mRNA-miRNA interactions under NaCl treatment and 139 mRNA-miRNA interactions under NaHCO3 treatment. This study provides new insights into the molecular mechanisms underlying the response of tobacco to salinity stress.


Introduction
Salinity stress is one of the most constraining environmental factors that threatens both the productivity and quality of crop plants worldwide [1]. In China, a large amount of land is unsuitable for agricultural production because of increasing soil alkalinity. Terrestrial soils affected by salt are divided mainly into saline soils, alkaline soils and salt-alkaline soils [2]. NaCl is a major neutral component that causes salt stress, and NaHCO 3 and Na 2 CO 3 play a key role in soil alkalization with increasing pH. By causing osmotic stress, ion toxicity, redox imbalance and metabolic damage, soil salinity and salinization can severely affect biological homeostasis in plants [3,4]. Plants have evolved strategies to cope with salt stress. These strategies include increased accumulation of compatible solutes such as proline, soluble sugars and betaine to adjust to osmotic injury [5]; activation of reactive oxygen species (ROS) scavenging systems to protect plants against oxidative damage triggered by salinity [6]; and regulation of the expression of salt-responsive genes to re-establish ion homeostasis

Differentially Expressed Genes under Salt-Alkali Stress
To identify genes that are differentially regulated under various salt-alkali conditions, we used the DESeq software to perform a differential expression analysis between the treatment and CK samples. With strict criteria of |log2FC|> 2 and p-value less than 0.05 (|log2FC|> 2, p < 0.05), a total of 2206 genes (Table S1) and 7468 genes (Table S2) were found to be differentially expressed in response to SS and AS, respectively ( Figure 2A). There were many more differentially expressed genes (DEGs) in the AS treatment than in the SS treatment, indicating that NaHCO3 treatment induced a significantly greater number of biological and metabolic changes than did the NaCl treatment. To better characterize the salt-alkali stress-induced changes in gene expression in tobacco, the overlapping DEGs between the SS vs. CK and AS vs. CK comparable groups are shown as a Venn diagram in Figure 2B. Detailed analysis of these 1610 overlapping DEGs revealed that 1578 of them exhibited the same trends under salt-alkali stress: 828 genes were upregulated and 750 genes were downregulated in both the SS vs. CK and AS vs. CK comparisons.

Differentially Expressed Genes under Salt-Alkali Stress
To identify genes that are differentially regulated under various salt-alkali conditions, we used the DESeq software to perform a differential expression analysis between the treatment and CK samples. With strict criteria of |log 2 FC|> 2 and p-value less than 0.05 (|log 2 FC|> 2, p < 0.05), a total of 2206 genes (Table S1) and 7468 genes (Table S2) were found to be differentially expressed in response to SS and AS, respectively ( Figure 2A). There were many more differentially expressed genes (DEGs) in the AS treatment than in the SS treatment, indicating that NaHCO 3 treatment induced a significantly greater number of biological and metabolic changes than did the NaCl treatment. To better characterize the salt-alkali stress-induced changes in gene expression in tobacco, the overlapping DEGs between the SS vs. CK and AS vs. CK comparable groups are shown as a Venn diagram in Figure 2B. Detailed analysis of these 1610 overlapping DEGs revealed that 1578 of them exhibited the same trends under salt-alkali stress: 828 genes were upregulated and 750 genes were downregulated in both the SS vs. CK and AS vs. CK comparisons.

Functional Classification of Salt-Alkali Responsive Genes
To better understand the functions of the DEGs, we mapped all of these salt-alkali responsive genes to the Gene Ontology (GO) database. In total, we identified 92 and 206 significantly enriched GO terms (false discovery rate (FDR) < 0.05) in the SS and AS treatments, respectively (Table S3). Among these GO terms, some crucial biological processes related to "response to water (GO: 0009415)," "antioxidant activity (GO: 0016209)" and "oxidoreductase activity (GO: 0016491)" were observed in both the SS and AS treatments, while genes involved in "photosystem (GO: 0009521)," "cellular carbohydrate metabolic process (GO: 0044262)," "glucan metabolic process (GO: 0006073)" and "defense response (GO: 0006952)" were enriched in only the AS treatment.
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis provided information concerning gene function. A total of 13 and 24 pathways were enriched in the SS and AS treatments, respectively ( Figure 3). The enriched KEGG pathways revealed that the DEGs detected under saltalkali stress were enriched in "photosynthesis-antenna proteins," "plant hormone signal transduction" and "glucosinolate biosynthesis". Fourteen KEGG pathways, including "photosynthesis," "starch and sucrose metabolism," "glutathione metabolism" and "carbon fixation in photosynthetic organisms," were specifically enriched in the AS treatment. These results suggest that the DEGs identified in this study may play important roles in salt and alkali stress in tobacco.

Functional Classification of Salt-Alkali Responsive Genes
To better understand the functions of the DEGs, we mapped all of these salt-alkali responsive genes to the Gene Ontology (GO) database. In total, we identified 92 and 206 significantly enriched GO terms (false discovery rate (FDR) < 0.05) in the SS and AS treatments, respectively (Table S3). Among these GO terms, some crucial biological processes related to "response to water (GO: 0009415)," "antioxidant activity (GO: 0016209)" and "oxidoreductase activity (GO: 0016491)" were observed in both the SS and AS treatments, while genes involved in "photosystem (GO: 0009521)," "cellular carbohydrate metabolic process (GO: 0044262)," "glucan metabolic process (GO: 0006073)" and "defense response (GO: 0006952)" were enriched in only the AS treatment.
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis provided information concerning gene function. A total of 13 and 24 pathways were enriched in the SS and AS treatments, respectively ( Figure 3). The enriched KEGG pathways revealed that the DEGs detected under salt-alkali stress were enriched in "photosynthesis-antenna proteins," "plant hormone signal transduction" and "glucosinolate biosynthesis". Fourteen KEGG pathways, including "photosynthesis," "starch and sucrose metabolism," "glutathione metabolism" and "carbon fixation in photosynthetic organisms," were specifically enriched in the AS treatment. These results suggest that the DEGs identified in this study may play important roles in salt and alkali stress in tobacco.

miRNA Sequencing Data Analysis
Nine small RNA (sRNA) libraries were constructed with three biological replicates for the CK, SS and AS treatments. A total of 44,642,448, 39,911,768 and 42,157,754 raw reads were obtained from the CK, SS and AS samples, respectively. After the adapter sequences, poly-A sequences, low-quality reads and reads shorter than 18 nt were discarded, 43,327,623, 38,719,988, and 40,868,990 clean reads were obtained from the CK, SS and AS samples, respectively. A majority of the sRNAs were 21-24 nt in length in all nine libraries, with 24 nt being the most frequent length ( Figure S1). On the basis of the method described above, a total of 88 mature miRNAs and 102 novel miRNAs were identified in these nine sRNA libraries (Table S4). Details concerning the distribution of the sRNAs are provided in Table 2.

miRNA Sequencing Data Analysis
Nine small RNA (sRNA) libraries were constructed with three biological replicates for the CK, SS and AS treatments. A total of 44,642,448, 39,911,768 and 42,157,754 raw reads were obtained from the CK, SS and AS samples, respectively. After the adapter sequences, poly-A sequences, low-quality reads and reads shorter than 18 nt were discarded, 43,327,623, 38,719,988, and 40,868,990 clean reads were obtained from the CK, SS and AS samples, respectively. A majority of the sRNAs were 21-24 nt in length in all nine libraries, with 24 nt being the most frequent length ( Figure S1). On the basis of the method described above, a total of 88 mature miRNAs and 102 novel miRNAs were identified in these nine sRNA libraries (Table S4). Details concerning the distribution of the sRNAs are provided in Table 2.

Differential Expression of miRNAs in Response to Salt-Alkali Stress
With respect to the miRNA-seq data analysis, a total of 33 differentially expressed miRNAs (DE miRNAs) were detected in the SS and AS treatments with the criteria |log 2 FC|> 0.5 and p-value < 0.05, 16 of which were identified as novel miRNAs (Table 3). Five conserved and 10 novel miRNAs were significantly differentially expressed under SS conditions, four of whose expression was upregulated and 11 of whose was downregulated. Fourteen conserved and eight novel miRNAs were differentially expressed in response to AS, and among those miRNAs, compared with that of the CK, the expression of 13 and nine was upregulated and downregulated, respectively. Of these 33 DE miRNAs, four, namely, nta-miR156a, nta-miR6149a, miR-novel_64 and miR-novel_216, that overlapped between the SS and AS samples were identified. As expected, with the exception of nta-miR156a, all of these overlapping DE miRNAs exhibited similar regulation of expression under the SS and AS treatments. Table 3. List of differentially expressed microRNAs (DE miRNAs) in response to SS and AS treatments.

Integration Analysis of the miRNAs and mRNAs
The interaction between DE miRNAs and their targets was investigated using Cystoscape (version 3.6.1). A total of 165 miRNA-mRNA interactions were identified among the SS and AS treatments, with the involvement of 24 DE miRNAs and 158 DE mRNAs ( Figure 4, Table S5). Figure 4 shows that a single miRNA can regulate multiple mRNAs and that a single mRNA can be targeted by more than one miRNA. These results suggest that the miRNA and mRNA interaction network involved in salt and alkali stress is highly complex, especially for different miRNAs with their targets under various conditions. The interaction between DE miRNAs and their targets was investigated using Cystoscape (version 3.6.1). A total of 165 miRNA-mRNA interactions were identified among the SS and AS treatments, with the involvement of 24 DE miRNAs and 158 DE mRNAs ( Figure 4, Table S5). Figure  4 shows that a single miRNA can regulate multiple mRNAs and that a single mRNA can be targeted by more than one miRNA. These results suggest that the miRNA and mRNA interaction network involved in salt and alkali stress is highly complex, especially for different miRNAs with their targets under various conditions.

RT-qPCR Validation of DE mRNAs and miRNAs
To validate the mRNA and miRNA sequencing data, the expression levels of four DE miRNAs (nta-miR156a, nta-miR6149a, nta-miR394, nta-miR398) and eight DEGs (CHX 15, AKT 1, ABCG 11, PIP 2-7, SOD, AMY, SUS2, RPP 13), which were determined to be important for salt and alkali stress in our study, were manually selected and analyzed by RT-qPCR ( Figure 5). As expected, with the exception of that of RPP 13 in the SS treatment, the expression patterns of most of the examined miRNAs and genes were similar, showing the accuracy and reliability of the RNA-seq data.

RT-qPCR Validation of DE mRNAs and miRNAs
To validate the mRNA and miRNA sequencing data, the expression levels of four DE miRNAs (nta-miR156a, nta-miR6149a, nta-miR394, nta-miR398) and eight DEGs (CHX 15, AKT 1, ABCG 11, PIP 2-7, SOD, AMY, SUS2, RPP 13), which were determined to be important for salt and alkali stress in our study, were manually selected and analyzed by RT-qPCR ( Figure 5). As expected, with the exception of that of RPP 13 in the SS treatment, the expression patterns of most of the examined miRNAs and genes were similar, showing the accuracy and reliability of the RNA-seq data.

Figure 5.
RT-qPCR analysis of miRNAs and mRNAs. The miRNAs and mRNAs were isolated from tobacco treated with NaCl (SS) and NaHCO3 (AS), respectively. The 5.8S rRNA and ef1α genes were chosen as the internal controls for miRNA and mRNA, respectively.

Discussion
As one of the most threatening environmental stress factors, salt-alkali stress can hinder plant growth and development in several ways. Accelerated NGS technology enables the investigation of a wide range of physiological and metabolic mechanisms of stress responses at the molecular level. In this study, we first conducted an integrative analysis of mRNA and miRNA expression levels in tobacco under salt and alkali stress. From this integrative analysis, we obtained a valuable set of saltalkali-responsive mRNAs and miRNAs and identified the interactions and potential roles of these mRNAs and miRNAs in biophysiological processes between NaCl and NaHCO3 treatments. These data provide a comprehensive understanding of the differences between salt and alkali stress response mechanisms.

mRNA Profiling in Leaves of Tobacco under Salt-Alkali Stress
Under salt-alkali stress, plants can accumulate essential substances for osmotic adjustment and can regulate gene expression levels to maintain ionic homeostasis [38]. Cation/H + antiporters (CHXs) operate as Na + /H + antiporters or K + /H + antiporters that help mediate biological responses to excess ions. For example, AtCHX17 is strongly induced by salt stress and is involved in K + acquisition [39]. Moreover, an appropriate K + :Na + ratio in the cytoplasm is important for salt tolerance in plants.
Overexpression of the potassium channel gene OsKAT1 has been shown to enhance rice Na + tolerance by regulating K + /Na + homeostasis [40]. In our study, the expression of three CHXs and AKT 1 significantly changed under SS and AS treatments ( Table 4), suggesting that these ion transporters are responsible for tobacco resistance to salt toxicity via enhanced efflux of Na + and absorption of K + . In addition, ABC transporters use energy from ATP hydrolysis for nutrient uptake and play an important role in ion homeostasis [41]. In our study, an increased abundance of ABCG11 but suppressed expression of ABCA7 were identified under both the SS and AS treatments (Table 4). Panikashvili et al. demonstrated that AtWBC11 (AtABCG11) is induced by salt and that wbc11 mutants are highly sensitive to salt stress [42], suggesting a potential role for this protein in plant salt tolerance. Figure 5. RT-qPCR analysis of miRNAs and mRNAs. The miRNAs and mRNAs were isolated from tobacco treated with NaCl (SS) and NaHCO3 (AS), respectively. The 5.8S rRNA and ef1α genes were chosen as the internal controls for miRNA and mRNA, respectively.

Discussion
As one of the most threatening environmental stress factors, salt-alkali stress can hinder plant growth and development in several ways. Accelerated NGS technology enables the investigation of a wide range of physiological and metabolic mechanisms of stress responses at the molecular level. In this study, we first conducted an integrative analysis of mRNA and miRNA expression levels in tobacco under salt and alkali stress. From this integrative analysis, we obtained a valuable set of salt-alkali-responsive mRNAs and miRNAs and identified the interactions and potential roles of these mRNAs and miRNAs in biophysiological processes between NaCl and NaHCO 3 treatments. These data provide a comprehensive understanding of the differences between salt and alkali stress response mechanisms.

mRNA Profiling in Leaves of Tobacco under Salt-Alkali Stress
Under salt-alkali stress, plants can accumulate essential substances for osmotic adjustment and can regulate gene expression levels to maintain ionic homeostasis [38]. Cation/H + antiporters (CHXs) operate as Na + /H + antiporters or K + /H + antiporters that help mediate biological responses to excess ions. For example, AtCHX17 is strongly induced by salt stress and is involved in K + acquisition [39]. Moreover, an appropriate K + :Na + ratio in the cytoplasm is important for salt tolerance in plants.
Overexpression of the potassium channel gene OsKAT1 has been shown to enhance rice Na + tolerance by regulating K + /Na + homeostasis [40]. In our study, the expression of three CHXs and AKT 1 significantly changed under SS and AS treatments ( Table 4), suggesting that these ion transporters are responsible for tobacco resistance to salt toxicity via enhanced efflux of Na + and absorption of K + . In addition, ABC transporters use energy from ATP hydrolysis for nutrient uptake and play an important role in ion homeostasis [41]. In our study, an increased abundance of ABCG11 but suppressed expression of ABCA7 were identified under both the SS and AS treatments (Table 4). Panikashvili et al. demonstrated that AtWBC11 (AtABCG11) is induced by salt and that wbc11 mutants are highly sensitive to salt stress [42], suggesting a potential role for this protein in plant salt tolerance. Aquaporins (AQPs) are integral membrane proteins that regulate the transport of water and small molecules. In Arabidopsis and Beta vulgaris, the transcript abundance of most of the AQPs decreased in response to salt [43,44]. Hu et al. demonstrated that overexpression of the wheat AQP gene TaAQP8 increased salt stress tolerance in tobacco [45]. The expression levels of eight AQP genes in this study were consistent with previous results (Table 4), indicating that salt-alkali stress disrupts the balance of water transport efficiency [46].
Plant exposure to salt-alkali stress often generates ROS, which cause lipid peroxidation and DNA damage [47]. Plants can mediate the removal of ROS via antioxidant defense enzymes, including superoxide dismutase (SOD), catalase (CAT), peroxidase (POD) and glutathione S-transferase (GST) [48]. Our GO enrichment results indicated that DEGs were significantly enriched in the "antioxidant activity" term under both the SS and AS treatments. As expected, our experimental data showed that the expression of most of the ROS-scavenging enzymes, including SOD, CAT, POD, glutathione peroxidase (GR) and GST, was upregulated under salt-alkali treatment ( Table 4), indicating that these scavenging enzymes have active functions that provide resistance to salt and alkali stress in tobacco.
Photosynthesis is a well-studied primary metabolic process that is affected by environmental stress. The response of plant photosynthesis to salinity is typically accompanied by stomal closure, reduced CO 2 diffusion into chloroplasts and alterations in leaf photochemistry [6]. KEGG pathway analysis indicated that DEGs were significantly enriched in "photosynthesis-antenna proteins" under salt-alkali treatment but were enriched in only "photosynthesis" in AS-treated samples. As shown in Figure 6, the transcript abundance of genes involved in the photosynthetic process decreased. Under salt-alkali stress, DEGs encoding light-harvesting complexes, such as Lhca in photosystem I (PSI) and Lhcb in photosystem II (PSII), was downregulated ( Figure 6A), suggesting that the capture of light energy was restricted by the SS and AS treatments. Photo-oxidation of water in plants is accomplished by the oxygen-evolving complex of PSII. The decreased abundance of the genes encoding this complex, including PsbO, PsbP, PsbQ, PsbS, PsbY and Psb27, may inhibit the oxygen-evolving capacity of PSII for water oxidation ( Figure 6B). Furthermore, the similar expression profiles of DEGs encoding electron transport proteins and ATP synthase indicated an inhibitory effect on the mechanism of electron transport and resulted in the inhibition of ATP synthesis under AS treatment. Taken together, our experimental data revealed that NaHCO 3 stress strongly affects the process of photosynthesis because of the marked decrease in the expression levels of a wide range of relevant genes. Soluble carbohydrates can promptly accumulate in plants to cope with various forms of abiotic stress. Starch is the predominant carbohydrate storage material that is synthesized by plants. Our KEGG pathway enrichment analysis revealed that the DEGs were markedly enriched in the "starch and sucrose metabolism" pathway under AS treatment but not under SS treatment. Granule-bound starch synthase (GBSS) catalyzes the synthesis of starch, while α-amylase (AMY) and β-amylase (BMY) are responsible for degrading plant starch. Numerous studies have revealed significant changes in the expression levels of GBSS, AMY and BMY under various forms of abiotic stress [49,50]. Similarly, the expression of GBSS 1 was downregulated by NaHCO3 treatment, but that of AMY and Soluble carbohydrates can promptly accumulate in plants to cope with various forms of abiotic stress. Starch is the predominant carbohydrate storage material that is synthesized by plants. Our KEGG pathway enrichment analysis revealed that the DEGs were markedly enriched in the "starch and sucrose metabolism" pathway under AS treatment but not under SS treatment. Granule-bound starch synthase (GBSS) catalyzes the synthesis of starch, while α-amylase (AMY) and β-amylase (BMY) are responsible for degrading plant starch. Numerous studies have revealed significant changes in the expression levels of GBSS, AMY and BMY under various forms of abiotic stress [49,50]. Similarly, the expression of GBSS 1 was downregulated by NaHCO 3 treatment, but that of AMY and BMY was upregulated (Figure 7). Zhang et al. reported that the transcript abundance of genes encoding sucrose synthase (Sus) increased under high-pH stress (NaOH treatment). Consistent with this observation, the induction of Sus2 genes by AS in our study may be responsible for high-pH-induced oxygen deficiency conditions, as Sus2 was induced specifically by oxygen deficiency [51]. Trehalose is a non-reducing disaccharide that serves as an energy source for plants subjected to stress [52]. Trehalose-6-phosphate synthase (TPS) is a key enzyme involved in the biosynthesis of trehalose. In our study, we found that the expression of three genes encoding the TPS enzyme was induced and that the expression of only one TPS was downregulated by AS (Figure 7). Taken together, the results show that enhancement of starch and sucrose metabolism may provide energy requirements for tobacco resistance to NaHCO 3 -induced osmotic stress. this observation, the induction of Sus2 genes by AS in our study may be responsible for high-pHinduced oxygen deficiency conditions, as Sus2 was induced specifically by oxygen deficiency [51].
Trehalose is a non-reducing disaccharide that serves as an energy source for plants subjected to stress [52]. Trehalose-6-phosphate synthase (TPS) is a key enzyme involved in the biosynthesis of trehalose. In our study, we found that the expression of three genes encoding the TPS enzyme was induced and that the expression of only one TPS was downregulated by AS (Figure 7). Taken together, the results show that enhancement of starch and sucrose metabolism may provide energy requirements for tobacco resistance to NaHCO3-induced osmotic stress.

miRNA Profiling of the Salt-Alkali Response in Tobacco
Numerous recent studies have demonstrated that miRNAs play crucial roles in plant growth and development and in the plant responses to adverse environmental stress. In our study, we performed sRNA sequencing on the CK, SS and AS samples and identified 89 known miRNAs from 52 families and 102 novel miRNAs. Consistent with the results of a previous sRNA sequencing study in tobacco [33], our results showed that the most abundant transcript was that of nta-miR166a, followed by nta-miR159, nta-miR6149a and nta-miR6021, across the nine libraries (Table S4).
Many miRNAs are closely associated with salt-alkali stress in various plant species, such as Arabidopsis [53], radish [54], flax [55] and Medicago truncatula [56]. In the present study, 17 known miRNAs and 16 novel miRNAs were found to respond to salt-alkali stress in tobacco. We detected significant regulation of the miR156, miR482, miR6145 and miR6149 families under both salt and alkali stress. The conserved miR156 family is one of the most well-studied plant miRNA families and plays a vital role in plant stress responses, including salt-alkali resistance in alfalfa [57], cotton [58] and Populus euphratica [59]. The roles of miR482 in drought and cold stress, miR6145 in the wound response and miR6149 in Cr stress were also investigated [60][61][62][63]. Our experimental data confirmed the significance of miR156 in tobacco under salt-alkali conditions and provide valuable information concerning the roles of miR482, miR6145 and miR6149 in the response to salt-alkali stress. Notably, two of these overlapping miRNA families, namely, miR156 (nta-miR156a) and miR6145 (nta-miR6145a/e), exhibited differential regulation upon exposure to NaCl and NaHCO3 stress in our experiment. A similar result for miR156 was obtained in a previous study, which revealed that the expression of miR156 was downregulated after 4 h [64] but upregulated after 72 h of salt treatment in Medicago [56].
In addition to these four overlapping miRNA families, eight upregulated miRNAs (including

miRNA Profiling of the Salt-Alkali Response in Tobacco
Numerous recent studies have demonstrated that miRNAs play crucial roles in plant growth and development and in the plant responses to adverse environmental stress. In our study, we performed sRNA sequencing on the CK, SS and AS samples and identified 89 known miRNAs from 52 families and 102 novel miRNAs. Consistent with the results of a previous sRNA sequencing study in tobacco [33], our results showed that the most abundant transcript was that of nta-miR166a, followed by nta-miR159, nta-miR6149a and nta-miR6021, across the nine libraries (Table S4).
Many miRNAs are closely associated with salt-alkali stress in various plant species, such as Arabidopsis [53], radish [54], flax [55] and Medicago truncatula [56]. In the present study, 17 known miRNAs and 16 novel miRNAs were found to respond to salt-alkali stress in tobacco. We detected significant regulation of the miR156, miR482, miR6145 and miR6149 families under both salt and alkali stress. The conserved miR156 family is one of the most well-studied plant miRNA families and plays a vital role in plant stress responses, including salt-alkali resistance in alfalfa [57], cotton [58] and Populus euphratica [59]. The roles of miR482 in drought and cold stress, miR6145 in the wound response and miR6149 in Cr stress were also investigated [60][61][62][63]. Our experimental data confirmed the significance of miR156 in tobacco under salt-alkali conditions and provide valuable information concerning the roles of miR482, miR6145 and miR6149 in the response to salt-alkali stress. Notably, two of these overlapping miRNA families, namely, miR156 (nta-miR156a) and miR6145 (nta-miR6145a/e), exhibited differential regulation upon exposure to NaCl and NaHCO 3 stress in our experiment. A similar result for miR156 was obtained in a previous study, which revealed that the expression of miR156 was downregulated after 4 h [64] but upregulated after 72 h of salt treatment in Medicago [56].
In addition to these four overlapping miRNA families, eight upregulated miRNAs (including nta-miR164a, nta-miR168a, nta-miR168d, nta-miR171c, nta-miR394) and two downregulated miRNAs (nta-miR398, nta-miR408) were markedly expressed under NaHCO 3 stress, and nta-miR396a was specifically upregulated by NaCl. Liu et al. reported that miR168, miR171 and miR394 were induced and that miR398 was suppressed by salt stress in Arabidopsis [53]. Sun et al. demonstrated that the expression of miR164a and miR396a-3p was upregulated and that the expression of miR398b-3p was downregulated under salt stress in Raphanus sativus L. [54]. As miRNAs are evolutionarily conserved between plant species, our work here indicated the occurrence of miRNA regulation in tobacco under salt-alkali stress. Moreover, 16 novel salt-alkali-responsive miRNAs were identified, and the structural precursors of these novel miRNAs are presented in Figure S2. However, further validation is needed to confirm the authenticity of such novel miRNAs accurately.
In general, miRNAs negatively regulate gene transcription levels by target RNA cleavage based on sequence complementarity with targets. Our integrated analysis of miRNA and mRNA transcript levels revealed inverse regulation of 26 and 139 miRNA-mRNA interactions in the SS and AS treatments, respectively. For instance, the miR156 target squamosa promoter-binding protein like (SPL), which plays a role in flower regulation and reproductive growth [65], was also identified in our study. Recently, a role for the miR156/SPL network in the regulation of the content of ROS and salicylic acid (SA) signal transduction was identified in Arabidopsis [66]. As expected, nt-miR156a-targeted SPL members exhibited strict negative regulation in both the SS and AS treatments, suggesting that nt-miR156a plays an important role in salt-alkali stress by regulating the expression of SPLs. Nucleotide-binding sites and leucine-rich repeat (NBS-LRR) regions constitute the main components of disease resistance proteins [67]. These NBS-LRR proteins function in the plant response to pathogenic infection. Interestingly, nt-miR482b-targeted NBS-LRR proteins were identified in our AS treatment, suggesting that this previously verified regulatory mode may also be involved in alkali stress because plants can employ multiple stress response pathways to cope with environmental stress [68][69][70].
Additionally, a large number of salt or alkali stress-responsive DE miRNAs and their targets, such as nta-miR396a/U-BOX, nta-miR408/NAC domain, nta-miR482b-3p/NB-ARC domain, nta-miR6145e/leucine-rich repeat protein kinase, nta-miR6149a/F-BOX and miR-novel_93/Dof domain protein, which may provide salt and alkali stress tolerance, were predicted in our study (Table S5); however, the interactions in this regulatory network need to be further studied.

Plant Materials and Salt-Alkali Treatment
N. tabacum ecotype K326 used in this study was obtained from Henan Agricultural University. Tobacco seeds were first sterilized in 75% ethanol for 10 min, rinsed with distilled water three times and then transferred to Murashige and Skoog (MS) medium. After approximately two weeks, germinated seedlings in the MS medium were transferred to plastic pots filled with Hoagland nutrient solution and cultivated in a phytotron under normal conditions (26 • C with 14 h days, 24 • C with 8 h nights, and 70% relative humidity). Tobacco plants grown to the five-leaf stage were used in this study. For salt stress treatment, seedlings were transferred to distilled water supplemented with 100 mM NaCl (pH = 7.0). For alkali stress, seedlings were transferred to distilled water supplemented with 100 mM NaHCO 3 (pH = 8.4). After they were exposed to the stress treatment for 24 h, the samples were harvested, frozen in liquid nitrogen and then stored at −80 • C.

Construction and Sequencing of mRNA-seq and sRNA Libraries
Total RNA was extracted from tobacco leaves using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The concentration and quality of the RNA samples were assessed using a Nanodrop 2000 instrument (Thermo, Waltham, MA, USA). mRNA-seq libraries were established using an Illumina TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocol [71]. Briefly, 5 µg of mRNA was isolated from the total RNA using oligo (dT) magnetic beads (Invitrogen, Carlsbadcity, CA, USA), fragmented and reverse transcribed into cDNA. Adapters with a hairpin loop structure were ligated to cDNA molecules and amplified by PCR. The RNA-seq library was sequenced using the Illumina HiSeq X Ten platform (Illumina, San Diego, CA, USA). For small RNA (sRNA) sequencing, RNA bands of approximately 18-30 nt were isolated, linked to 5 and 3 adapters and reverse transcribed to cDNA. The final PCR products were purified and sequenced with the Illumina HiSeq 2500 sequencing platform (Illumina, San Diego, CA, USA) with 50-bp single-end reads.

Mapping of mRNA-seq Reads
To obtain clean reads, adapters and low-quality sequence reads were removed from the raw data using Trimmomatic (version 0.30) [72]. These clean reads were mapped to the tobacco genome (available online: ftp://ftp.solgenomics.net/genomes/Nicotiana_tabacum/) using HISAT2 (version 2.0.5) [73]. The gene expression levels in all nine libraries were presented as fragments per kilobase of transcript per million (FPKM) [74].

Analysis of Differential Gene Expression in Response to Salt-Alkali Stress
DEG analysis was carried out using the DESeq R package (1.10.1). Using a model based on the negative binomial distribution, this tool provides statistical routines for determining differential expression within digital gene expression data. On the basis of this method, genes with |log 2 FC|>2 and a p-value < 0.05 were classified as DEGs between two comparisons among the three experimental conditions.

GO and KEGG Pathway Enrichment Analysis
GO enrichment analysis of the DEGs was carried out using the GOseq R package with the default parameters, by which gene length bias was corrected [75]. The KEGG database is a widely used resource for understanding high-level functions and utilities of biological systems (available online: http://www.genome.jp/kegg/). The DEGs in the KEGG pathways were analyzed using KOBAS (version 2.0) [76].

Identification of Known and Novel miRNAs and Prediction of miRNA Targets
After the adapter sequences, poly-A sequences and low-quality bases were removed from the raw reads, the remaining unique sequences were mapped to the tobacco genome using Bowtie2 (version 2.0.5) [73]. The sRNA sequences that matched those of rRNA, tRNA, small nuclear RNA (snRNA) and small nucleolar RNA (snoRNA) were excluded. The remaining unique sequences were searched against miRBase 20.0 to identify known miRNAs in tobacco. The software programs miREvo [77] and mirdeep2 [78] were used to predict novel miRNAs and to explore secondary structures, and the miRNA targets were predicted using psRobot [79].

Detection of the Differential Expression of miRNAs in Response to Salt-Alkali Stress
The evaluation of the differential expression of miRNAs between two comparable groups was based on a previously established model [80]. Using |log 2 FC|>0.5 and p-value < 0.05 as criteria, we adjusted for multiple testing to classify the DE miRNAs among the three experimental conditions.

Validation of mRNA and miRNA Expression
Total RNA was extracted from each treatment using TRIzol reagent (Invitrogen, USA) as described above. cDNA was generated using a PrimeScript cDNA Synthesis Kit (Takara, Dalian, China) according to the manufacturer's instructions and then diluted 10 times for use as a qPCR template. sRNA was extracted using a Small RNA Isolation Kit (Takara, Dalian, China). RT-qPCR was performed using SYBR Green Supermix (Takara, Dalian, China) on an ABI StepOne instrument. The tobacco ef1α gene was used as an internal control for mRNA analysis, and 5.8S rRNA was used as the reference for stem-loop miRNA analysis. Three independent biological replicates were included in the qRT-PCR analyses. The relative expression level of each sample was calculated using the 2 −∆∆CT method [81]. The primers used for this RT-qPCR experiment are listed in Table S6.

Conclusions
In the present study, we comprehensively analyzed the changes in the mRNA and sRNA profiles of tobacco under NaCl and NaHCO 3 treatment. Overall, 8064 mRNAs and 33 miRNAs were differentially expressed in response to salt and alkali stress. The overlapping DEGs between the two treatments suggested that plants could utilize similar strategies, such as ion channels, AQPs and antioxidant activity, to adapt to both NaCl and NaHCO 3 stress. However, our GO and KEGG pathway enrichments revealed that many essential biological processes, including "photosynthesis" and "starch and sucrose metabolism," were specifically enriched in the NaHCO 3 treatment. In addition, our analysis also revealed many miRNAs, including four overlapping miRNA families, that might play important roles in the response of tobacco to salt and alkali stress. Several inverse regulatory processes of miRNAs with their target transcript genes were also predicted for manipulation. These findings improve our understanding of the underlying salt and alkali responses in tobacco at the molecular level and provide a valuable foundation for future uses of the miRNA-based genetic regulatory network in breeding for salinity tolerance in crops.