Transcriptome analyses revealed molecular responses of Cynanchum auriculatum leaves to saline stress

Cynanchum auriculatum is a traditional herbal medicine in China and can grow in saline soils. However, little is known in relation to the underlying molecular mechanisms. In the present study, C. auriculatum seedlings were exposed to 3.75‰ and 7.5‰ salinity. Next, transcriptome profiles of leaves were compared. Transcriptome sequencing showed 35,593 and 58,046 differentially expressed genes (DEGs) in treatments with 3.75‰ and 7.5‰, compared with the control, respectively. Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of these DEGs enriched various defense-related biological pathways, including ROS scavenging, ion transportation, lipid metabolism and plant hormone signaling. Further analyses suggested that C. auriculatum up-regulated Na+/H+ exchanger and V-type proton ATPase to avoid accumulation of Na+. The flavonoid and phenylpropanoids biosynthesis pathways were activated, which might increase antioxidant capacity in response to saline stress. The auxin and ethylene signaling pathways were upregulated in response to saline treatments, both of which are important plant hormones. Overall, these results raised new insights to further investigate molecular mechanisms underlying resistance of C. auriculatum to saline stress.

Tuberous root of Cynanchum auriculatum is a traditional herbal medicine in Asian countries. It has been listed in Chinese Pharmacopoeia and Korean Pharmacopoeia, namely "Baishouwu" and "Baekshuoh", respectively 1 . Pharmacological researches revealed that application of C. auriculatum could enhance immunity 2 , display activities of anti-tumor 3 , antioxidant 4 and gastro protection 2 . Thus, C. auriculatum has high economic values.
In China, C. auriculatum is cultivated in Shandong, Jiangsu, and Anhui provinces. Approximately 95% of C. auriculatum are produced at Binhai County, Yancheng City (China), which is a coastal city. Due to human activities and natural seawater erosion, most agricultural lands there show severe salinization. Cultivation on saline soils and/or using brackish water resources has attracted widespread attentions in recent years. Planting C. auriculatum on saline soils may be an alternative approach to solve the problem of insufficient agricultural lands, thus increasing the total production and decreasing the unit price of C. auriculatum. However, there are no reports investigating responses of C. auriculatum to saline stress.
The genus Cynanchum (Linn.) includes approximately 200 species in the world, which are widely distributed in eastern Africa, Mediterranean region, tropical, subtropical and temperate regions of Eurasia 5 . There are 53 species and 12 varieties of this genus in China, which are mainly distributed in southwest provinces 6 . The complete chloroplast genome sequence has been sequenced for C. auriculatum and C. wilfordii, and phylogenetic analysis revealed that these two species are evolutionarily close 7,8 . According to previous studies, plants in Cynanchum showed strong tolerance to salinity. It was found that the seed germination and radicle length of swallowwort (Cynanchum acutum L.) decreased with increased saline stress, which is more tolerant to salinity than common milkweed, hairy beggarticks, and scotch thistle 9 . As a halophyte speices, Cynanchum chinense increases the content of osmotic regulators such as soluble sugar, betaine and organic acids in saline-alkali environments, and its osmotic regulators content is higher than plants in Poaceae, Leguminosae, Cyperaceae and Boraginaceae families 10 . C. auriculatum is closely related to C. acutum and C. chinense 11 , and may also have similar tolerance to saline stress.
Besides antioxidant regulation, plants have also evolved other physiological and molecular mechanisms to resist saline stress 12 . For example, plant hormones, such as abscisic acid (ABA), auxin and ethylene, could function as essential endogenous signal molecules and regulate plant development and tolerance to salinity 13,14 . The treatments with salinity. Based on the preliminary results, saline treatments included two salinities, 3.75‰ and 7.5‰ (dissolving commercial sea salts in 1/2 Hoagland's solution). The experiments were carried out in plastic boxes (20 cm * 10 cm * 8 cm). In each box, 1 L of solution was added. For saline treatments, suitable amount of sea salts were added in 1/2 Hoagland's solution to achieve salinities of 3.75‰ and 7.5‰. 1/2 Hoagland's solution without extra addition of sea salts was also prepared as the control (CT). Each treatment consisted of five individuals and repeated three times independently. To avoid effects of evaporation on salinity, the culture media were renewed every three days. After 15 days, the top four leaves from five individuals in each treatment were collected and mixed as one sample. The plants were taken out, and quickly dried using sterile filter paper. After freezing with liquid nitrogen, leaves were stored at −80 °C until RNA extraction. transcriptome sequencing. The Biozol reagent produced by Bioer Technology (Hangzhou, China) was used to extract total RNA according to the manufacture's protocol. The integrity of RNA was visualized by 1% agarose gel electrophoresis, and its concentration and purity were detected using a NanoDrop Microvolume Spectrophotometers and Fluorometer (ThermoFisher, Shanghai, China). Oligo(dT) magnetic beads (NEB, USA) were used to enrich mRNA in the samples and then sequencing libraries were prepared using NEBNext mRNA Library Prep Master Mix Set for Illumina (NEB, USA) as described in the manufacture's protocol. AMPure XP beads (Beckman, Germany) were used to purify the DNA fragments with 250-300 bp length, which were further amplified by PCR using Phusion High-Fidelity DNA polymerase. Qubit 2.0 software was used for preliminary quantitative analysis of the libraries. The library was diluted to 1.5 ng/μL to check the distribution of insert size using an Agilent Bioanalyzer 2100 system. Real-time quantitative PCR (RT-qPCR) was applied to determine the accurate concentration of library (concentration > 2 nmol was the effective concentration of the library). Finally, Illumina HiSeq. 4000 platform was used to sequence the DNA libraries.

Bioinformatics analyses.
To determine whether sequencing data were suitable for subsequent analyses, quality control (QC) of raw reads was performed first. Reads with low quality (with >20% bases having Phred quality score ≤10), joint contamination and/or high ratio of unknown base (ratio ≥5%) were removed. Clean reads were assembled using the Trinity program 25 , and the longest transcript of each sequence was selected for the follow-up analysis of unigene. After assembly, the gene expression levels were compared between treatments and the control by calculating FPKM (expected number of fragments per kilobase of transcript sequence per million base pairs sequenced) values of each unigene using RSEM v1.2.8 26 . Genes with fold change ≥ 2 and Q-value (adjusted P-value) ≤ 0.001 were considered as significantly differentially expressed genes (DEGs).
TransDecoder software was used to identify candidate coding sequence (CDS). The longest open reading frame was extracted for each unigene, and then the Pfam protein homologous sequences were searched by blasting against SwissProt and Hmmscan database to predict CDS 27 . For gene ontology (GO) annotations, the DEGs were blasted against the GO database using Blast2Go program (http://www.blast2go.com/Ver.2.3.5) 28 . For Kyoto Encyclopedia of Genes and Genomics (KEGG) annotations, all the DEGs were mapped to the KEGG database (https://www.genome.jp/kegg/pathway.html) using BLASTX 29 . Comparisons with Q-value ≤ 0.05 were considered statistically significant.

Real-time quantitative PCR.
To verify the reliability of RNA-seq results, nine unigenes were randomly selected for RT-qPCR. The total RNA samples used for qPCR were the same for RNA-seq. The first strand of cDNA was synthesized by reverse transcription using the BioRT cDNA first strand synthesis kit (Bioer, Hangzhou, China) with oligo(dT) primer (Supplementary Table S1). Beta-actin, which showed similar FPKM values among treatments, was also included as the internal control for qPCR. The qPCR experiments were performed in 20 µl of reaction system using BioEasy master mix (Bioer, Hangzhou, China). The thermal cycles were performed and the signals of SYBR Green were detected using a Gene9600 Plus qPCR machine (Bioer, Hangzhou, China). Melting curves were prepared to reveal the specificity of primers. Expression levels of each gene were compared by calculating their relative fold change using the 2 −ΔΔCt method 30 .

Results and Discussion
Preliminary experiments were conducted to test the tolerance of C. auriculatum to saline stress. Seedlings were treated with salinities ranging from 7.5‰ to 15‰. After treatment for one week, the leaves of C. auriculatum turned yellow, the whole plant was shorter and the bottom leaves withered in the 7.5‰ treatment, compared with the control. In treatments with 10‰, 12‰ and 15‰, the plants died (Fig. S1). These results indicated that C. auriculatum could tolerate 7.5‰ salinity. To explore the underlying molecular mechanisms, treatments with 3.75‰ and 7.5‰ were re-prepared for transcriptome sequencing.
illumina sequencing. The sequencing data have been deposited in National Center for Biotechnology Information (NCBI) with the bioproject number of PRJNA558052. RNA sequencing resulted in 65 M to 73 M of clean reads for each leaf sample. Q20 and Q30 values were all higher than 96.92% and 88.19%, respectively (Table S2). These results suggested that the quality of sequencing data could satisfy further transcriptome analyses.
De novo assembly. Separate assembly of each sample was tried first. The results showed that the total number of unigenes ranged from 43,123 to 76,090, and the mean length of unigenes ranged from 1,257 bp to 1,521 bp. Alternatively, all sequencing data were assembled together, finally producing 196,199 unigenes with the mean length of 1,673 bp (Table S3). Among them, 15.89% unigenes were shorter than 300 bp and 17.07% were longer than 3,000 bp (Fig. S2). Benchmarking Universal Single-Copy Orthologs (BUSCO) 31 analyses revealed that 99% BUSCOs were complete and 1% BUSCOs were missing (Fig. S3). Due to the longer mean length and higher ratio of BUSCOs, the unigene library assembled from nice samples together was used for the subsequent analyses. prediction of coding sequences and functional annotation. Prediction of coding sequences detected 122,855 CDSs in total. Their total length and N50 length was 137.21 M bp and 1,431 bp, respectively (Table S4). Among them, 77.44% CDSs were longer than 500 bp, and 42.61% CDSs were longer than 1,000 bp (Fig. S4).
Differentially expressed genes and qPCR validation. Pairwise comparisons revealed that 16,849 and 19,404 unigenes were significantly up-regulated, 18,744 and 38,642 were down-regulated in treatment with 3.75‰ and 7.5‰, respectively, compared with the control. After removing unigenes with FPKM values lower than 2 in all samples, 4,975 and 6,566 unigenes were significantly up-regulated, 4,820 and 6,484 were significantly down-regulated in treatment with 3.75‰ and 7.5‰, respectively (Fig. 1a). Among them, 2,897 up-regulated DEGs and 2,956 down-regulated DEGs were commonly shared between treatments with 3.75‰ and 7.5‰ (Fig. 1b). The top 30 DEGs are listed in Table S5.
To validate the DEGs identified by RNA-Seq, nine randomly selected DEGs were determined by qPCR. The results showed similar changing trends between qPCR validation and FPKM calculation for all selected genes ( Fig. S5), indicating that the transcriptome data were reliable.

Enrichment of GO categories for DEGs.
In both treatments with 3.75‰ and 7.5‰, the top four GO categories were 'membrane part' , 'cell' , 'catalytic activity' and 'cellular process' (Fig. 2). In these categories, more DEGs were included in treatment with 7.5‰ than 3.75‰, suggesting that higher salinity caused more severe effects on C. auriculatum leaves. Besides, the category of 'binding' was only highly expressed in treatment with 3.75‰ but not in treatment with 7.5‰.
Enrichment of KEGG pathways for DEGs. KEGG enrichment analysis of down-regulated genes in treatment with 3.75‰ and 7.5‰ revealed 19 and 28 pathways, respectively (Tables 2 and 3). These pathways might represent the inhibitory effects on the molecular processes in C. auriculatum leaves. Analysis of up-regulated genes, which represented mechanisms underlying resistance to saline stress, showed 19 and 23 significantly enriched pathways at 3.75‰ and 7.5‰, respectively, which shared 12 KEGG pathways (ko00500, ko00940, ko04075, ko00563, ko00941, ko00330, ko00511, ko00591, ko04016, ko00270, ko00460 and ko00904; Tables 2 and 3). These enriched pathways mainly involved in oxidative stress-related function, signal transduction, carbohydrate metabolism and lipid metabolism.  www.nature.com/scientificreports www.nature.com/scientificreports/ ion transportation. Under saline conditions, the first response of plants is ion toxicity, and generally Na + displays more severe influence than Cl − , since Na + would reach toxic concentrations earlier than Cl − . Thus, most studies focused on Na + uptake, transportation and exclusion in plants 28 . In the present study, compared with the control, Na + /H + exchanger was significantly up-regulated in treatment with 3.75‰ (P < 0.05), but K + efflux antiporter, cation/H + antiporter and sodium/proton antiporter were significantly down-regulated at 7.5‰ salinity (P < 0.05). V-type proton ATPase (V-ATPase) was significantly up-regulated, while cyclic nucleotide-gated ion channel (CNGC) were significantly down-regulated at 3.75‰ salinity (P < 0.05) ( Table 4). The up-regulation of Na + /H + exchanger might alleviate cytosolic salt accumulation in the vacuole and other cellular compartments 32,33 . Similar results were observed in Arabidopsis 34 , which assisted Na + extrusion from plant cells. Activation of V-ATPase in plants might reduce Na + accumulation by interacting with H + pumps to resist saline stress 35 . CNGCs are non-selective ion channels involved in mediation of Na + transport. The expression level of CNGC significantly decreased at 3.75‰ salinity compared with the control, which was in agreement with the findings reported in Arabidopsis that CNGC10 expression was inhibited when exposed to high salinity or low salinity for long time 36 . K + efflux antiporter was significantly down-regulated at 7.5‰ salinity, similar to the transcriptome analysis of grapevine (Vitis vinifera L.) under saline conditions 37 , implying the replacement of K + by Na + .   38 . Among them, SOD is an important antioxidant enzyme 39 . SODs have three isozymes including copper-zinc SOD (Cu/Zn-SOD), manganese SOD (Mn-SOD) and iron SOD (Fe-SOD). In the present study, three isozymes were all detected. In comparison to the control, Mn-SOD was significantly up-regulated for 2.15 times in 7.5‰ treatment, while Fe-SOD increased significantly for 1.91 times at 3.75‰ salinity (P < 0.05) ( Table 4). These changes might enhance tolerance of C. auriculatum to saline stress, since it has been reported that overexpression of Mn-SOD enhanced salt-tolerance in Arabidopsis 40 . In pea, the activity of Mn-SOD in salt-tolerant plants is higher than that in salt-sensitive plants under NaCl treatment. Further analysis found that NaCl-tolerant plants had protective mechanisms against salt-induced O 2 − production by increasing the mitochondrial Mn-SOD activity 41 . According to these results, C. auriculatum might be salt-tolerant, but the specific functions of SODs need further study.
Flavonoid biosynthesis pathway plays a pivotal role and protects the plant cells from oxidative damage by scavenging free radicals 42,43 . KEGG enrichment analyses of up-regulated DEGs significantly enriched flavonoid biosynthesis pathway. Within this pathway, trans-cinnamate 4-monooxygenase (C4H) and chalcone isomerase (CHI) are key enzymes, directly related to the synthesis of flavonoids 44   www.nature.com/scientificreports www.nature.com/scientificreports/ control, CHI and C4H were significantly up-regulated in treatments with 3.75‰ and/or 7.5‰ (Table 4). These results implied that more flavonoids were synthesized, which might contribute to the total antioxidant capacity in response to saline stress in C. auriculatum. Similarly, Walia et al. 45 reported that a large number of genes were up-regulated in the flavonoid biosynthesis pathway under salinity stress, which played an important protective role against saline stress.
Phenylpropanoids are involved in antioxidant activity in cell walls and lignin biosynthesis in secondary metabolites, and play an important role in plant response to abiotic stress, which connects to flavonoid biosynthesis 46 . Phenylalanine ammonia-lyase (PAL) is the first enzyme that converts phenylalanine into cinnamic acid and is transcriptionally regulated in response to many environmental cues 47 . PAL was a marker during molecular breeding to improve resistance of plants to injury and saline stress 48 . In the present study, PAL was significantly up-regulated for 1.47 and 2.13 times at 3.75‰ and 7.5‰ salinity (P < 0.05), respectively, compared with the control (Table 4). Similar findings were observed in salinity-stressed Olea europaea 49 .
Carotenoids reveal antioxidant activity in plants. In the carotenoid biosynthesis pathways, phytoene synthase (PSY) is a rate-limiting enzyme, related to abiotic stress and mainly involved in plant stress response 50 . In the present study, compared with the control, the expression level of phytoene desaturase (PDS) was significantly  www.nature.com/scientificreports www.nature.com/scientificreports/ lower in treatment with 7.5‰ and that of nine-cis-epoxycarotenoid dioxygenase (NCED) was significantly lower in treatment with 3.75‰. Besides, the expression level of phytoene synthase (PSY) was lower in saline treatments than that in the control, although no statistical significance was detected (P > 0.05; Table 4). Similarly, PSY1 gene was reduced in response to enhanced NaCl level 51 , and NCED gene was clearly repressed under drought and saline stress in cotton 52 . The decreased expression of key enzymes in carotenoid biosynthesis pathway revealed that carotenoid biosynthesis might not be the antioxidant mechanism in salinity-stressed C. auriculatum.

Saline treatments regulated lipid metabolism in C. auriculatum. Six lipid metabolism-related
pathways were significantly enriched when up-regulated or down-regulated genes were subjected to KEGG enrichment analyses, including alpha-linolenic acid (LnA) metabolism, linoleic acid (LA) metabolism, glycerophospholipid metabolism, fatty acid biosynthesis, fatty acid degradation, and sphingolipid metabolism. In response to adverse factors, alpha-linolenic acid and linoleic acid were released and became the initial step of unsaturated fatty acid metabolism 53,54 . LnA was catalyzed by different functional lipoxygenase (LOX) to produce hydroperxy linolenic acid (HPOT), and then metabolized into a series of compounds with certain biological activities, such as jasmonic acid (JA), 3-hexenol and traumatic acid, which related to plant defense responses to abiotic stresses [54][55][56] . Previous studies have shown that JA biosynthesis genes were up-regulated in Arabidopsis and sunflower roots in saline treatments 57,58 . In the present study, compared with the control, LOX and hydroperoxide dehydratase were significantly down-regulated at 3.75‰ (P < 0.05) ( Table 4), revealing that unsaturated fatty acid metabolism might be inhibited by saline stress in C. auriculatum. changes of plant hormone signaling transduction pathways under saline stress. Plants respond to abiotic stress through a complex coordination of various phytohormone signaling pathways, mainly including auxin (AUX), cytokinine (CTK), gibberellin (GA), abscisic acid (ABA), ethylene (ETH), brassinosteroid (BR), jasmonic acid (JA) and salicylic acid (SA) 59 . In the present study, the functional analysis of DEGs by KEGG enrichment analysis showed that a large number of genes were involved in plant hormone signaling transduction and hormone responses under saline stress, highlighting the key roles of plant hormones in regulating C. auriculatum response to saline stress. In the JA signal pathway, the expression level of jasmonate ZIM domain-containing protein (JAZ) was significantly down-regulated at 3.75‰ salinity (Table S6), indicating that JA signaling transduction was inhibited by saline stress. This result was consistent with the inhibition of unsaturated fatty acid metabolism.
As previously reported, the BR signaling pathway was activated under saline stress in Arabidopsis 60 , and the exogenous application of BR effectively alleviated the adverse effects of abiotic stress or increased plant resistance to stresses 61,62 . In the present study, compared with the control, brassinosteroid insensitive 1 (BRI1), BR-signaling kinase (BSK), protein brassinosteroid insensitive 2 (BIN2), brassinosteroid resistant 1/2 (BZR1/2), xyloglucan: xyloglucosyl transferase (TCH4) and cyclin-D3 (CYCD3) were significantly up-regulated in treatments with 3.75‰ and/or 7.5‰ (Table S6), suggesting that activation of BR signaling pathway might participate in resistance to saline stress.
Auxin can induce expression of auxin-response factor gene (ARF), auxin early response gene (Aux/IAA), gretchen hagen3 (GH3) and small auxin-up RNAs (SAUR) rapidly and instantaneously 63 . In the present study, compared with the control, GH3 and SAUR were significantly up-regulated in response to saline treatments. However, the expression levels of protein transport inhibitor response (TIR1) was significantly down-regulated under saline stress (P < 0.05; Fig. 3). Since GH3 and SAUR were downstream genes of auxin signaling pathway, the present results suggested that auxin signaling pathway was activated in C. auriculatum in response to saline stress.
The abscisic (ABA) signaling pathway is well known to be associated with responses to saline stress in plants and the increase of ABA content could help plants adapt and survive under saline condition by reducing the accumulation of Na + and improving osmotic adjustment 15 . Besides, carotenoid biosynthesis pathway provides precursors for plant hormone abscisic acid (ABA) synthesis 64 . In the present study, ABA receptors (PYR/PYL) were significantly down-regulated in treatment with 7.5‰ (Fig. 3), which might be resulted from the inhibition of carotenoid biosynthesis. Similar results were also detected in Corchorus spp. and Vitis vinifera L 37,65 . Besides, downstream of PYR/PYL, serine/threonine-protein kinase 2 (SNRK2) was also significantly down-regulated for 2.03 and 1.97 times in treatment with 3.75‰ and 7.5‰, respectively (P < 0.05; Fig. 3). Overall, these results suggested that ABA signaling pathway was inhibited by saline treatments in C. auriculatum.
Ethylene is a stress hormone regulating myriad stress responses 66 . Ethylene receptors (ETR) serve as negative regulators of the ethylene signaling transduction pathway, while ethylene-insensitive protein 2 (EIN2) and ethylene-insensitive protein 3 (EIN3) are positive regulators. Thus, a decrease in receptor levels is predicted to sensitize the plant such that it responds to lower levels of ethylene than usual 67 . In the present study, compared with the control, mitogen-activated protein kinase kinase (SIMKK) and ethylene-responsive transcription factor 1/2 (ERF1/2) were significantly up-regulated in treatment with 7.5‰, while serine/threonine-protein kinase CTR1 (CTR1) and EIN3-binding F-box protein 1/2 (EBF1/2) were significantly down-regulated in treatment with 3.75‰ and 7.5‰, respectively (P < 0.05; Fig. 3, Table S6). No significant differences in expression levels of ETR, EIN2 and EIN3 were detected between saline treatments and the control. Considering ERF1/2 is a downstream gene of CTR1 and EBF1/2 (Fig. 3), the present results implied that ethylene signaling pathway might be activated in C. auriculatum under saline conditions.
In summary, the changes of gene expression profile in C. auriculatum leaves under saline stress were explored using the transcriptome sequencing to investigate regulatory mechanisms in response to saline stress. Compared with the control, C. auriculatum showed suppressed carotenoid biosynthesis, lipid metabolisms, abscisic acid and jasmonic acid signaling transduction, which might damage plant cells. Expression levels of Na + /H + exchanger and V-type proton ATPase, the flavonoid and phenylpropanoids biosynthesis pathways, the auxin and ethylene signaling pathways were upregulated to improve the antioxidant capacity of C. auriculatum in response to saline treatments. These changes might facilitate resistance of C. auriculatum to saline stress. Flavonoid and phenylpropanoids appear to play roles in resistance to salt tolerance in C. auriculatum, possibly due to its pharmacodynamics.

Data availability
The original sequencing files have been deposited in the National Center for Biotechnology Information (NCBI) with the bioproject number of PRJNA558052.