Physiological Differences and Transcriptome Analysis Reveal That High Enzyme Activity Significantly Enhances Drought Tolerance in Chinese Fir (Cunninghamia lanceolata)

Abstract

Chinese fir (Cunninghamia lanceolata) is the most cultivated timber species in China, with a plantation area of 11 million ha. Due to its extensive geographical distribution, drought stress caused by the spatial and seasonal heterogeneity of precipitation has limited its survival and productivity. To facilitate the breeding of drought-tolerant clones and understand the inter-response mechanisms to drought stress, we screened two drought-tolerant (DT) clones and evaluated their differences in physiological and molecular response to drought. The results showed that the No. 228 clone (high-DT ability) had higher antioxidant enzyme abilities than the No. 026 clone (low-DT ability) under drought stress, e.g., peroxidase (POD), polyphenol oxidase (PPO), superoxide dismutase (SOD), and catalase (CAT). Transcriptome analyses revealed that 6637 genes and 1168 genes were up-regulated in No. 228 and No. 026 under drought stress, respectively, when compared to the control (CK). The genes may participate in response to drought-stimulated signal transduction, water/oxygen-containing compound synthesis, photosynthesis, and transmembrane transport functions. Particularly, under drought stress, 14,213 up-regulated and differentially expressed genes (DEGs) were observed in the No. 228 clone compared with the No. 026 clone, and 4274 up-regulated genes were differentially expressed (15-fold difference). These significant DEGs were involved in plant hormone signal transduction, flavonoid biosynthesis, peroxisomes, and other key pathways related to drought. Interestingly, under drought stress, two Chitinases (ClCHIs) and four POD genes (ClPERs) were induced to express in No. 228, which was consistent with the higher antioxidant enzyme activities in No. 228. A heat map of 49 DEGs revealed that dehydrin family genes, ion binding/transmembrane proteins, auxin receptor proteins, and ethylene-responsive transcription factors were significantly up-regulated under drought stress. The results can enhance our understanding of drought tolerance mechanisms and provide a guideline for screening DT genes and breeding drought-tolerant Chinese fir clones.

1. Introduction

Water is a key abiotic factor affecting the survival, growth, and production of plants. Global warming has increased the frequency of extreme weather events, which leads to increased drought stress within plants due to increasing spatiotemporal heterogeneity of precipitation globally and regionally [1]. Drought stress severely hampers development and productivity of agriculture and forest ecosystems [2]. Many recent studies have focused on physiological and biochemical responses to drought stress in plants [3,4]. However, few studies have explored how plants respond to drought stress at the cellular level.

Plants can maintain normal physiological activity by regulating the activities of thousands of genes and metabolic pathways to reduce or repair drought stress damage [4]. Plants would produce reactive oxygen species (ROS) when exposed to drought, and the ROS accumulation can injure biological molecules and cellular organelles and ultimately lead to cell death [5]. To protect cells against the negative effects of excessive ROS, plants have evolved enzymatic antioxidant defense systems (such as superoxide dismutase (SOD), Peroxidase (POD), polyphenol oxidase (PPO), catalase (CAT), etc.)) [6]. Antioxidant enzyme activities jointly establish ROS-scavenging pathways to maintain Hydrogen peroxide (H₂O₂) homeostasis in plants [7]. Some studies have reported that drought tolerance can be enhanced by eliminating ROS accumulation in cells through ROS-mediated signal transduction in plants [8,9]. However, the drought response mechanisms of plants are complicated, and drought tolerance ability varies among plants. More investigations are required to find and characterize the genes, which regulate antioxidant defense systems using transcriptomic approaches.

Some researchers have reported molecular response mechanisms of plants under drought stress. Some genes encoding antioxidant enzymes have been observed, and these genes have been exploited to improve the plants’ drought resistance through transgenic technologies [4,10]. There are eight Ascorbate peroxidase (APX) gene family members in Arabidopsis. AtAPX3 and AtAPX5 are up-regulated and expressed markedly under drought stress, and their mutants are more sensitive to drought [11]. Transgenic lines harboring the AgAPX1 gene exhibit significant increases in APX activity under drought stress, and net photosynthetic rates and relative leaf water contents are enhanced [12]. Overexpression of glutathione peroxidases 5 (GPX5) can decrease the water loss rate in Salvia miltiorrhiza seedlings caused by drought and simultaneously alter the expression levels of transcription factor genes and ROS and ABA pathway genes [13]. AtRCI3, which encodes an active-cationic peroxidase, has been reported to positively enhance tolerance to drought stress in Arabidopsis [14].

Numerous Cis-regulatory elements that respond to drought have been found in the promoters of nine SlSOD genes, and most of them were altered under drought stress in Solanum lycopersicum L. [15]. PtFSD, which belongs to the Populus tremula SOD family, has been reported to be induced and to sustain plant photosynthesis under drought stress [16]. Compared with wild-type leaves, the OsMn-SOD transgenic leaves exhibited lower electrolyte extravasation and higher net photosynthetic rate after drought stress, indicating that OsMn-SOD expression could enhance drought tolerance in rice [17]. In addition, transgenic rice lines harboring OsPPO have been reported to exhibit increased drought tolerance because of higher water potential in their shoots, lower oxidative damage under drought stress, and superior redox balance [18].

Chinese fir (Cunninghamia lanceolata [Lamb.] Hook) is one of the most important evergreen conifer species in subtropical China, with fast growth and high-quality commercial characteristics [19]. It is distributed across 17 provinces and autonomous regions, ranging from 21°31′ to 34°03′ N and 101°30′ to 121°53′ E. The species is used extensively in building, furniture, and bridge construction. It also has numerous medicinal applications; for instance, essential oil from the plant is used to treat pain, wounds, and rheumatism [20]. In recent years, Chinese fir has suffered drought stress caused by increasing spatial and seasonal heterogeneity of precipitation due to global climate warming, especially associated with the subtropical high-pressure area in the Pacific Ocean [21]. Drought stress often occurs in summer and autumn in subtropical China, when Chinese fir is in a rapid growth stage, which requires relatively high water resources [22]. Drought stress has severely impaired the survival and growth of Chinese fir [23]. Previous studies have largely focused on the physiological responses of Chinese fir plantlets to drought stress [24,25]. Transcriptome studies on Chinese fir are mainly focused on lignin and cellulose synthesis, xylem layer formation regulation, and response to low phosphorus stress [23,26,27]. However, our understanding of transcriptomic profiles of tissues in Chinese fir under drought stress remains poor [28]. Consequently, it is essential to investigate, which genes may participate in drought responses for Chinese fir and whether the molecular mechanisms response to drought is different among different drought tolerant clones. These findingss would facilitate marker-assisted selection and the development of drought-resistant cultivars in the near future.

