1 INTRODUCTION

Dunaliella salina is a halophilic model alga with an excellent ability to adapt to the environment (Chen et al., 2020). Given its inherent advantages, D. salina has been widely used in various fields, such as food processing, medical treatment, and biodiesel. For example, valuable chemicals, including carotenoids, glycerol, lipids, vitamins, minerals, and proteins, have been accumulated in different strains of D. salina (Tafreshi and Shariati, 2009; Rammuni et al., 2019). Additionally, this species can be used as an immunostimulant and recombinant vaccine for the prevention and control of aquaculture diseases (Ma et al., 2020). D. salina also has great value in biotechnology applications, including natural material preparation and biofuel production (Feng et al., 2020b). In addition to its broad application prospects, D. salina is an important type of economic microalgae that has characteristics that make it resistant to abiotic stresses (Sun et al., 2018; Moayedi et al., 2019).

Salt stress is an important abiotic stress that affects the physiology, biochemistry, and metabolism of many plants. Under extremely high salinity, the biomass production, chlorophyll, and carbohydrate contents of D. salina decrease (Tammam et al., 2011). In addition, at the biochemical level, the levels of antioxidant enzymes, antioxidant substances, and extracellular polymer (EPS) show varying degrees of change under salt stress (Mishra and Jha, 2009). Moreover, salt stress induces the activity of some metabolic enzymes in D. salina, including glycerolmetabolizing enzymes and starch-metabolizing enzymes (Panahi et al., 2019). Under salt stress, plants tend to adapt by modulating gene expression, and the regulation of gene expression mediated by microRNAs (miRNAs) is one of the main ways of adapting to this stress (Yu et al., 2016; Naik and Varadahalli, 2020). MiRNAs are endogenous non-coding RNAs that participate in post-transcriptional gene regulation by inhibiting mRNA translation or inducing mRNA degradation (Alzahrani et al., 2019). Studies have confirmed that cre-mir906-3p, cre-mir910-3p, and their targets ATP4 (d-subunit of ATP synthase) and NCR2 (NADPH: cytochrome P450 reductase) are essential for the survival of the green alga Chlamydomonas reinhardtii under abiotic stress (Gao et al., 2016). In addition, computer simulations and comparative genomics have been used to predict miRNAs across the whole genome of the brown alga Ectocarpus siliculosus, along with single-cell miRNAs and target genes in the red alga Porphyridium cruentum (Billoud et al., 2014; Barozai et al., 2018). However, research on miRNAs and salt stress has largely been focused on higher terrestrial plants (Tian et al., 2014; Wu et al., 2016; Feng et al., 2017; Yin et al., 2018; Joshi et al., 2021). Reports on lower plants, including algae, are relatively rare, leading to a lack of miRNA data in D. salina. Therefore, it is essential to discover the regulatory factors and key genes related to the salt stress response in algae, study the mechanisms of salt tolerance, and analyze the unique regulatory pathways in algae.

The emergence of high-throughput sequencing platforms, such as Roche 454, SOLiD, and HiSeq, was prompted by the rapid development of high-throughput sequencing technologies that provided the requisite tools for small RNA (sRNA) sequencing in D. salina. It is worth noting that BGISEQ-500 is a new generation high-throughput sequencing platform that is independently produced in China for small RNA and transcriptome sequencing. Related research has been reported in many plant species, including algae (Zhu et al., 2018; Lin et al., 2019). In this study, a large number of miRNAs and their corresponding target genes that were differentially expressed in response to salt stress were identified by sRNA and transcriptome sequencing analysis on a BGISEQ-500. In addition, the metabolic pathways in which miRNA-targeted genes may participate were subjected to KEGG pathway enrichment analysis, and a unique D. salina metabolic network map under salt stress was constructed. The results of this study provide insight into the salt stress regulation mechanisms of Dunaliella and related metabolic pathways and have implications for the artificial molecular improvement of algae and the in-depth study of salt stress mechanisms in lower plants.

2 MATERIAL AND METHOD 2.1 Dunaliella salina strain culture and salt stress treatment

