Physiological changes and transcript identification in Coreopsis tinctoria Nutt. in early stages of salt stress

Plant materials and treatment

C. tinctoria seeds were collected from Hotan, China (altitude, 2,196 m, N: 37°37′17.00″E: 78°16′58.80″). The seeds were planted in an artificial-lighting hydroponic plant factory under 10 h of artificial light (8,100 1×  light intensity, 25–28 °C), 14 h of darkness (20–22 °C), and a relative humidity of 60–70%. A total of 300 undamaged seeds of the same size were selected, sterilized and sown in Petri dishes. After the seeds germinated for 15 days (2 leaves), seedlings with consistent growth were transferred to a 6-L hydroponic container containing 1/2 Hoagland nutrient solution for continuous cultivation. During this period, the nutrient solution was changed every 5 days, and an oxygen pump was used to supply oxygen regularly for 6 h a day. On the 30th day (when the plants had 4–6 leaves), 0 (control (CK)), 50 (FS), and 200 (TS) mM NaCl solutions were used as salt stress treatments. Leaves of C. tinctoria were quickly cut off at 12 h, 24 h, and 48 h (the second pair of functional leaves from bottom to top) and put into liquid nitrogen at −80 °C for cold storage. For each time point, three biological repetitions were conducted, with 30 seedlings per repetition.

Physiological index and growth determination

The MDA, Pro, soluble sugar and soluble protein contents were determined by the thiobarbituric acid (TBA) colorimetric method, indene triketone staining method, anthrone colorimetric method and Coomassie brilliant blue GMel 250 colorimetric method (Gao, 2006). Superoxide dismutase, POD, catalase (CAT), and ascorbate peroxidase (APX) were detected by Rayto RT-6100 enzyme colorimetry and an Elabscience plant kit (Wuhan, China).

The fresh weight and the root/shoot ratio of the plants were determined by the following method. The seedlings were removed from the nutrient solution, washed clean with deionized water, and quickly dried with absorbent paper. The above- and belowground fresh weights were measured, and the root/shoot ratio (underground fresh weight/aboveground fresh weight) was calculated (Rahimzadeh Karvansara & Razavi, 2019).

RNA Extraction, Library Establishment and Sequencing

RNA extraction, library establishment and sequencing were performed similarly to previously reported methods (Xue et al., 2020). The total RNA of 21 leaf samples of C. tinctoria was extracted by the TRIzol (Invitrogen) method, and the quality of RNA samples was detected by agarose gel electrophoresis, a Qubit 2.0 Fluorometer (Life Technologies, CA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) to ensure that qualified samples were used for transcriptome sequencing.

The cDNA library was established by a TruSeq^TM RNA sample preparation kit (Illumina, San Diego, CA). mRNA with a polyA tail was enriched with oligo(dT) magnetic beads, and the mRNA was broken into short fragments with the addition of fragmentation buffer. Using the short fragment mRNA as the template, first-strand cDNA was synthesized by six-base random primers (random hexamers). Buffer, dNTPs (dUTP, dATP, dGTP and dCTP) and DNA polymerase I were then added to synthesize double-stranded cDNA. The double-stranded cDNA was then purified by AMPure XP beads. The double-stranded cDNA structure had a sticky end that was filled until blunt via the addition of End Repair Mix, and an A base was added at the 3′ end to overlap with a T base overhang. The fragment size was then selected by AMPure XP beads, and the final cDNA library was enriched by PCR as described in the instructions. After cDNA enrichment by PCR, a Qubit 2.0 instrument was used for preliminary quantification, and an Agilent 2100 instrument was used to detect the insert size of the library. The effective concentration of the library was subsequently determined by Q-PCR (effective concentration of the library > 2 nM), and the Illumina HiSeq platform was used for sequencing (Wuhan Metware).

Data filtering and reassembly of the transcriptome