Therefore, we first proposed screening out two drought tolerant (DT) clones, i.e., high DT capacity and low DT capacity using the comprehensive subordinate function method. Subsequently, we investigated the response of antioxidant enzyme activity (POD, SOD, PPO, and CAT) to different drought stress treatments. Finally, leaf samples from different DT clones under moderate drought stress and normal water conditions were selected for transcriptome sequencing to identify drought response genes. To the best of our knowledge, this is the first study to conduct a transcriptomic analysis of Chinese fir based on antioxidant enzyme activity under drought stress, which may serve as a gene expression profile in conifers under drought stress.

2. Materials and Methods

2.1. Screening of Study Material with Different Drought-Tolerance Capacities

2.1.1. Screening Experimental Design

The planting materials used in the experiment were seedlings of 50 Chinese fir clones from the Youxi State Forest Farm in Fujian Province, China. The seedlings were transferred and cultivated in the Intensive Breeding Greenhouse of Fujian Agriculture and Forestry University (FAFU). Healthy and straight seedlings were selected with similar initial heights. To standardize the planting material, six individuals of each clone were selected, and their roots were rinsed with water and then rinsed again three times with distilled water. The plantlets were grown in a liquid culture medium with a ½-strength Hoagland nutrient solution for 5 d under controlled conditions: 16-h/8-h light:dark regime, 25–30 °C, and 50%–60% relative humidity. Based on our previous results and other studies [25,29], 10%~20% PEG addition was suitable for dehydration stress for plants in the seedling stage. Therefore, after acclimatization, 15% polyethylene glycol (PEG) 6000 + ½-strength Hoagland nutrient solution was used to simulate drought stress (PEG 6000, 25322-68-3, Solarbio, Beijing, China), and PEG was not added in the CK treatments. Each treatment had three replicates, and ventilation was performed at 8:00, 12:00, and 18:00 every day for 1 h. After 72 h of simulated drought stress in PEG, leaves from the same parts of seedlings (from the top of the plant to the second branch level) were sampled for analyses to compare the physiological response among high DT and low DT clones [30].

2.1.2. Determination of Leaf Relative Water Content, Plasma Membrane Permeability, and Lipid Peroxidation

Leaf relative water content (RWC) was measured according to the saturation weighing method [31]. Fresh, healthy, and undamaged leaves of the plants were cut, and the fresh weights (FW; M0) were determined. Subsequently, the leaves were soaked in ultrapure water for 24 h, and saturated FW (M1) was weighed again. Finally, the leaves were oven-dried at 80 °C to constant weights (DW: M2). The RWC (%) was calculated as follows: RWC = ([M0−M2]/[M1−M2]) × 100%, where M0 is leaf fresh weight, M1 is leaf turgid weight, and M2 is leaf dry weight.

Leaf plasma membrane permeability was assessed using leaf relative conductivity (REC) [32]. Approximately 0.5 g of fresh leaf samples were initially homogenized in 25 mL of deionized water. Subsequently, the leaves were soaked for an hour at room temperature (25 °C), and the initial conductivity (S1) was determined using a conductivity meter (Starter 3100C/B, Ohaus, Parsippany, NJ, USA). The mixture was heated at 100 °C for 20 min and then cooled at room temperature. The final conductivity (S2) was measured using the procedure described above. The conductivity of the deionized water (S0) was also determined, and the leaf REC was calculated as follows: REC = ([S1−S0]/[S2−S0]) × 100%, where S0 is the conductivity of deionized water, S1 is initial conductivity, and S2 is final conductivity.

Lipid peroxidation was estimated by measuring malondialdehyde (MDA, a product of lipid peroxidation) concentrations in leaves using the thiobarbituric acid (TBA) method [33]. Fresh leaf samples (0.2 g) were homogenized in 4 mL of 0.5 mol L⁻¹ phosphate-buffered solution (PBS, pH = 7; 1% polyvinylpyrrolidone, PVP), and then the mixture was centrifuged at 10,000 rpm for 25 min (Eppendorf 5810R). Afterward, 1 mL of 0.67% (w/v) TBA in 5% (w/v) trichloroacetic acid was added to 3 mL of the supernatant. The mixture was heated at 100 °C for 30 min and then quickly cooled in an ice bath. After cooling, the mixture was centrifuged at 4000 rpm for 10 min, and the absorbances of the supernatants were recorded at 450 nm, 532 nm, and 600 nm. MDA concentration was calculated from the following equation:

$$M=\left[6.45\times \left(A532-A600\right)-0.56A450\right]\times \frac{V}{W}$$

(1)

where M is MDA concentration, A is absorbance, V is the total volume of extract, and W is the fresh weight of samples. Results were expressed as μmol g⁻¹ FW.

2.1.3. Data Analysis

One way-Analysis of Variance (ANOVA) was used to test for differences among clones with the least significant difference (LSD) post-hoc tests (p < 0.05). Data are presented as mean ± standard error. All statistical analyses were conducted using IBM SPSS 25 (IBM Corp., Armonk, NY, USA). To comprehensively evaluate the drought tolerance capacity of each clone, the membership function method in fuzzy mathematics was used [34]. The following equation was used:

$$U\left(X_j\right)=\frac{X_j-X_{min}}{X_{max}-X_{min}}(j=\mathrm{1,2},3\dots n)$$

(2)