Four strains of D. salina, namely, DS-CN1, DS-CN2, DS-C43, and DS-C48 (the first two strains were obtained from the Yellow Sea Fisheries Research Institute of Chinese Academy of Fishery Sciences and the latter two strains were obtained from the Institute of Hydrobiology, Chinese Academy of Sciences), were used to analyze algal salinity tolerance during different growth periods. Here, salinity refers to the concentration of NaCl. Dunaliella medium (DM) was used as the algal material culture medium (Supplementary Table S1), and the NaCl concentration in the culture solution was either 0.05 mol/L (basic physiological salinity for algae cultivation indoors, control group: CO) or 3.42 mol/L (high salinity for algae cultivation indoors, high salt stress group: SS) based on our previous research results (Gao et al., 2021). The light intensity for indoor algae cultivation was 10 000 lx, the temperature was 25±2 ℃, and the photoperiod was 12 h꞉12 h.

2.2 Screening for algae strains with high salt tolerance

The following index values were used to screen comprehensively for algae strains with high salt tolerance. The experimental flow chart is shown in Supplementary Fig.S1a.

2.2.1 Morphological observation and density determination of algae cells

Fifty milliliters each of algal culture was drawn from the CO and SS groups at different growth stages, and then the corresponding algal cells were centrifuged at 5 000 r/min for 3 min and immobilized using a small amount of Lugol's solution (diluted iodine and potassium iodide in 100-mL ultrapure water were prepared at a mass ratio of 1꞉2). An electron microscope (BX-51, Olympus Corporation, Tokyo, Japan) was used to observe the morphology of the algal cells at different growth stages (50× magnification). One milliliter of algal culture was taken every three days, and a DR6000 UV-visible spectrophotometer (HACH, Colorado, USA) was used to measure its absorbance at 686 nm based on our previous research (Xu et al., 2022). A cell counter (Shanghai Optical Instrument Factory, Shanghai, China) was used to measure the algal cell density at different growth stages, and each count was repeated five times.

2.2.2 Determination of chlorophyll and carotenoid content

Algal culture (50 mL) was drawn from the CO and SS groups and centrifuged at 5 000 r/min for 8 min each. After the supernatant was removed, an equal volume of 95% ethanol was added for pigment extraction. The reaction solution was placed in a thermostat at 4 ℃ for 24 h and then centrifuged at 10 000 r/min for 10 min. Next, 3 mL of the supernatant was removed, and the absorbance values at 665 nm, 649 nm, and 470 nm were measured with a DR6000 UV-visible spectrophotometer. Each measurement was repeated three times. The contents of chlorophyll a (Chl a), chlorophyll b (Chl b), and carotenoids were calculated using the following formulas Chlorophyll a, chlorophyll b, and carotenoid have the maximum absorption peaks within these wavelength ranges. Their contents in the algae could be calculated through the peak wavelength differences. (Mera et al., 2016; Liu et al., 2019; Yang and Hu, 2020):

Chl a (mg/L)=13.95×OD₆₆₅–668×OD₆₄₉,

Chl b (mg/L)=24.96×OD₆₄₉–7.32×OD₆₆₅,

Carotenoid (mg/L)=[1000×OD₄₇₀–2.05×Chl a–114.8×Chl b]/245.

2.2.3 Measurement of osmotic potential and plasma membrane relative conductivity in algal cells

In accordance with our previous research methods (Gao et al., 2022), a Vapro 5600 Dew Point Osmometer (Wescor, Missouri, USA) was used to measure the osmotic potential of algal cells from the CO and SS groups, and each test was repeated three times. A DDs-11 conductivity meter (INESA Scientific Instrument Co., Ltd., Shanghai, China) was used to measure the background conductivity of the algal and culture solutions at room temperature (25 ℃) and high temperature (95 ℃ water bath for 15 min), and each test was repeated three times. The relative conductivity of the algal cells was calculated using the following formula (Zhang et al., 1989):

Relative conductivity (%)=(conductivity of the algae solution at room temperature–conductivity of the culture solution at room temperature)/(conductivity of the algae solution at high temperature–conductivity of the culture solution at high temperature)×100%.

2.2.4 Determination of total lipid and polysaccharide contents