The raw data of 21 samples were obtained by de novo transcriptome sequencing, and the raw data underwent filtering and quality control before reassembly. SeqPrep (https://github) was used to remove the adaptor sequences and low-quality sequences (Q <20), sequences containing more than 10% N bases, and sequences less than 20 bases long to obtain clean reads. C. tinctoria has no reference genome, and the clean reads were spliced by Trinity (http://trinityrnaseq.sourceforge.net/) to obtain the transcript sequence. The transcript sequence was used as a reference sequence for subsequent analyses. Hierarchical clustering of transcripts was performed by comparing the read number and expression pattern of transcripts using Corset (https://code.google.com/p/corset-project/). Hierarchical clustering was conducted with Corset to ultimately obtain the longest cluster sequence as a follow-up analysis on unigenes (Conesa et al., 2016).

Functional annotation and classification of unigenes

To annotate and infer the functions of the transcripts, we used BLAST software to compare the unigene sequences with 7 public databases, including the NCBI nonredundant protein (NR, http://www.ncbi.nlm.nih.gov), TrEMBL (a computer annotated supplement to SwissProt), Pfam (http://xfam.org/), SwissProt (http://www.expasy.ch/sprot), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (http://www.genome.jp/kegg). Gene Ontology (GO, http://www.geneontology.org) classification based on molecular function (MF), biological process (BP), and cell component (CC) and euKaryotic Ortholog Groups (KOG) classification using the KOG database (http://www.ncbi.nlm.nih.gov/COG/) were also performed.

Quantitative gene expression and differentially expressed gene (DEG) analysis

The clean reads of 21 samples were mapped back to a reference sequence using Bowtie 2 (http://bowtie-bio.sourceforge.net), and the fragments per kilobase of transcript per million fragments mapped (FPKM) method was used to normalize the expression of mapped reads. The read count of genes was realized using feature counts (Liao, Smyth & Wei, 2014). DESeq2 was used to analyze the difference in read count data. The multiple hypothesis test correction of the P value (probability) was carried out using the Benjamini–Hochberg method, and the false discovery rate (FDR) was obtained. Genes that satisfied the conditions of — log₂FC — ≥ 1 and FDR < 0.05 were differentially expressed. GO enrichment analyses of DEGs were performed using the GO database (Wang, Gerstein & Snyder, 2009; Love, Huber & Anders, 2014).

Weighted gene coexpression network analysis (WGCNA)

WGCNA is used to find highly related and hub gene modules that play important roles. The present study analyzed the physiological indexes and all DEGs of C. tinctoria seedlings under 21 low-salt and high-salt treatments using WGCNA. The first 100 groups of genes with the highest edge weight of WGCNA in the gene module with high correlation with physiological traits were selected and visually analyzed by Cytoscape 3.7.1 to determine the hub gene (Lou et al., 2017).

QRT-PCR analysis

Total RNA was extracted from leaves using TRIzol reagent (Invitrogen). A NanoDrop 2000 was used to test the concentration and purity of the RNA. A Servicebio RT First-Strand cDNA Synthesis Kit for reverse transcription was used to synthesize cDNA to serve as the template for qRT-PCR. Then, 2 × SYBR Green qPCR Master Mix (High ROX) was used for quantitative PCR (Wuhan Servicebio Technology Co., Ltd). Primers targeting the reference gene GAPDH were designed according to the conserved region in a GAPDH gene sequence alignment of Helianthus annuus. The 2^−ΔΔCT method was used to calculate the relative expression of genes. Each gene was examined in three biological replicates. Primer information is shown in the attached Table S1.

Statistical analysis

An analysis of variance was carried out using SPSS 23.0, and the mean values were separated using the least significant difference (LSD) test at P = 0.05. All figures were drawn using Origin 2018. Pearson’s correlation and DEG heat map analyses were achieved via the functions “cor.test” and “heatmap”, respectively.

Changes in the physiology and growth of C. tinctoria under salt stress

Salt stress changed the contents of MDA, proline, soluble sugar and soluble protein in C. tinctoria seedlings

To understand the physiological and growth changes of C. tinctoria under the 50 and 200 mM NaCl treatments, the osmotic regulation substances, antioxidant enzymes and biomass of seedling leaves were quantified within 72 h. The MDA content in the 50 and 200 mM NaCl samples first increased then stabilized over time. The maximum values appeared at 48 h and were 1.05 and 1.87 times higher than that of the CK, respectively (Fig. 1A). The overall change trend of Pro and soluble sugar contents from 0 to 72 h after 50 mM NaCl treatment was the same, with an initial increase followed by a decreasing trend. Their maximum values were reached at 48 h and were 1.40 and 1.26 times higher than CK, respectively. However, the Pro and soluble sugar contents increased gradually after 200 mM NaCl treatments and reached maximum levels at 72 h, which were 1.11 and 0.41 times higher than CK, respectively. The soluble protein content increased at different times after 50 and 200 mM NaCl treatments. Soluble protein increased significantly at 48 h after 50 mM NaCl treatments and was 1.62 times higher than that of the CK. No obvious change in the content of soluble protein was found after 200 mM NaCl treatments, and the content was highest at 12 h and 48 h, which was not significantly different from 50 mM NaCl treatments for 48 h (Figs. 1B–1D). The results showed that the Pro, soluble sugar and soluble protein contents in early C. tinctoria seedlings increased significantly under low and high salt stress but increased more under high salt stress, which indicates the positive roles of osmotic and organic regulators in salt stress, especially under high salt stress (200 mM NaCl).

Salt stress changed the antioxidant enzyme content of C. tinctoria seedlings

The change trends of POD in the 50 and 200 mM NaCl treatments increased initially then stabilized with the increase for the 200 mM NaCl samples and was slightly higher than the 50 mM NaCl samples. The change trends of CAT were exactly the opposite. CAT in the 50 mM NaCl treatments decreased at 12 h then gradually increased. However, CAT in the 200 mM NaCl treatment increased at 12 h to a level 1.32 times higher than that of the CK and then decreased to a stable level. Under 50 and 200 mM NaCl treatments, the SOD content increased significantly at 12 h, but the content of other treatments was lower than that of the CK. The APX content increased significantly from 12 to 24 h under 50 mM NaCl treatments and was significantly lower than that of the CK after 48 h. Under 200 mM NaCl treatments, the content of APX was significantly lower than that of the CK at different times (Figs. 1E–1H). Taken together, these results showed that low and high salt stress promoted the activities of POD and CAT to different degrees, and high salt stress significantly inhibited the activities of SOD and APX.

Salt stress changed the biomass of C. tinctoria seedlings

The aboveground fresh weight, underground fresh weight and root/shoot ratio of C. tinctoria seedlings were measured 72 h after salt stress. The biomass and root/shoot ratio of C. tinctoria seedlings increased to varying degrees under 50 mM NaCl treatments compared to the CK treatment. The aboveground fresh weight, underground fresh weight and root/shoot ratio increased 31.44%, 47.22% and 10.13%, respectively. Under the 200 mM NaCl treatments, the biomass decreased 48.26%, 46.87% and 0.27%, respectively. These results showed that low and high salt significantly affected the growth of C. tinctoria. The effect of 50 mM NaCl treatments on the aboveground plant part was less than the underground plant part, and the effect of 200 mM NaCl treatments on the aboveground and underground parts was the same.Low salt (50 mM NaCl) promoted the growth of C. tinctoria seedlings, and high salt (200 mM NaCl) significantly inhibited the growth of C. tinctoria seedlings (Fig. S1).

Transcriptome sequencing analysis

To further understand the response of C. tinctoria morifolium seedlings to low and high salt treatments, de novo transcriptome sequencing was performed on the leaves of C. tinctoria after 50 (FS) and 200 (TS) mM NaCl treatments for 0 (CK), 12 h, 24 h and 48 h. cDNA libraries of 21 transcriptome sequences were constructed. For each sample, 41,630,370–50,805,280 raw reads were obtained from high-throughput sequencing, and after splices and low-quality sequences were excluded, 40,125,080–49,057,376 clean reads were obtained. Q30 bases accounted for greater than 93% of the total, and the average GC content was 45.49%. The data for each sample are summarized in Table S2. After assembly, 507,230 transcripts with an average length of 678 bp and 428,638 unigenes with an average length of 750 bp were identified. The N50 values were 864 bp and 922 bp, respectively (Table S3). A total of 14,588 transcripts and unigenes, with an average sequence length of more than 2,000 bp, were found (Fig. S2).

Comments and functional classification of unigenes

We compared a total of 428,638 unigenes with information from 7 public databases (the NR, TrEMBL, SwissProt, Pfam, KOG, GO, and KEGG databases) and found that all unigenes shared a similarity of 63.76% in at least one database (Table S4). A comparison of 270,860 genes with the NR database showed that 76.31% (206,689) of the genes had the greatest sequence similarity with Helianthus annuus and low homology with other species (Fig. S3A). Based on the classification of 227,019 genes successfully annotated to one or more GO terms in the GO database (Table S3), all genes were divided into 3 GO categories and 59 groups. The most obvious matching GO terms in the “biological process” category were cellular process (134,866), metabolic process (115,375) and response to stimulus (60,937). The main terms in “cell components” were cell (154,999), cell part (154,612), and organelle (117,183). For the category of “molecular function”, GO binding (136,084), catalytic activity (121,727), and transporter activity (14,544) had the largest number of genes (Fig. S3B). The most common homologous groups were general function prediction only (33,064), signal transduction mechanisms (17429), posttranslational modification, protein turnover, chaperones (15,911) and transcription (8,049) in the KOG database for linear gene classification and functional prediction of 146589 unigenes (34.20%) (Fig. S3C). We further mapped unigenes (203,045) to the KEGG database for annotation and enrichment analysis. All genes were divided into five biochemical pathways, cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems, with the largest proportion of unigenes in the metabolic pathway. The metabolic pathway terms metabolic pathways (38,578), biosynthesis of secondary metabolites (18,406) and carbon metabolism (5,186) involved the largest number of genes (Fig. S3D).

Gene expression and DEG analysis

To evaluate the overall trends of gene expression, we used a box plot with the expression level in terms of FPKM on the horizontal axis. The expression levels of most genes were consistent, and only a few genes were distinct (Fig. S4A). The Pearson’s correlation coefficients of the expression levels of each library showed a high correlation among biological repeats in each group of 50 and 200 mM NaCl treatment treatments (R² > 0.9) and a significant difference in gene expression between groups (Fig. S4B).

The selected DEGs (— log₂FC) — ≥ 1 and FDR < 0.05) were then compared and analyzed. In total, 3,679, 1,761, 3,772, 9,055, 6,533, and 4,110 upregulated genes and 4,905, 1,999, 4,061, 10,286, 6,700 and 5,114 downregulated genes at 12 h, 24 h and 48 h for the 50 and 200 mM NaCl treatments compared with the CK treatmen were detected (Fig. 2). The number of DEGs under 200 mM NaCl treatments decreased gradually over time, and the number of DEGs under 50 mM NaCl treatments decreased initially then increased (Fig. 2). The number of DEGs peaked at 12 h in both treatments. DEG analysis showed that the DEGs were downregulated slightly under 50 and 200 mM NaCl treatments, and the gene expression changed more actively within 24 h under the 200 mM NaCl treatment (Fig. 2). A Venn diagram showed that 3,603 DEGs were expressed under 200 mM NaCl treatments at different times (Fig. 3A), which was much higher than the coexpressed DEGs in the 50 mM NaCl treatment groups (Fig. 3B). A total of 306 DEGs were expressed in all the comparison groups (Fig. 3C).

Identification and functional annotation of the physiological indicators related to DEGs by WGCNA

A total of 34 DEG modules were identified by WGCNA, four of which (blue2, Darkolivegreen, Lightskyblue4 and Darkseagreen) highly correlated with physiological indicators (Pearson’s r > 0.7 or r < −0.7, p value < 0.01) (Figs. 4A–4B). To determine the transcriptional expression of these modules, heat maps of the genes in these four modules were drawn (Fig. S5). Under salt stress, the genes in module blue2 had important effects on Pro (r = −0.74, p value = 1.0 ×10⁻⁴), soluble sugar (r = - 0.88, p value = 2.0 ×10⁻⁷), soluble protein (r = −0.78, p value = 3.0 ×10⁻⁵) and POD (r = −0.79, p value = 2.0 ×10 ⁻⁵). The heat map of the gene expression in the blue2 module showed that gene expression in the CK and 50 mM NaCl treatments was significantly higher than that in the 200 mM NaCl treatment (Fig. 4B, Fig. S5A), which indicated that blue2 was a specific coexpression module in response to low salt stress. The Darkolivegreen module highly negatively correlated with soluble sugar (r = −0.77, p value = 5.0 ×10⁻⁵) and POD (r = −0.80, p value = 1.0 ×10⁻⁵). The expression of genes in this module was the highest in the CK and 50 mM NaCl treatment groups (Fig. 4B, Fig. S5B). The Darkseagreen module highly positively correlated with MDA (r = 0.78, p value = 3.0 ×10 ⁻⁵). The heat map of gene expression indicated that the gene expression levels in 200 mM NaCl groups after 24 h and 48 h of treatment were higher than the other treatments (Fig. 4B, Fig. S5C). The genes in the Lightskyblue4 module had an important influence on CAT (r = 0.75, p value = 9.0 ×10⁻⁵). The expression of genes in Lightskyblue4 under 200 mM NaCl treatments at different times was significantly higher than that under 50 mM NaCl treatments (Fig. 4B, Fig. S5D), which indicates that Lightskyblue4 was a specific coexpression module in response to high salt stress.

The interaction network of the top 100 genes with the highest weight of WGCNA that were specific to the blue2 and Lightskyblue4 modules was visualized using Cytoscape. Twenty-five genes with WGCNA weights > 0.19 were highly connected in the blue2 module, and 9 genes had more than 10 edges (Fig. 5A). Thirty-seven genes with WGCNA weights > 0.31 were highly connected in the Lightskyblue4 module, and 6 genes had more than 10 edges (Fig. 5B).

To explore the biological functions of these genes, the most significantly enriched GO terms in the blue2 and Lightskyblue4 modules were first identified (Table S5). Five BP GO terms, including reductive pentose-phosphate cycle, photosynthesis, dark reaction, photosynthetic electron transport chain, photosynthetic electron transport in photosystem I and plastid membrane organization, were significantly enriched in the blue2 module. The annotation of the 25 genes is shown in Table S6, and it mainly includes 15 enzyme genes, 5 protein-related genes, and 5 genes with unknown functions. The genes in the Lightskyblue4 module were significantly enriched in proline biosynthetic process, L-proline biosynthetic process, proline metabolic process, peptidyl-threonine phosphorylation, and peptidyl-threonine modification. The annotations of 37 genes, including 22 enzyme genes, 8 protein-related genes and 7 genes with unknown functions, are shown in Table S7.

Related genes responding to physiological changes in C. tinctoria under salt stress

Twenty and 30 coexpressed hub genes related to physiological changes in the blue2 and Lightsskyblue4 modules were screened, respectively. In the blue2 module, 9 RuBPcase, 1 RpiA, 1 GS, 1 CHY, 1 GAPD, and 1 MDH genes were found to be involved in the regulation of C. tinctoria metabolism and other biological processes. In the Lightskyblue4 module, 4 P5CS encoding key enzymes, 3 PSA, 2 P5CDH, 2 RS, 2 AASS, 2 AGT, 1 UGE, 1 G6PDH, 1 INF2, 1 PP2C, 1 IRAK-4, and 1 E3-UBPL genes were detected. These DEGs played important regulatory roles in the antioxidant system under high salt (200 mM NaCl) conditions. The hub gene expression of 20 genes in blue2 and 30 genes in Lightskyblue4 that were related to changes in physiological indicators are shown in Fig. 6. The genes in the blue2 module were highly expressed in CK and under the 50 mM NaCl treatment and included 4 RuBPcase, 1 RpiA and 1 CHY genes located in the core position of the network diagram. These genes were expressed at a higher level after 48 h of treatment. The accumulation of physiological indicators (except SOD) that were highly related to the blue2 module in the 200 mM NaCl treatment group was higher than that in the 50 mM NaCl treatment group. However, the total expression of the hub genes was contrary to the gene accumulation, which indicates a negative regulatory relationship between them. Genes in Lightskyblue4 were highly expressed under 200 mM NaCl treatments, especially after 12 h of treatment. Among these genes, 2 PSA, 2 P5CDH and 1 E3UBP genes were in the core position of the network diagram. The accumulation of the total expression of these key genes in the 200 mM NaCl treatments was significantly higher than that in the 50 mM NaCl treatments, which was consistent with the content changes in CAT (Fig. 6).

Verification of gene expression by qRT-PCR

To confirm the expression patterns of DEGs obtained from RNA-seq, 12 DEGs involved in flavonoid biosynthesis were randomly selected for analysis by real-time quantitative reverse transcription PCR (qRT-PCR) (Fig. S6). For all selected genes, the expression trend from qRT-PCR was consistent with that from the RNA-seq data, which indicates that the RNA-seq data re reliable.