where U (X_j) is the membership function value, X_j is the value of the j index, and X_max and X_min are the maximum and minimum values of the j index, respectively.

2.2. Physiological Response of Different DT Clones to Drought Stress

2.2.1. Drought Stress Experimental Design

Based on the screening experiments, two DT types (No. 228 with high drought tolerance and No. 026 with low drought tolerance) were selected as the study materials. The drought stress treatments were pot-water controlled artificially. It included four levels of drought stress: normal water (N, 75%–80% of maximum field water capacity [FWC]), light drought stress (L, 55%–60% FWC), moderate drought stress (M, 45%–50% FWC), and severe drought stress (S, 35%–40% FWC) [35]. The soil moisture content was measured by weighing method at 08:00 every morning [17]. These two DT clones were asexually reproduced, and healthy, uniform seedlings (height: 65.0 ± 5.0 cm; basal diameter: 1.10 ± 0.1 cm) were selected for planting in pots. The cultivation substrate was sand and nutrient soil (1:1). The soil bulk density was 1.74 g cm⁻³, and the maximum FWC was 15.63%. One seedling was planted in each pot, and each treatment had three replicates. Leaf samples for use in the analysis of antioxidant activity (POD, PPO, SOD, and CAT) were collected after drought tolerance experiments for 20 days under varying drought levels. The sampling procedures were similar to those applied for screening drought-tolerance in clones, i.e., all leaves were obtained from the top of the plant to the second branch. To maintain similar experimental conditions, such as light, temperature, and humidity, physiological response experiments were carried out in the same greenhouse at FAFU.

2.2.2. Determinations of Antioxidant Enzyme Activity

Enzyme liquid extraction: 0.2 g fresh leaves from different drought treatments and CK were homogenized in liquid nitrogen in centrifuge tubes (5 mL) for 15 min, and 4 mL PBS (0.5 mol L⁻¹, pH = 7) was added after grinding. The homogenates were centrifuged at 10,000× g for 25 min at 4 °C. The supernatants were used to evaluate POD, PPO, SOD, and CAT activity.

POD activity was evaluated according to Elmekca and Terzioglut [36]. The reaction mixture contained 2.9 mL of 0.05 mol L⁻¹ PBS (pH = 6.0), 1 mL of 0.05 mol/L guaiacol (Number: YZ-111510, Solarbio, Beijing, China), 1 mL of 2% H₂O₂, and 0.1 mL enzyme extract. Absorbance was recorded at 470 nm.

PPO activity was determined using the catechol method [37]. The reaction mixture contained 3 mL of 50 mmol L⁻¹ PBS (pH = 6.0), 1 mL of 0.08 mmol L⁻¹ catechol (Sigma-Aldrich, 430749, Burlington, MA, USA), and 0.3 mL of enzyme extract. After the reaction was complete, absorbance was recorded every 1 min, with four absorbance readings in total obtained. One unit of enzyme activity was defined as the amount of enzyme required to cause a rate of change of 0.001 absorption units per min at 420 nm.

Superoxide dismutase activity (SOD) was determined using nitrogen blue tetrazole (NBT) colorimetry [38]. The reaction mixture contained 1.7 mL of 0.05 mmol L⁻¹ PBS (pH = 7.8), 0.3 mL of 750 μmol L⁻¹ NBT, 0.3 mL of 20 μmol L⁻¹ riboflavin, and 0.1 mL of enzyme extract. The reaction mixture was treated with strong light (4000 lux) for 15 min, placed in the dark to terminate the reaction, and then absorbance was recorded at 560 nm.

CAT was evaluated using the ultraviolet absorption method [39]. The reaction mixture contained 3 mL of 50 mmol L⁻¹ PBS (pH = 7.0), 0.9 mL of 0.01 mmol L⁻¹ H₂O₂, and 0.1 mL of enzyme extract. The absorbance of the measured tube is measured once every 1 min at 240 nm, with a total of four measurements. The absorbance of the control tube with ultrapure water was set at zero.

2.2.3. Data Statistics

One way-ANOVA was used to test the differences in antioxidant enzymy activities among different drought treatments with the least significant difference (LSD) tests (p < 0.05). The differences between two DT clones within the same treatment were compared by independent t-tests. Data is average ±SE. All statistical analyses were conducted using IBM SPSS 25 (IBM Corp., Armonk, NY, USA).

2.3. Transcriptome Analysis

2.3.1. RNA Extraction and Inspection

Leaf samples for transcriptome analysis were collected from moderate drought stress and CK at 20 d for two DT clones. The samplings procedures were similar to those applied for screening experiments and physiological response experiments. Leaf samples (approx. 5 g) (12 samples in total) were frozen in liquid nitrogen, wrapped in tin foil, and total RNA was extracted using an RNA extraction kit (DP441, Tiangen, Beijing, China) according to the manufacturer’s instructions. RNA concentrations were measured using a Qubit^® RNA Assay Kit in a Qubit^® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). RNA integrity was assessed using an RNA Nano 6000 Assay Kit in an Agilent Bioanalyzer 2100 system (Agilent Technologies Inc., Santa Clara, CA, USA). The average RNA concentration was 171.2 ng μL⁻¹, and the average RNA integrity number (RIN) value was 8.35.

2.3.2. cDNA Library Preparation and Transcriptome Sequencing