To determine the total lipid content of a single cell of each algae strain, 50 mL of a mixed solvent consisting of 8꞉4꞉3 chloroform꞉methanol꞉water by volume was used, and 50 mL of clean algal cell solution was washed with ddH₂O and was then subjected to full fat extraction (which is referred to as the Folch method) (Folch et al., 1957; Ametaj et al., 2003). Each sample was subjected to three replicate measurements. The polysaccharide content was determined as follows: first, 50 mL of algal cell solution was centrifuged at 6 000 r/min for 5 min, and the supernatant was discarded. Next, the algae were incubated at -80 ℃ for 30 min and placed in a YZG-600 vacuum dryer (Yutong drying equipment Co. Ltd., Changzhou, China) for 24 h, and the weight of the algae powder was then measured. Additionally, 30 mL of ultrapure water was added to the algae powder, and the algae solution was concentrated using a RE-2000A rotary evaporator (YaRong Biochemistry Instrument Factory, Shanghai, China). The supernatant was collected after centrifugation at 10 000 r/min for 10 min. After 1 mL of supernatant was removed, 5 mL of concentrated sulfuric acid and 1 mL of 6% phenol solution (phenol-sulfuric acid method) were added (Rasouli et al., 2014; Zavřel et al., 2018). All the reaction solutions were subsequently boiled in a water bath at 95 ℃ for 15 min and then cooled to room temperature. Lastly, the absorbance was measured at 480 nm. The polysaccharide concentration of the algal cells was calculated according to the standard curve formula shown in Supplementary Fig.S1b. Each test was repeated three times.

2.3 Small RNA sequencing

