Effects of NaCl stress on growth indexes
G. divaricata shoots had a high rooting percentage (98.89%) within 7 d without NaCl, and showed 28.89–83.33% rooting percentage at low NaCl concentrations (50–100 mM). However, at high concentrations (≥ 200 mM), there was basically no rooting within 7 d, when cultured to 21 d, 100% rooting at 0–50 mM NaCl, and no rooting at a high concentration (400 mM NaCl), the abaxial surface of leaves were dark, while leaves themselves were curled and brittle (Fig. 1). Generally, as NaCl concentration increased, rooting percentage decreased (Table 1).
When cultured at 0–300 mM NaCl, there was no plantlet death, but when the NaCl concentration was 400 mM, a small number of plantlets died. The average number of roots and root length decreased as NaCl concentration increased, but average root length increased when the NaCl concentration was 50 mM. Plantlet height decreased from 12.27 cm to 7.89 cm with increasing NaCl concentration. Above-ground biomass per plantlet reached a maximum (245.89 mg) at 200 mM NaCl, and plantlet water content was highest (95.30%) at 50 mM NaCl but gradually decreased to 92.4% at 400 mM NaCl (Table 1).
Antioxidant enzyme activity, and MDA, proline, and soluble sugar content
CAT activity of 50–400 mM NaCl-treated G. divaricata leaves was significantly higher than in control leaves, and highest at 50 mM (Fig. 2a). POD activity in the 200-400 mM NaCl treatments was significantly higher than in control leaves, and highest at 400 mM (Fig. 2b). SOD activity in all NaCl treatments was significantly higher than in leaves of the control group, and highest at 300 mM (Fig. 2c). MDA content in all NaCl treatments was significantly higher than in the control group, and highest at 200 mM (Fig. 2d). Except for 50 mM NaCl, proline content in all other NaCl treatments was significantly higher than in the control group, and highest at 200 mM (Fig. 2e). Except for 50 mM, soluble sugar content in all NaCl treatments was significantly lower than in the control group, and lowest at 300 mM (Fig 2f).
Na+, K+ and Mg2+ content
Na+ content accumulated significantly in leaves relative to control leaves, reaching a maximum at 200 mM (Fig. 3a). Conversely, K+ decreased significantly in leaves relative to control leaves, reaching a minimum at 50 and 200 mM (Fig. 3b). Mg2+ content decreased significantly at 100, 300 and 400 mM relative to control leaves (Fig. 3c). Na+ / K+ content was significantly higher in all NaCl treatments relative to the control, and maximum at 200 mM (Fig. 3d).
Correlation analysis of growth and physiological indicators under NaCl stress
The correlation analysis of growth and physiological indicators of G. divaricata leaves in response to different NaCl treatments showed that NaCl concentration was positively correlated with POD activity, as well as Pro, Na+, and Na+/K+ content, but negatively correlated with soluble sugar and Mg2+ content, as well as root number, root length, plant height and water content (Table 2).
De novo assembly and sequencing quality analysis
In 0, 50 and 200 mM NaCl stress treatments, 43405450–57982988 clean reads were obtained after filtering the original sequencing data, Q30 exceeded 92.58%, and GC content lay in the 42.91–43.53% range (Table 3). The length distribution of the assembled transcripts and unigenes are shown in Fig. 4.
Total number of annotated unigenes
The total number of unigenes functionally annotated by six databases, NR, GO, KEGG, clusters of orthologous groups (COG), swiss prot protein database (Swiss Prot) and the protein family database (Pfam), was 57297, and 23395 (40.83%), 17564 (30.65%), 11060 (19.30%), 25070 (43.75%), 22,062 (38.50%), and 20982 (36.62%), respectively, with 27588 (48.15%) unigenes annotated to at least one database (Supplementary Table 1).
Among the genes annotated in the NR database, 76.71% (66.45% + 10.26%) had an E value < 1e-20 (Fig. 5a), and 69.40% (49.48% + 19.92%) had a genetic similarity > 60% (Fig. 5b). 52.32% of the genes were annotated in known species, the most similar species being Selaginella moellendorffii (16.53%), followed by Marchantia polymorpha (9.64%) (Fig. 5c).
In the GO database, a total of 73360 unigenes were annotated into three categories: molecular functions, biological processes and cellular components, accounting for 29.82%, 34.06% and 36.12%, respectively. In the molecular function category, they were mainly associated with binding and catalytic activity, among biological processes, they were mainly associated with cellular and metabolic processes, and in the cell components category, they were mainly associated with cell and membrane parts, and organelles (Fig. 6a); 73.54% of unigenes were associated with seven categories of GO terms. A total of 7729 unigenes were annotated in in five branches of the KEGG database: metabolism, genetic information processing, environmental information processing, cellular processes and organism systems, with metabolism accounting for 54.03%; most genes were associated with carbohydrate metabolism, translation, folding, sorting and degradation, amino acid metabolism, and transport and catabolism (Fig. 6b). A total of 13233 unigenes were associated with unknown functions in the COG database while 13399 unigenes were associated with 22 known categories, with the largest number of annotations related to posttranslational modification, protein turnover, chaperones (1790), followed by replication, recombination and repair (1753), signal transduction mechanisms (1648), transcription (1418), and others (Fig. 6c).
Analysis of DEGs
NaCl stress induced the differential expression of genes. As NaCl concentration increased, the number of DEGs in tissue-cultured G. divaricata leaves also increased. In the L50 vs CK, L200 vs CK, and L200 vs L50 comparisons, 177 and 128, 1396 and 1158, and 1024 and 940 DEGs were up- and downregulated, respectively (Fig. 7a). The number of unique DEGs in common to all three groups was 80, while 48, 1243 and 772 DEGs were unique to L50 vs CK, L200 vs CK, and L200 vs L50 comparisons, respectively (Fig. 7b).
In the L50 vs CK comparison, mevalonate kinase activity was highest, followed by response to hydrogen peroxide while most DEGs were related to oxidoreductase activity, followed by lipid metabolic processes (Fig. 8a). In the L200 vs CK comparison, highest expression was related to rhythmic processes, followed by photomorphogenesis, while most DEGs were related to DNA binding, followed by transporter activity (Fig. 8b). In the L200 vs L50 comparison, highest expression was related to the ethylene-activated signaling pathway followed by rhythmic processes, while most DEGs were related to DNA binding, followed by transcription regulator activity (Fig. 8c).
KEGG enrichment analysis showed that most DEGs in the L50 vs CK comparison were related to linoleic acid metabolism followed by zeatin biosynthesis (Fig. 9a). In the L200 vs CK and L200 vs L50 comparisons, most DEGs were related to plant hormone signal transduction, MAPK signaling pathway-plant, and starch and sucrose metabolism (Fig. 9b, 9c).
Genetic analysis related to plant hormone signaling, MAPK signaling and starch, sucrose metabolic and ion transport pathways
PYL, PP2Cs and SnRK2s are important components of the abscisic acid pathway. PYL expression was downregulated under salt stress, resulting in the activation of PP2C expression, which inhibited SnRK2s expression. ARF in the auxin pathway and some IAAs and SAURs, CRE1 and A-ARR in the cytokinin pathway, GID1 in the gibberellin pathway and TGA in the salicylic acid (SA) pathway were significantly upregulated (Fig. 10). Seven genes, including TRINITY_DN6033_c0_g1, TRINITY_DN27156_c0_g1 and TRINITY_DN4509_c0_g1 in the MAPK signaling pathway were highly expressed in the 200 mM NaCl treatment while eight genes, including TRINITY_DN2820_c0_g2, TRINITY_DN2815_c0_g2 and TRINITY_DN29407_c0_ g2, were highly expressed only in the control group or in the 50 mM NaCl-treated group (Fig. 11). Under salt stress, the expression of sucrose invertase (INV, sacA), sucrose synthase (SUS), and TRINITY_DN15485_c0_g2 (E2.7.1.4, scrK) was upregulated (Fig. 12).
KEGG enrichment showed that 10 genes associated with the plant hormone signaling pathway, like TRINITY_DN2820_c0_g2, TRINITY_DN2815_c0_g2, and TRINITY_DN29407_c0_g2, were only highly expressed in the control group or in the 50 mM NaCl-treated group, while 15 genes, like TRINITY_DN6033_c0_g1, were highly expressed in the 200 mM NaCl-treated group (Fig. 13a). Nine genes (TRINITY_DN29407_c0_g2, TRINITY_DN2815_c0_g2, TRINITY_DN9614_c0_g1, TRINITY_DN8933_c0_g1, TRINITY_DN4509_c0_g1, TRINITY_DN6033_c0_g1, TRINITY_DN2820_c0_g2, TRINITY_DN14265_c0_g1, TRINITY_DN12556_c0_g1) were co-mediated two pathways, plant hormone signaling and the MAPK pathway (Supplementary Tables 1 and 2). The expression of four PP2Cs in the abscisic acid pathway were up-regulated in response to salt stress, as was one ERF1 in the ethylene pathway, and WRKY22, related to ethylene and hydrogen peroxide responses (Fig. 13b). KEGG enrichment analysis showed that six genes in the starch and sucrose metabolic pathway, including TRINITY_DN4486_c1_g1, TRINITY_DN6333_c0_g1, and TRINITY_DN32855_c1_g2, were highly expressed in the 200 mM NaCl treatment while five genes, including TRINITY_DN209_c0_g1, TRINITY_DN4277_c0_g1, and TRINITY_DN6608_c0_g1, were highly expressed only in the control group (Supplementary Table 3; Fig. 13c). KEGG enrichment analysis showed that eight genes, like TRINITY_DN38284_c0_g1, related to ion transport, were highly expressed in the 200 mM NaCl treatment while five genes, such as TRINITY_DN30952_c0_g2 and TRINITY_DN16375_c0_g1, were highly expressed in the control group (Supplementary Table 4). DEGs included cation transporters (TRK, HKT), sodium-related transporters (chaA, CAX), potassium transporters (kup), magnesium transporters (MRS2, MFM1), and others (Fig. 13d).
Transcription factor analysis
According to the PlantTFDB database, a total of 972 TFs were annotated. These belonged to 33 TF families, mainly MYB (157; 16.15%), AP2/ERF (102; 11.01%), bHLH (78; 8.02%), NAC (70; 7.20%), C2C2 (70; 7.20%) and others (Fig. 14).
Verification of target gene qRT-PCR
To verify the reliability of RNA-seq results, eight DEGs were screened with TRINITY_DN1091_c0_g1 as internal reference gene. TRINITY_DN22392_c0_g1, TRINITY_DN24475_c0_g2, TRINITY_DN43752_c0_g1 and TRINITY_DN4509_c0_g1 were those genes displaying upregulated expression while TRINITY_DN2820_c0_g2, TRINITY_DN8933_c0_g1, TRINITY_DN4277_c0_g1, and TRINITY_DN29407_c0_g2 were genes displaying downregulated expression (Supplementary Table 5). These genes were validated by qRT-PCR analysis. The expression of these eight genes based on qRT-PCR corresponded to their trends in expression as shown by RNA-seq sequencing, indicating that the RNA-seq sequencing results were valid and reliable (Fig. 15).