RNA samples were pooled into 1.5 µg RNA for cDNA preparation. Sequencing was performed at Novogene in Beijing, China (http://www.novogene.com/, accessed on 1 October 2019), using a NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, USA) according to the manufacturer’s instructions. Briefly, poly-T oligo-attached magnetic beads were used to purify mRNA from Chinese fir total RNA, which were then cracked into short fragments at high temperatures. First-strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNaseH⁻) using the short fragments as templates. Afterward, second-strand cDNA was synthesized under the action of DNA polymerase I and RNase H. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptors with hairpin loop structures were ligated to prepare for hybridization. After purifying the library fragments using an AMPure XP system (Beckman Coulter, Beverly, USA), 250–300-bp cDNA fragments were selected. The short-fragment cDNA was connected to the sequencing adapter, and was PCR amplified with 3 µL of user enzyme (NEB, USA). Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. After the PCR products were purified using the AMPure XP system, the library quality was evaluated on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). After the library was qualified, the different libraries were pooled according to the effective concentration and the target data volume requirements, and then Illumina HiSeq sequencing was performed.

2.3.3. Analysis of Illumina Transcriptome Sequencing Results

The original sequences (raw reads) obtained were filtered to remove low-quality and unclear base information reads. De novo assembly was accomplished using Trinity (version r20140413p1) with min_kmer_cov set to 2 by default and all other parameters set to default [40]. The assembled transcripts were hierarchically clustered using Corset (https://code.google.com/p/corset-project/, accessed on 1 October 2019) to obtain unigene sequences for subsequent analyses [41]. The unigene sequences were used for blast search and annotation against an NCBI non-redundant protein (NR) database (http://www.ncbi.nlm.nih.gov, accessed on 1 October 2019) and a Swiss-Prot protein database (http://www.ebi.ac.uk/uniprot/, accessed on 1 October 2019) using Diamond (version 0.8.22) (E-value = 10⁻⁵). Diamond v0.8.22 (E-value = 10⁻³) was also used to compare and classify unigene sequences with homologous proteins in the Clusters of Orthologous Groups (COG) database (http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/COG/, accessed on 1 October 2019). The ‘hmmscan’ homology search algorithm in HMMER3 (E-value = 10⁻²) was used to search the established HMM model of the protein domain to annotate the protein families of genes. NCBI blast (version 2.2.28+) (E-value = 10⁻⁵) was used to compare annotations of unigenes and NCBI nucleotide (NT) databases, and Blast2GO (version 2.5) (E-value = 10⁻⁶) was used for functional annotation of the sequences to facilitate assignment of gene ontology (GO) terms (http://www.geneontology.org/, accessed on 1 October 2019) [42]. KAAS vr140224 (E-value = 10⁻¹⁰) was used to complete the comparison and annotation with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/, accessed on 1 October 2019). iTAK sofeware (version 1.2) was used for transcription factor identification and classification. The number of transcription factors of other species was obtained from the plant transcription factor database (http://planttfdb.cbi.pku.edu.cn/, accessed on 1 October 2019).

2.3.4. Bioinformatics for Functional Annotation of Differentially Expressed Genes

Using the transcriptome stitched from Trinity as the reference sequence, RSEM (version 1.2.15) (default parameter settings [mismatch = 0]) was used to map the clean reads of each sample to the reference sequence [43]. Because each sample had three replicates, differential expression analyses of two conditions/groups were performed using DESeq (version 1.10.1). The filtering threshold was padj < 0.05 and |log2FoldChange| > 1. Controlling the false discovery rate using the method with the p-Value of the result [44]. Genes with adjusted p-Values < 0.05 identified by DESeq were classified as differentially expressed genes (DEGs). Graph analysis of DEGs was conducted using TBtools.

2.3.5. GO Functional and KEGG Pathway Enrichment Analysis for DEGs

GO enrichment uses GOseq (version 1.10.0) (Corrected p-Value < 0.05) based on Wallenius non-central hyper-geometric distribution [45]. First, all the DEGs are mapped to each term of the GO database, and the number of genes for each term is calculated. Subsequently, the difference with the entire genome background is determined. KEGG is the main public database about Pathway. We used KOBAS (version 2.0) (Corrected p-Value < 0.05) to determine the pathways in which the DEGs were enriched significantly relative to all the annotated genes.

3. Results

3.1. Screening of Different Drought-Tolerant Chinese Fir Clones

Under normal moisture conditions, leaf RWC among 50 Chinese clones varied from 54.17%–76.50%, while leaf RWC decreased at different degrees under drought stress, ranging from 35.14%–70.03%. Under normal moisture conditions, leaf relative conductivity among the 50 clones was 9.60%–33.87%, while it ranged from 35.58%–92.94% under drought stress. Under normal water conditions, leaf MDA concentrations were 1.97 μmol g⁻¹–9.82 μmol g⁻¹, which increased to 3.48 μmol g⁻¹–26.16 μmol g⁻¹ under drought stress (Table S1). The comprehensive subordinate function values indicated that there were major differences in drought tolerance among the 50 tested clones, with the highest drought tolerance observed in No. 228 and the lowest drought tolerance observed in No. 026 (Table S2).

3.2. Effects of Drought Stress on Antioxidant Enzyme Activity among High and Low DT Clones

Under normal moisture conditions, all the antioxidant enzyme activities of leaves were significantly higher in No. 228 than in No. 026 (p < 0.05) (Figure 1). Under all drought stress levels, the POD, PPO, SOD, and CAT activities of leaves were significantly higher in No. 228 than in No. 026 (p < 0.05). Under light drought stress, the POD and PPO activities of leaves in No. 026 were only 11.70% and 22.22% of those in No. 228 (Figure 1A,B). With drought stress increasing, POD and PPO activities in No. 228 increased from light drought to moderate drought and then decreased under severe drought stress, while these two activities showed an increasing trend in No. 026 (Figure 1A,B). For both two clones, SOD and CAT activities first increased and then decreased with drought stress increasing, SOD activity was highest under light drought stress, and CAT activity was highest under moderate drought stress (Figure 1C,D).

3.3. RNA-Sequencing and Read Assembly

To maximize the number of genes in the transcriptome, a cDNA sample was prepared from an equal mixture of total RNA isolated from leaves for four libraries, corresponding to the moderate drought stress-028 (A), CK-028 (B), moderate drought stress-026 (C), and CK-026 (D) (Table S3), and which were sequenced on an Illumina high-throughput sequencing platform (NEB, Ipswich, MA, USA). A total of 880,448,444 raw reads were obtained from the above four treatments (A, B, C, and D). After stringent quality checks and data cleaning, we obtained 863,141,178 clean reads, which comprised 129.47 Gb of nucleotides. The average Q20, Q30, and GC (guanine + cytosine) percentages were 96.90%, 91.95%, and 44.49%, respectively. We used the transcript sequence spliced by Trinity as the reference sequence, and corset hierarchical clustering obtained 380,835 unigenes with an average length of 1125 bp (

Table 1. Summary for the Chinese fir transcriptome. Note: Q20 and Q30 are the proportions of bases accounting for the total bases with Phred values >20 or 30. GC percentage is the proportion of guanine and cytosine numbers accounting for the total bases.

Items	Values
Total number of raw reads	880,448,444
Total number of clean reads	863,141,178
Total nucleotide length	129.47 Gb
Total number of unigenes	380,835
Mean length of unigenes	1125 bp
Q20 percentage (%)	96.90
Q30 percentage (%)	91.95
GC percentage (%)	44.49

). The length distributions of spliced transcripts and unigene sequences are illustrated in Figure S1.

3.4. Gene Annotation and Functional Classification

To validate and annotate the assembled unigenes, we searched similar sequences against the NR, NT, protein family (Pfam), euKaryotic orthologous groups (KOG), Swiss-Prot, KEGG, and GO databases using various software and parameters. According to the results, among 380,335 unigenes, 188,114 (49.39%) had significant similarities to known proteins in the NR database. In the NT database, 73,490 (19.29%) DNA sequences were matched. In the Swiss-Prot database, 153,508 (40.3%) DNA sequences matched significantly with known proteins, and 163,178 (42.84%) were matched in the Pfam database. About 164,930 (43.30%) and 49,395 (12.97%) unigenes were matched in the GO and KOG databases, respectively (

Table 2. Annotation of unigene sequences in Chinese fir among different databases.

Sequence Database	Number of Annotated Unigenes	Percentage of Annotated Unigene (%)
NR	188,114	49.39
NT	73,490	19.29
KO	70,927	18.62
Swiss-Prot	153,508	40.30
PFAM	163,178	42.84
GO	164,930	43.3
KOG	49,395	12.97
Annotated in all Databases	19,497	5.11
Annotated in at least one Database	227,019	59.61
Total unigenes	380,835	100

NR: NCBI non-redundant protein sequences; NT: NCBI nucleotide sequences; K: KEGG Orthology; Swiss-Prot: a manually annotated and reviewed protein sequence database; PFAM: Protein family; GO: gene ontology; KOG: euKaryotic Ortholog Groups.

). Because Chinese fir genome sequence information has not been released, 50.61% of the unigenes could not be matched to known genes. Similarly, up to 215,905 unigenes (56.69%) had no GO annotations (Table 2).

To further evaluate the completeness and validity of our transcriptome data, we searched for genes involved in KOG classification in annotated sequences. According to the results, 49,395 unigenes could be annotated into the 25 KOG categories. Most unigenes were categorized into “posttranslational modification, protein turnover, chambers” (7345, 14.87%), “general function prediction only” (6703, 13.57%), and “translation, ribosomal structure, and biogenesis” (6667, 13.51%). Some processes that are associated with drought stress, such as “signal transmission mechanisms” (3317, 6.72%), “intracellular trafficking, secretion, and vesicular transport” (2904, 5.88%), and “carbohydrate transport and metabolism” (2478, 5.02%), provided critical information, which facilitated the uncovering of the molecular mechanisms of drought stress tolerance in Chinese fir (Figure 2).

GO assignments were used to classify the functions of the predicted Chinese fir genes. Based on sequence homology, 164,930 genes could be categorized into three categories, including biological process, cellular component, and molecular function, and could be assigned into 57 functional groups (Figure S2). In the “biological process” category, the “cellular process” (88,387, 53.59%) was dominant, followed by the “metabolic process” (85,098, 51.60%) and “single-organism process” (69,256, 41.99%). In the “cellular component” category, “cell”, “cell part”, and “organelle”, with percentages of 26.91%, 26.89%, and 18.18%, respectively, were the dominant groups. In the “Molecular function” category, the dominant groups were “binding”, “catalytic activity”, and “transporter activity”, accounting for 57.44%, 45.25%, and 6.67% of the groups, respectively, while “metallochaperone activity”, “nucleotide-binding transcription factor activity”, “molecular function regulation”, “antioxidant activity”, “protein binding”, and “transport activity” groups were less than 5000. The most annotated sequences belonged to the ‘cellular process’ category, which are potential novel genes involved in secondary metabolism pathways in drought-stress tolerance (Figure S2).

The KEGG Pathway database records the networks of molecular interactions at the cell, organism, and ecosystem levels. Based on comparisons against the KEGG database using KAAS with an E cutoff value of <10−10, the 70,927 unigenes had significant matches in the database and were assigned to five major branches and 19 secondary branches, with a total of 132 KEGG pathways (Figure S3). The most abundant pathways were “carbohydrate metabolism” (6670, 9.40%), “biosynthesis of other secondary metabolites” (3710, 5.23%), “metabolism of terpenoids and polyketides” (2826, 3.98%), “signal transduction” (1651, 2.33%), and “membrane transport” (611, 0.86%) (Figure S3).

Plant transcription factors control plant development and environmental resistance by regulating the transcription of downstream genes. In the Chinese fir drought stress transcriptome data, 7151 transcription factors were obtained after prediction and classification, which could be divided into 82 categories (Figure S4). Among them were MYB (647, 9.05%), accounting for the highest number, followed by AP2-EREBP (470, 6.57%), Orphans (416, 5.82%), bHLH (368, 5.15%), C2H2 (343, 4.8%), mTERF (338, 4.73%), C3H (303, 4.24%), HB (267, 3.73%), bZIP (237, 3.32%), NAC (229, 3.2%), WRKY (224, 3.13%) and GRAS (194, 2.71%) (Figure S4). The 12 transcription factors above accounted for a total of 4035 transcription factors (56.43%), and the rest of the 70 categories, including LOB, AUX/IAA, ARF, TUB, and HMG, accounted for 3116 transcription factors (43.57%) (Figure S4). The number of transcription factors observed in the present study is much higher than the numbers reported in other species (Table S4), indicating that the genetic information in Chinese fir is much more complex.

3.5. Differently Expressed Genes (DEGs) under Drought Stress

To investigate the genes involved in drought stress responses in Chinese fir, DEGs were screened and analyzed in the leaf transcriptome under drought stress and normal culture in the No. 228 and No. 026 clones, respectively (A vs. B, C vs. D). Filtered with a padj < 0.05 and |log2FoldChange| > 1, 6637 DEGs were up-regulated, and 5124 DEGs were down-regulated between drought stress and normal cultivation in No. 228 clones (A vs. B) (Figure S5A); 1168 DEGs were up-regulated, and 838 DEGs were down-regulated between drought stress and normal cultivation in the No. 026 clones (B vs. D) (Figure S5B) Furthermore, to study the significantly DEGs between No. 228 and No. 026 clones under drought stress and normal culture, we filtered data based on |log2FoldChange| ≥ 12.5, and a total of 296 DEGs (A vs. C) were found, with 176 DEGs up-regulated and 120 DEGs down-regulated (Figure 3A). Two-hundred and fifty genes were differentially expressed (B vs. D), with 116 DEGs up-regulated and 134 DEGs down-regulated (Figure 3B). Higher numbers of DEGs were up-regulated under drought stress than under normal water stress, indicating that the expression of the genes shifts dramatically between the varying drought stress conditions in Chinese fir.

MYB, WPKY, DREB, and other transcription factors were also significant participants in responses to drought stress. There were 61 transcription factor families and 673 DEGs under drought stress (A vs. C), with 375 DEGs up-regulated and 298 DEGs down-regulated (Figure S6A). Among them, key drought-related transcription factor families such as AP2-EREBP, AUX/IAA, HB, GRAS, MYB, NAC, and WRKY had 53, 5, 18, 8, 12, 23, and 14 up-regulated DEGs, respectively (Figure S6A). Under the normal water treatment (B vs. D), there were 55 transcription factor families and 513 DEGs. Among them, 263 DEGs were up-regulated, and 250 DEGs were down-regulated. In the AP2-EREBP, AUX/IAA, HB, GRAS, MYB, NAC, and WRKY transcription factor families, there were only 17, 3, 15, 6, 8, 10, and 7 DEGs, respectively (Figure S6B). The data indicated considerable differences in transcription factor expression between the two Chinese fir clones. Overall, more transcription factors were up-regulated in the No. 228 clone under drought stress, which indicated that the No. 228 clone had the potential to improve drought tolerance.

3.6. Gene Enrichment Analysis for DEGs among Four Libraries

To investigate the biological function of DEGs, they were analyzed and annotated. Under drought stress (A vs. C), 32 and 9 DEGs associated with “response to water” were up-regulated and down-regulated, respectively; 10 and 3 DEGs associated with “photosynthesis” were up-regulated and down-regulated, respectively; 761 and 650 DEGs associated with “signal transduction” were up-regulated and down-regulated, respectively; 53 and 22 DEGs associated with “response to abiotic stimulus” were up-regulated and down-regulated, respectively; and 30 and 15 DEGs associated with “ion migration and transmembrane transport” were up-regulated and down-regulated, respectively (Figure 4A). Under normal water treatments (B vs. D), 585 and 612 DEGs associated with “signal transduction” were up-regulated and down-regulated, respectively. Most DEGs were associated with “adenyl and pure ribonucleotide binding”, while no DEGs were associated with “water shortage, transmembrane transport, and response to abiotic stimuli” (Figure 4B), which shows that under drought stress, the stimulus signals in No. 228 clone leaves were transduced more rapidly, which promoted more active water and nutrient transport processes, resulting in greater drought tolerance in No. 228 clone than in No. 026 clone.

Based on the different biological pathways, a DEG KEGG pathway enrichment analysis was conducted. Under drought stress (A vs. C), the KEGG pathways with most DEGs were “plant-pathogen interaction” (247), “protein processing in endoplasmic reticulum” (205), “plant hormone signal transduction” (131), and others such as “amino sugar and nucleotide sugar metabolism” (116), “flavonoid biosynthesis” (101), “peroxisome” (96), and “terpenoid backbone biosynthesis” (53), “photosynthesis” (51) (Figure S7A). Under normal culture conditions (B vs. D), the DEG KEGG pathway enrichment trends changed, and even on the pathways also enriched under drought stress, the number of DEGs decreased significantly. For example, the most enriched pathways were “plant-pathogen interaction” (210), “Spliceosome” (177), “Phenylpropanoid biosynthesis” (164), “Starch and sucrose metabolism” (151), and “Basal transcription factors” (36) (Figure S7B). Overall, under drought stress, most of the DEGs were enriched on pathways related to drought, and the number of DEGs enriched in the No. 228 clone was greater, indicating increased physiological activities.

3.7. Key DEGs in Response to Drought Stress in Chinese Fir

All the DEGs and key functions such as “response to water”, “signal transduction”, “response to auxin”, “response to stress”, and “transcription factor”, etc., were ranked by filtering based on |log2FoldChange| ≥ 12 (Figure 5). There were nine stress-related protein kinases, two Chitinases, four Peroxidases, two Glycosyltransferase family members, and one Malate dehydrogenase; the two genes up-regulated in response to water belonged to the Dehydrin family; the 11 highly expressed signal transduction genes belonged to Small GTPase, NB-ARC domain, Carbonic anhydrase, LOB domain-containing protein 41, Calcium-binding protein, RRM superfamily, ABC transporter, GTP-binding protein SEC4, Calcium-dependent protein kinase, and 2 Metal ion binding; among the four genes up-regulated in response to auxin, two were auxin-responsive protein, and two were SAUR family protein. The GO annotations for the 12 up-regulated “response to stress” genes were Histidine kinase complex, Cysteine proteinase inhibitor 6, Zinc finger CCCH domain-containing protein 63, L-ascorbate peroxidase, RNA polymerase ll-associated factor 1, Calcium-binding protein CML, Biofilm development protein YmgB, Peroxisome biogenesis protein 6, Damaged DNA binding, Transmembrane receptor Secretin family, Cellular response to DNA damage stimulus, and Hydrolase activity. In addition, 11 transcription factors were up-regulated, among which five belonged to AP2-EREBP and one each of HB, MYB, NAC, PHD, and WRKY families. The 49 key genes were seldom expressed in No. 026 clone leaves under drought stress but were expressed highly in No. 228 clone leaves, indicating that they participated in the response of No. 228 clones to drought stress (Figure 5).

3.8. ClCHIs and ClPERs Highly Expressed in No. 228 Clone

After differential analysis and annotation of the transcriptome sequencing results, two chitinase genes, ClCHI1 and ClCHI2, and four peroxidase genes, ClPER12-1, ClPER12-2, ClPER2B, and ClPER25, were obtained. Their expression levels under drought stress (A, C) and normal culture (B, D) in No. 228 and No. 026 clones were evaluated. According to the results, ClCHI1 was not expressed in leaves (B) under normal conditions in the No. 228 clone. However, the expression levels of ClCHI1 in No. 228 after drought stress (A) increased sharply, while there was no expression in No. 026 clone (C, D) (Figure 6A); ClCHI2 was not expressed under normal conditions (B, D), and it exhibited high expression in No. 228 (A) (Figure 6B). ClCHI1, ClPER12-1, ClPER12-2, ClPER2B, and ClPER25 were also expressed at low levels in No. 228 clone under normal conditions (B); the expressions of these five genes increased suddenly after drought stress (A), while the expressions could not be detected in No. 026 clone (C,D) (Figure 6C–F), which indicated that No. 228 clone had a strong response to drought stress.

4. Discussion

4.1. Differences Responses of Antioxidant Enzyme among High and Low DT Chinese Fir Clones

We screened out two distinct DT Chinese fir clones, i.e., No. 228 and No. 026, by adjusting physiological activities (RWC, REC, and MDA contents) under drought stress (Tables S1 and S2). The method has also been used to screen for drought-tolerance types in Barley and Fig cultivars [46,47]. Subsequently, in the present study, we observed that POD, PPO, SOD, and CAT activities were higher under drought stress than the control for the No. 228 clone. These four enzyme activities were all significantly higher in No. 228 than in No. 026 clones under different drought levels, which may be the key factors enhancing drought tolerance in No. 228 (Figure 1). Drought stress caused the excessive accumulation of ROS, which severely damaged cells. The improvement of antioxidant enzyme activities was beneficial to scavenge ROS to maintain redox balance [48]. However, the plant antioxidant system is limited. We observed that these four antioxidant enzyme activities in No. 228 under severe drought stress were decreased compared to those under moderate drought stress (Figure 1). This was also confirmed in Quercus variabilis seedlings [49].

4.2. Transcriptome Sequencing and Annotation

Transcriptome analysis is crucial for understanding gene function in plant development and the tolerance to adverse biotic or abiotic factors at the genome-wide scale. In the present study, an Illumina HiSeq sequencing platform was used to establish a large-scale database from leaf samples among two different DT Chinese fir clones. We totally obtained 129.47 Gb of data, and a total of 380,335 unigenes with an average length of 1125 bp were successfully spliced (Table 1). The unigenes are much longer than those reported in other tree species (497 bp, 522 bp) [50,51], which indicates that the transcriptomes obtained from the de novo assembly cover adequate width and depth, which can improve the accuracy significantly. In the NR database, approximately 49.39% of the unigenes (188,114 of sequencing 380,335) had significant homologs, whereas only 36.48%, 36.6%, and 41.29% of the unigenes matches in Cornus officinalis [52], Dioscorea composita [53], and Phenacoccus solenopsis [54]. Based on the annotated data in the NT, KO, KOG, etc., databases, the least number of genes successfully annotated in one database was 227,019 (Table 2). Among the unigene sequences of Chinese fir, 132,324 sequences were longer than 1 kb (Figure S1). Such longer sequences provide a higher hit rate, which is consistent with the findings reported in Salvia miltiorrhiza [55], Carrot [56], and Hulless barley [57]. Such large amounts of unigenes can be broadly attributed to “biological regulation”, “signaling”, “cell wall/membrane”, “transcription factor activity”, and “antioxidant activity”, etc., under GO classification (Figure S2). Most of the representative unigenes were annotated to specific pathways, such as “carbohydrate metabolism”, “lipid metabolism”, “amino acid metabolism”, “transport and catabolism”, and “Energy metabolism” in the KEGG database (Figure S3). Such annotations have also been retrieved in transcriptome sequences of Pinus and Paulownia [58,59], indicating that most of the genes participated in the responses of Chinese fir to drought stress.

4.3. Key Enriched DEG Categories in GO and KEGG Pathways

In the GO functional enrichment, we found DEGs differently expressed genes (DEGs) were enriched in the “binding”, “cellular response to stimulus”, “signaling”, and “response to oxygen-containing compound” categories that are associated with drought (Figure 4). In the KEGG pathway enrichment, DEGs were enriched in the “Plant hormone signal transduction”, “Flavonoid biosynthesis”, and “Peroxisome” categories (Figure S7), which is consistent with the findings of previous studies on Potato [60], Oryza sativa [61], and Zea mays [62]. The DEG changes in GO and KEGG enrichments revealed that Chinese fir would respond to drought stress from physiological response to signal transduction. The different enriched DEG categories were not only attributed to species-specific differences but also to the simulated drought methods. Most results were from the polyethylene glycol (PEG)-induced drought, which is significantly different methodologically from soil water drought stress. PEG can quickly decrease water potential in the solution, resulting in the plant being unable to absorb water and dying in several hours or days, depending on PEG stress degrees [25]. However, water stress in the soil is a slow process for plants, which causes dominant differences in physiological response and molecular regulation compared to quick-devastating drought stress. The DEGs enriched in the GO and KEGG pathways indicated that Chinese fir might enhance signal transduction, plant hormone transport, and kinase activity to increase its tolerance under drought stress.

4.4. Transcription Factors Participate in Plant Responses to Drought Stress

Transcription factors are key participants and regulators of plant drought stress. In the present study, a large number of transcription factors were up-regulated under drought stress. Between No. 228 and No. 026 (A vs. C), 21 DEGs encoding the bZIP TFs were detected, and 10 DEGs exhibited significant up-regulation under drought stress, while 11 DEGs were down-regulated (Figure S6). In a previous study, PtabZIP1 in poplar regulated the development of lateral root primordia and enhanced drought resistance by modulating several flavonoid metabolites [63]. In addition, OsbZIP62 overexpression actively regulates rice drought tolerance by modulating antioxidant capacity by regulating the expression of stress-related genes [31]. Similarly, we observed that 14 and 12 MYB and WRKY transcription factors, respectively, were up-regulated, while 23 and 11, respectively, were down-regulated in No. 228 and No. 026 clones under drought stress (A vs. C) (Figure S6). AtMYB2, AtMYB15, and AtMYB96 expression are induced by drought, which also activates the transcription of genes such as dehydration response genes, AtRD22 and AtRD29B, ABA biosynthesis genes, AtABA1 and AtABA2, which, in turn, positively regulates the drought tolerance of plants [64,65,66]. OsWRKY11 and OsWRKY45 overexpression can induce stomatal closure and the expression of other genes and, in turn, decrease water loss rates in leaves while significantly enhancing drought tolerance in rice [67,68]. DREB belongs to a subgroup of AP2/EREBP transcription factors. After drought stress, PeDREB2a expression in Populus euphratica was rapidly and greatly induced, and PeDREB2a overexpression significantly enhanced the length of the root system and plant height significantly, indicating that PeDREB2a participated in the regulation of drought responses [69]. In addition, the expression levels of genes related to drought stress response, photosynthesis, signal transduction, carbohydrate metabolism, and protein protection in Arabidopsis overexpressing AtDREB1A were significantly increased [70]. In the present study, 68 Chinese fir AP2-EREBP transcription factors were differentially expressed, out of which 53 were significantly up-regulated (A vs. C) (Figure S6). A total of 7151 transcription factors were identified, and only 375 were significantly up-regulated after drought stress treatment (Figures S4 and S6), indicating that the up-regulated transcription factors participated in the enhancement of drought tolerance in Chinese fir.

4.5. Expression of Auxin and ABA Response Genes under Drought Stress

Auxin, a major phytohormone, participates in plant growth and development and enhances plant drought tolerance [28]. Plants respond to drought stress by altering auxin synthesis, transportation, and signal transduction processes to regulate root structure development, stomata closure, and cell wall morphology [71,72]. In the present study, under drought stress, 27 auxin-responsive genes exhibited differential expression, and 17 DEGs were up-regulated, with the differential expression in four genes (2 SAUR family proteins) being significantly different (over 20-fold difference) (Figure 4 and Figure 5). SAURs are key auxin transport carriers. AtSAUR40, AtSAUR41, AtSAUR71, and AtSAUR72 belong to the same subfamily and are expressed specifically during the roots development, hypocotyls, vascular development, and stomatal formation, which enhances stress tolerance [73,74]. This implied that the significant upregulation of SAUR may influence the drought tolerance of Chinese fir through the establishment of leaf auxin concentration gradients and stomata closure. ABA induction, synthesis, transport, and signal transduction are key elements in promoting plant responses to drought stress by regulating the stomatal closure and Ca²⁺ concentrations [75]; GBFs, AZI, and LWTs are induced by ABA participation in drought stress resistance [76,77,78]. However, interestingly, during drought stress, a high number of genes associated with auxin, jasmonic acid, and cytokinin were up-regulated, and only one was involved in ABA synthesis, which could be unique in the drought resistance mechanism in Chinese fir.

4.6. Expression of Gene Encoding Protein Kinase under Drought Stress

Protein kinases are involved in signal sensing, transduction, and expression of genes in plants in response to drought stress [78]. In this study, ClCHI1 and ClCHI2 were only expressed in the No. 228 clone under drought stress and not expressed in the No. 026 clone (Figure 6A,B). Some studies found that after drought stress treatment in tomatoes, the inductive degree of CHIs in Lycopersicon chilense Dun (drought-tolerant) is greater than that in L. esculentum Mill (drought-sensitive) [79]. Similarly, FaChi2 and FaChi4 also showed drastically up-regulated expression after 12 h of drought treatment in strawberries [53]. This indicates that Chitinase genes play a role in plant drought resistance. The plant peroxidase gene family has multiple members, such as APXs and PERs. Under drought stress, AgAPX1 overexpression by Apium graveolens L enhanced the activity of peroxidase, reduced the rate of loss of water content in leaves, and significantly improved its drought resistance in Arabidopsis [12]. And the absence of any one of the eight members of AtAPXs will reduce the drought tolerance of Arabidopsis [26]. We found that ClPER12-2 and ClPER2B were highly expressed in No. 228, and the same situation occurred in ClPER12-1 and ClPER25, which was strongly induced by drought (Figure 6C–F). This finding was consistent with the higher peroxidase activity and leaf water content of No. 228 in the physiological response to drought. CfCHIs and CfPERs participate in the positive regulation of drought tolerance in Chinese fir, which may be key genes regulating drought tolerance of Chinese fir.

5. Conclusions

In addition to physiological and biochemical analyses, using RNA-seq-based transcriptomic technologies, we systematically studied the physiological and molecular responses of Chinese fir to drought stress. We have demonstrated that the responses to acid chemical/water/abiotic stimulus/stress/redox reactions, signal transduction, photosynthesis, ion binding, transmembrane transport, flavonoid biosynthesis, and peroxisome were involved in drought stress for Chinese fir. Six expressed genes associated with antioxidant enzyme systems under drought stress were identified. The high levels of expression of the genes in a drought-tolerant Chinese fir clone under drought stress is consistent with higher antioxidant enzyme activity and could facilitate studies on the underlying molecular mechanisms of drought tolerance and lay a foundation for the breeding of drought-tolerant Chinese fir clones.