Three algal cultures (200 mL for each) of the screened unique halophilic Dunaliella strain in the log phase above were taken from the SS and CO groups: SS1, SS2, SS3, CO1, CO2, and CO3. Six samples were centrifuged, and the corresponding algal cells were collected (approximately 80 μg for each). These cells were flash-frozen with liquid nitrogen and ground for 20 s. A TRIzol^® kit (Thermo Fisher Scientific, Massachusetts, USA) was used for total RNA extraction, and the RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., California, USA). High-quality RNA samples (OD260/280 ratio: 1.8–2.0, RNA integrity number (RIN) value >7.5) were used to construct sRNA libraries (MGIEasy Small RNA library kit, MGI Technology Co., Ltd., Shenzhen, China). High-throughput sequencing was performed using a BGISEQ-500 sequencing platform (BGI Co., Ltd., Shenzhen, China). After the contaminated fragments were removed (the same base number as 5ʹ or 3ʹ end joint >5, non-insert fragments, N base of fragment ratio >5%, fragments smaller than 18 nt, and poly A fragments) by running SOAPnuke (https://github.com/BGI-flexlab/SOAPnuke) (Chen et al., 2018) and low-quality fragments (fragment base quality score < 20) were checked by performing the sequencing data quality control of FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) (Lokhande, 2023); sequenced fragments with the same sequence information were combined and recorded as unique fragments. In addition, the total number of sequenced and unique fragments in the CO and SS libraries were counted separately, and the length of high-quality fragments was calculated by running Trimmomatic (http://www.usadellab.org/cms/index.php?page=trimmomatic) (Sewe et al., 2022). The raw sequencing data from the sRNA libraries constructed from all six algal samples were uploaded to the NCBI Sequence Read Archive (SRA) database in 2020 (accession numbers: SRR8389145, SRR8389147, SRR8389148, SRR8389149, SRR8389150, and SRR8389153) and the CNCB Genome Sequence Archive (GSA) database in 2022 (accession numbers: CRA006740, CRA006741, CRA006742, CRA006743, CRA006744, and CRA006745).

2.4 Classification and annotation of sRNAs

High-quality sRNA fragments were mapped to the D. salina reference genome (Dunsal1 v.2, INSDC: NSFN00000000.2) using Bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) (Langdon, 2015) to analyze sRNA expression and distribution with the R language. sRNA fragments were annotated in the miRBase (release v22.0, http://www.mirbase.org), Rfam (release v14.0) (Kalvari et al., 2018), siRNA (http://sirna.cgb.ki.se), piRNA (http://www.regulatoryrna.org/database/piRNA), and snoRNA (https://www-snorna.biotoul.fr) databases. In the order of miRNA>piRNA>snoRNA>Rfam>other sRNAs, each unique small RNA had a unique annotation. Moreover, high-quality sequencing fragments with annotated information were classified and counted, and conserved D. salina miRNA sequences were obtained.

2.5 Prediction of novel miRNAs

Novel miRNA prediction was performed on highquality sRNA fragments that were mapped to the reference genome by running SOAPaligner (http://soap.genomics.org.cn/soapaligner.html) (Garcia-Seco et al., 2015) with no annotations. Sequences that met the screening criteria of Axtell and Meyers (2018) were predicted to be novel miRNAs. MiRA software (https://github.com/mhuttner/miRA) (Evers et al., 2015) and Origin software (https://www.originlab.com) were used to construct a miRNA precursor (premiRNA) secondary structure and count the length and quantity, respectively, of mature miRNA sequences.

2.6 Screening of miRNAs in response to salt stress

Based on the expression of miRNAs in the CO libraries (CO1, CO2, and CO3), the fold-change value of miRNAs in the SS (SS1, SS2, and SS3) libraries [fold-change=log2(SS/CO)] was determined using DESeq2 (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html) (Love et al., 2014), and the significance of miRNA differential expression was detected using the corrected P value (FDR) (P value correction was performed using R). Then, significantly differentially expressed miRNAs (DEMs) were identified with the threshold of log2|fold-change| >4.5 and -log10FDR >50. In addition, common DEMs were identified among the six libraries (SS and CO). In addition, a hierarchical clustering of DEMs (based on their relative expression levels followed by FPKM value calculation) was performed using MeV software (http://mev.tm4.org) (Saeed et al., 2003).

2.7 Prediction and functional annotation of miRNA target genes

Based on TargetFinder (http://targetfinder.org), TAPIR (http://bioinformatics.psb.ugent.be/webtools/tapir), and psRobot (http://omicslab.genetics.ac.cn/psRobot), the miRNA sequences responsive to salt stress and previously published transcriptomic data (Gao et al., 2021) (NCBI SRA: SRR8393721, SRR8393728, SRR8393729, SRR8393722, SRR8393723, and SRR8393725) were used for BlastN mapping to screen the miRNA target genes. In addition, Gene Ontology (GO, http://geneontology.org) and Web Gene Ontology Annotation Plot (WEGO, https://wego.genomics.cn) were used to annotate and classify the predicted target genes. The enrichment of the classified GO terms was assessed using GOEAST software (http://omicslab.genetics.ac.cn/GOEAST) (E-value cut-off ≤10^-5). Using the entire transcriptome data as the background, the prediction and enrichment of Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/kegg/pathway.html) pathways were performed on the target genes that were assigned KEGG ontology (KO) terms through the KEGG Automatic Annotation Server (KAAS, https://www.genome.jp/tools/kaas).

2.8 Quantitative real-time PCR of key miRNAs and target genes

The key miRNAs and corresponding target genes involved in the enrichment analysis were validated by quantitative real-time PCR (qPCR). RNA samples from the SS and CO groups were reverse transcribed according to the instructions of the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Massachusetts, USA). Specific primers based on miRNAs and corresponding target sequences (targets) were designed using Primer Premier 5.0, and qPCR amplification was performed (NCode^TM EXPRESS SYBR^® GreenER^TM miRNA qRT-PCR Kit, TaKaRa, Japan). SPSS software and the Student's t-test were used to calculate the 2^-ΔΔCt value to obtain the relative expression levels of miRNAs and their target genes. The D. salina 18S rRNA gene was used as an internal control. All samples were subjected to three biological and technical replicates.

2.9 Construction of a metabolic network involving key miRNAs

The significantly enriched KEGG pathways (enrichment factor value >0.10), the corresponding target genes and miRNAs, and their functional annotations were screened. Cytoscape network software (https://cytoscape.org) was used to construct a metabolic network diagram of the D. salina response to salt stress with the key involved miRNAs.

3 RESULT 3.1 Screening for a unique halophilic D. salina strain

In both the SS and CO groups, the algal cells were both ellipsoidal. Most of the algal cells in the SS group showed aggregate growth, and the algal cell solution was dark green. By contrast, the algal cells of the CO group grew more dispersed and the algal cell solution was light green, which was investigated using all four Dunaliella strains with DS-CN1 as a typical example (Fig. 1a). Under both physiological salinity and high salinity stress, the cell density of the algae strain DS-CN1 was the highest (Fig. 1b). In this strain, the average cell densities of the SS and CO groups at the 30-day maturity stage were (3.89±0.35)×10⁶ cells/mL and (6.46±0.54)×10⁶ cells/mL, respectively. The average absolute values of the single-cell osmotic potential for the 30-day maturity algae strain DS-CN1 in the SS and CO groups were lower than those of other algae strains at (3.05±0.62)×10^-6 MPa/cell and (1.65± 0.35)×10^-6 MPa/cell, respectively (Fig. 1c). Similarly, the average percentage of the plasma membrane relative conductivity for the DS-CN1 strain was lower than that of other algal strains in both the SS and CO groups, with average percentages of 2.36%± 0.15% and 1.42%±0.08%, respectively (Fig. 1c). Moreover, compared with other strains, the chlorophyll a, chlorophyll b, and carotenoid contents were the highest in the DS-CN1 strain, indicating that the levels of photosynthesis and carotenoid metabolism remained high under high salt stress (Fig. 1d). In addition, the total lipids and polysaccharides in single cells was highest in the 30-day maturity DS-CN1 strain, revealing that the lipid and sugar metabolism levels were greatly enhanced under high salt stress (Fig. 1e). Based on the above results, DS-CN1 was selected as the experimental algae strain for the following further study.

3.2 Variation in sRNA expression in response to salt stress

The SS and CO libraries from DS-CN1 generated 36.93-M and 34.98-M fragments, respectively. After the contaminated and low-quality fragments were filtered out, 33.83-M (91.62%) and 30.70-M (87.75%) fragments were retained in the SS and CO libraries, respectively (Supplementary Table S2a). Compared with the CO library, a higher percentage of sequenced fragments (on average) in the SS library were mapped to the reference genome (32.78% vs. 23.82%, Supplementary Table S2b). The unique total sRNA fragments in the SS and CO libraries were 21.95-M and 28.92-M, respectively (Fig. 2a). Additionally, the numbers of unique sRNA fragments were 21.07-M and 27.63-M in the SS and CO libraries, respectively (Fig. 2b), and the average number of RNA fragments per unique sequence was 1.1. Overall, 2.63% and 2.54% of the total number of sRNAs and unique sRNAs, respectively, were shared between the SS and CO libraries. These results revealed that the expression levels of some sRNAs were inhibited by high salt stress. In addition, although the length of sRNA fragments in the SS and CO libraries ranged from 10 nt to 44 nt (Fig. 2c), the lengths of most sRNA fragments were concentrated at 21 nt and 22 nt (28.05% and 25.13% in the SS libraries and 28.62% and 27.44% in the CO libraries, respectively). Thus, the number of sRNAs with lengths of 21 nt and 22 nt was higher in the SS group than in the CO group, suggesting that more sRNAs within this length range were stimulated by high salt stress in DS-CN1. The flow chart of bioinformatic analysis is shown in Supplementary Fig.S2.

3.3 sRNA classification and annotation

The percentages of all the annotated sRNAs in the SS and CO libraries were 5.94% and 6.18%, respectively, while the percentages of all the annotated unique sRNAs were 5.72% and 6.06%, respectively. In addition, the number and proportion of sRNA sites that were annotated as siRNAs in the CO library were lower than those in the SS library (Table 1), suggesting that more siRNAs were stimulated in response to salt stress. Information regarding the siRNA sequences identified in DS-CN1 is shown in Supplementary Table S3, and the base preference results for each position are shown in Supplementary Fig.S3a. Although no conserved miRNAs in DS-CN1 were annotated in the database, a total of 44 novel miRNAs were identified in this species through a predictive analysis of unannotated sRNA sequences (Table 2). These novel miRNAs were located across different positions in the reference genome, and most were located on the sense strand (Table 2). In addition, the base ratio of each position differed, and the length was 19–24 nt. Among the miRNAs, the percentages of the base C were relatively high (Fig. 3a), and the 24-nt miRNAs had the most fragments in the SS and CO libraries (the average fragments accounted for 30.21% of all fragments) (Fig. 3b). The precursors of novel miRNAs had stable secondary structures (MFE values less than -20 kcal/mol). Based on the stem loop characteristics of their precursors, these miRNAs can be divided into three categories: (1) a typical stem-loop structure with a ring-shaped branch or bubble-like structure; (2) a typical stem-loop structure without branches; and (3) an atypical stem loop with more than two ring branches or bubble structures (Axtell, 2013). The representative secondary structures of these three categories of novel miRNA precursors are shown in Fig. 3c, and the structures of the remaining novel miRNA precursors are shown in Supplementary Fig.S3b.

3.4 Identification of salt stress-responsive novel miRNAs

By comparing the expression levels of novel miRNAs in the SS and CO libraries (Supplementary Table S4), we found that 14 novel miRNAs were downregulated and 16 novel miRNAs were upregulated in response to salt stress. In particular, dsa-mir40 showed obvious downregulation (fold change < -24.5, P < 1.00E^-50), while dsa-mir3, dsa-mir16, dsa-mir17, and dsa-mir26 were significantly upregulated (fold change >24.5, P < 1.00E^-50). There were an additional 14 novel miRNAs that showed minimal differences in expression (|fold change| < 2.0, P < 0.05) between the SS and CO libraries (Fig. 4a). Through a Venn diagram analysis of novel miRNAs in the six sRNA sequencing libraries, a total of 21 common differentially expressed miRNAs (DEMs) were identified (Fig. 4b). Furthermore, a cluster analysis was performed to determine the expression levels of these DEMs. The results show that these DEMs could be clustered into two large groups of those that were upregulated and downregulated according to their expression levels in the SS and CO libraries of DS-CN1 (Fig. 4c).

3.5 Prediction and functional annotation of novel miRNA-regulated genes related to the salt stress response

After the identified novel miRNAs in DS-CN1 mapped to previously published D. salina transcriptome sequencing data (Gao et al., 2021). 319 potential miRNA target genes were jointly analyzed and predicted using TAPIR, TargetFinder, and psRobot (Supplementary Table S5a). A total of 104 overlapping genes were identified (Fig. 5a). In addition, based on the unique gene function annotation, 253 target genes were subjected to GO functional annotation and classification (Supplementary Table S5b). This analysis included 97 genes annotated as biological processes, 99 genes annotated as cellular components, and 57 genes annotated as molecular functions (Fig. 5b). Furthermore, KEGG pathway enrichment analysis was performed on the 104 target genes (Supplementary Table S5c), and a total of 57 KEGG pathways were predicted. In particular, the six most significantly enriched pathways (enrichment factor >0.01) from the top 20 pathways included amino acid metabolism, secondary metabolic biosynthesis, energy metabolism, carbon anabolism, overall metabolic pathways, translation and transcription (Fig. 5c). Moreover, through in the detection of the 6 key KEGG pathways, five novel miRNAs, namely, dsa-mir3, dsa-mir16, dsa-mir17, dsa-mir26, and dsa-mir40, were found to be significantly related to salt tolerance metabolism. These miRNAs were found to target Unigene51665, Unigene32647; CL371. contig2, and BM448588.1; CX160958.1, CX160930.1; Unigene43639, Unigene1468; and Unigene13268, Unigene13536, respectively (Table 3). The functions of these target genes with putative encoding proteins were annotated to 52 GO terms and 31 KO terms. These protein functions and corresponding miRNA-target gene pairs are as follows: metabolism of aspartate kinase (dsa-mir3-Unigene5 1665), genetic information processing of serine/arginine repetitive matrix protein 1 (dsa-mir3-Unigene32647), genetic information processing of 60S ribosomal protein L27a-3 (dsa-mir16-CL371. Contig2), genetic information processing of large subunit ribosomal protein L27Ae (dsa-mir16-BM448588.1), metabolism of photosystem Ⅱ 10-kDa protein (dsa-mir17-CX160958.1), metabolism of hypothetical protein GPECTOR_10g759 (dsa-mir17-CX160930.1), metabolism of ribulose bisphosphate carboxylase small chain 1 (dsa-mir26-Unigene43639), metabolism of chloroplast ribulose-1, 5-bisphosphate carboxylase, oxygenase small subunit (dsa-mir26-Unigene1468), metabolism of DNA-directed RNA polymerase subunit beta (dsa-mir40-Unigene13268), and genetic information processing of small subunit ribosomal protein S11 (dsa-mir40-Unigene13536).

3.6 qPCR analysis of novel miRNAs and target genes in response to salt stress

The relative expression of five novel miRNAs and their corresponding ten target genes, which were significantly related to salt stress in the SS and CO groups, were analyzed by qPCR (the detection gels of novel miRNAs and their target genes are shown in Supplementary Fig.S4a–d. The qPCR primer sequences are listed in Supplementary Table S6). Compared with the CO group, dsa-mir3, dsa-mir16, dsa-mir17, and dsa-mir26 were significantly upregulated (P < 0.01) in the SS group, while dsa-mir40 was significantly downregulated (P < 0.01) in the SS group (Fig. 6a). By contrast, compared with the expression patterns of novel miRNAs in the two groups, their corresponding target genes showed opposing differential expression patterns (Fig. 6b); that is, in the SS group, the target genes of dsa-mir3, dsa-mir16, dsa-mir17, and dsa-mir26 were significantly downregulated (P < 0.01), while the target genes of dsa-mir40 were significantly upregulated (P < 0.01).

3.7 Analysis of key metabolic networks mediated by novel miRNAs in D. salina in response to salt stress

Following the integrated analysis of five key novel miRNAs, target genes, coding products, and main KEGG pathways, the key metabolic network of D. salina in response to salt stress was constructed (Fig. 7). In the SS group, dsa-mir3 was upregulated, and the synthesis of the coding products of Unigene51665 and Unigene32647 were inhibited, which affected translation, amino acid metabolism, secondary metabolism, and global and overview maps. Dsa-mir16 was upregulated in D. salina, and the synthesis of large subunit ribosomal proteins L27Ae and 60S ribosomal protein L27a-3, encoded by CL371 and BM448588.1, respectively, was inhibited, thereby regulating the translation process. dsa-mir17 was upregulated, which inhibited the synthesis of CX160958.1 and CX160930.1, and in turn affected energy metabolism. dsa-mir26 expression was upregulated and inhibited the synthesis of the Unigene43639 and Unigene1468 gene products, which affected the carbohydrate metabolism process in this alga.

By contrast, in the SS group, dsa-mir40 was downregulated and stimulated the synthesis of RNA polymerase β subunit (rpoB) and the small subunit ribosomal protein S11 via the target genes Unigene13268 and Unigene13536, which affected the transcription process and global and overview maps. In short, as part of the overall metabolic network, there was a close interaction among carbohydrate metabolism, energy metabolism, transcription and translation, amino acid metabolism, and secondary metabolism in D. salina under high salt stress.

4 DISCUSSION

The aim of this work was to study the molecular mechanisms underlying metabolism in mature D. salina in response to salt stress through sRNA sequencing and bioinformatic analysis. In this study, we identified five key novel and differentially expressed dsa-miRNAs and ten miRNA target genes involved in the metabolic regulation pathways of salt tolerance in D. salina. In addition, combined with the construction of a metabolic network for D. salina in response to salt stress and recent literature reports on genes related to salt stress, four miRNA-target gene pairs were further emphasized: dsa-mir26-Unigene43639, dsa-mir26-Unigene1468, dsa-mir17-CX160958.1, and dsa-mir40-Unigene13536.

In plants, the salt stress response is highly complex and involves signal transduction pathways, the activation, and synthesis of stress-regulated genes, and multiple regulatory proteins (Yang and Guo, 2018). Here, we identified pathways that may be involved in the salt stress response, such as genetic information processing (ko03010, ko02900, and ko03010) and metabolism (ko00195). These pathways are involved in metabolic pathways, amino acid metabolism, secondary metabolite synthesis, and translation. These results showed that salt stress triggers a complex response in D. salina.

Unigene43639 and Unigene1468 were identified in the constructed metabolic network and encode ribulose bisphosphate carboxylase/oxygenase (RuBisCO), which catalyzes the first major carbon fixation reaction in photosynthesis and converts free carbon dioxide in the atmosphere into organismal energy storage molecules (Feng et al., 2020a). Therefore, based on the metabolic network diagram, dsa-mir26-Unigene43639 and dsa-mir26-Unigene1468 were thought to participate in sugar metabolism and photosynthesis mediated by RuBisCO in D. salina under salt stress. However, the downregulated expression of the target genes Unigene43639 and Unigene1468 in DS-CN1 suggested that the salt tolerance of DS-CN1 might be triggered by other mechanisms. The increased polysaccharide content in DS-CN1 may be regulated by other genes. It has been confirmed in land plants such as rice and alfalfa that RuBisCO participates in salt tolerance in plants by affecting photosynthesis (Jin et al., 2010; Frukh et al., 2020). It is worth noting that the chlorophyll a and b contents in the salt-tolerant strain DS-CN1 were significantly higher than those in the other salt-intolerant strains under salt stress. We speculate that this result was related to the significant upregulation of dsa-mir26 and the downregulation of Unigene43639 and Unigene1468 expression in DS-CN1, which reduced the damaging effects of salt stress on photosynthesis.

The photosystem Ⅱ 10-kDa protein is a transcriptional regulatory node identified by the microarray analysis of salt stress tolerance in rice and was significantly downregulated in salt-tolerant rice varieties (Pandit et al., 2011). However, there have been few studies on salt tolerance in other plants. Similarly, the expression level of CX160958.1, which encodes the photosystem Ⅱ 10-kDa protein, was significantly downregulated in the D. salina DS-CN1 strain. Therefore, dsa-mir17-CX160958.1 may serve as a new regulatory node that participates in the salt-tolerant metabolic pathways of D. salina.

In the metabolic network, the dsa-mir16-CL371, Contig2, dsa-mir16-BM448588.1, and dsa-mir40-Unigene13536 pairs were involved in the genetic information processing of ribosomal proteins. In addition, many studies have confirmed that some ribosomal protein genes are induced in response to salt stress in some lower and higher organisms (Duché et al., 2002; Kanesaki et al., 2002; Omidbakhshfard et al., 2012). Under salt stress, an analysis of the expression level of five novel dsa-miRNAs revealed that only dsa-mir40 expression decreased significantly, resulting in a significant increase in the expression level of its target gene Unigene13536, which encodes the small subunit ribosomal protein S11. However, changes in the expression of dsa-mir16-CL371. Contig2 and dsa-mir16-BM448588.1 were the opposite of that of dsa-mir40-Unigene13536. Among the above three miRNA-target gene pairs, only dsa-mir40-Unigene13536 exhibited a change in expression that was consistent with the tendency of ribosomal proteins to be induced under salt stress. We speculated that dsa-mir40-Unigene13536 might be involved as a regulatory factor in the salt stress tolerance of D. salina. However, supporting biochemical evidence is still lacking, and additional experimental evidence, such as detecting the expression of ribosomal proteins in-vivo and in-vitro, regarding the role of dsa-mir40-Unigene13536 in response to salt stress is needed.

The analysis of the metabolic network diagram demonstrates that we lack qualitative or quantitative physiological and biochemical data on the SS and CO groups, such as data on ribosomal proteins. Combined with the miRNA-target gene relationship pairs we obtained, a deeper analysis of the D. salina metabolic pathway is needed. In addition, the screening and identification of upstream regulatory factors, such as long non-coding RNAs (lncRNAs) and an analysis of lncRNA-miRNA-mRNA interactions in response to salt stress can further improve our understanding of the salt tolerance mechanism of D. salina and potentially expand to applications in other plants.

5 CONCLUSION

In this study, by using small RNA sequencing, transcriptome sequencing, and bioinformatic analysis methods, we constructed a metabolic network diagram of D. salina (DS-CN1) under salt stress to investigate the molecular pathways associated with the salt stress response. The miRNA-target gene pairs, dsa-mir26-Unigene436391, dsa-mir26-Unigene1468, dsa-mir17-CX160958.1, and dsa-mir40-Unigene13536 may be closely associated with salt tolerance in D. salina. The genes involved in these miRNA-mRNA pairs may provide a basis for in-depth research on the salt tolerance mechanisms of additional species.

6 DATA AVAILABILITY STATEMENT

The raw sRNA-seq datasets generated in this study are available from the NCBI SRA database under the accession numbers SRR8389145, SRR8389147, SRR8389148, SRR8389149, SRR8389150, and SRR8389153 and the CNCB GSA database under the accession numbers CRA006740, CRA006741, CRA006742, CRA006743, CRA006744, and CRA006745. The authors declare that all the other data supporting the findings of this study are available within the article and its supplementary files.

7 AUTHOR CONTRIBUTION

Fan GAO performed the experiments, analyzed the data, and co-wrote the manuscript; Fangru NAN analyzed the data and co-wrote the manuscript; Jia FENG, Junping LÜ, Qi LIU, and Xudong LIU analyzed the data; and Shulian XIE supervised the project and revised the manuscript. All the authors read and approved the final manuscript.

All the samples were obtained from the wellknown algae culture and collection, and the sampling procedures in this experiment were performed in accordance with the guidelines for the care and use of laboratory plants of Shanxi University and were approved by the Wildlife Care and Use Committee of Taiyuan Forestry Bureau, Shanxi Province, China.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S6 and Figs.S1–S4) is available in the online version of this article at https://doi.org/10.1007/s00343-022-2130-1.