Transcriptome Analysis of Drought Stress and Rehydration in the Trachycarpus fortunei Seedling Stage

Background: Trachycarpus fortunei has broad economic benefits and excellent drought resistance; however, its drought response, adaptation, and recovery processes remain unclear. In this study, the response, tolerance, and recovery processes of T. fortunei leaves and roots under drought stress were determined by Illumina sequencing. Results: Under drought stress, T. fortunei reduced its light-capturing ability and composition of its photosynthetic apparatus, thereby reducing photosynthesis to prevent photo-induced chloroplast reactive oxygen damage during dehydration. The phenylpropanoid biosynthesis process in the roots was contained and ABA induced the accumulation of intracellular osmoprotectants, as well as DHNs , LEA, Annexin D2 , NAC , and other genes, which may play important roles in protecting the cell membrane’s permeability in T. fortunei root tissues. During the rehydration phase, fatty acid biosynthesis in T. fortunei roots was repressed. Weighted correlation network analysis (WGCNA) screened modules that were positively or negatively correlated with physiological traits. The Real-time quantitative PCR (RT-qPCR) results indicated the reliability of the transcriptomic data. Conclusions: These findings provide valuable information for identifying important components in the T. fortunei drought signaling network and enhances our understanding of the molecular mechanisms by which T. fortunei responds to drought stress.

, and Rv (Fig. 1I) and area (Fig. 1C) decreased to their lowest values. After 12 and 15 d, samples clustered into one group (Fig. 1K). After rehydration, most T. fortunei gradually stretched at Re0.5d. Each indicator showed a correlation with the positive/negative modes of stress (Fig. 1L). The PCA results revealed that Re6d and 0 d were closer and coercion was lifted (Fig. 1H).
Data quality control, splicing, comparison, expression quantification, and differential analysis The raw data of each sample (Q30) ranged from 96.88-97.57%. The effective data ranged from 6.98-7.18 G. The average GC content was 47.94%. The four base content distribution of each sample was relatively uniform. A total of 66,270 unigenes were spliced with a total length of 55,575,529 bp and an average length of 840 bp. In total, 43,238 (65.25%) unigenes were annotated in the NCBI NR database, while 30,054 (45.35%) unigenes were annotated in the Swiss-Prot database. A total of 51,381 CDS sequences were predicted, of which, 43,373 were predicted by the database comparison method and 8,008 were predicted by ESTScan. Among them, 8.08% of the unigene TFs were compared to oil palm (Elaeis guineensis) and 5.95% were compared to apples (Malus domestica); a total of 15,696 unigene annotations were identified in the TF database and were distributed across 58 families. The raw data were stored in the NCBI/SRA database (BioProject accession No.: PRJNA598974).
After 15 d, the number of DEGs reached its peak value compared to 0 d. The five differential combinations between the L and R sites had a total difference of 11,510 DEGs ( Fig. 2E), which were associated with tissue specificity.

Functional enrichment analysis of DEGs
The GO similarity enrichment (BP category) revealed that the leaf responses to hormones were active in the L1 group (Fig. 2H). In the L2 group, photosynthesis (ko00195), photosynthesis-antenna proteins (ko00196), and porphyrin and chlorophyll metabolism (ko00860) were the most enriched terms (Fig. 2I). Most of the DEGs in the pathway were downregulated in 15d_L. The light-harvesting complex I chlorophyll a/b binding protein 1-5 (Lcha1-5) and Lhcb1-7 were lowered in 15d_L.
The biosynthesis of metabolites in the roots was enriched in the R2 group (Fig. 2I). In the R2 group, the most significantly enriched pathway included phenylpropanoid biosynthesis (ko00940), and most of the ko00940 DEGs were downregulated in 15d_R. In the R1, R3, and R4 groups, fatty acid biosynthesis (ko00061) and fatty acid metabolism (ko01212) were enriched.

Weighted correlation network analysis
The Weighted correlation network analysis (WGCNA) revealed that the leaves divided into 14 modules and the roots divided into 11 modules. In the leaf co-expression analysis network (Fig. 4A), the lightcyan1 and dark turquoise modules were significantly positively correlated with multiple traits, the black module was significantly negatively correlated with multiple traits, and the brown4 module was correlated with rehydration stress expression after 15 d in Re0.5d (Fig. 4C). The GO enrichment analysis (CC category) revealed that it was associated with thylakoid (GO: 0009579) and chloroplast thylakoid able to survive 95% of their cell water loss, one tolerance mechanism is to reversibly shut down photosynthesis, for example, in Xerophyta humilis, psbR, psbA, and psbP were downregulated during dehydration, and complex water regulation expression trends were exhibited [14][15]. In this study, under drought stress, the expression of chlorophyll a-b binding proteins decreased and the synthesis of photosynthesis-related factors decreased.
After rehydration, RAC and other photosynthetic-related genes were activated and recovered to a relatively consistent level after 0 d in Re1d. Gene expression of the thylakoid-associated cellular components proliferated after rehydration (Fig. 4C, brown4).
T. fortunei Psb O, Psb P, Psb Q, Psb R, and other PSII subunits were downregulated under drought stress and gradually increased after rehydration. T. fortunei also reduced its lighttrapping ability and the composition of the photosynthetic apparatus, thereby reducing photosynthesis and increasing drought resistance by leaf folding to prevent light-induced chloroplast ROS damage due to dehydration. accumulated in the cytoplasm and chloroplasts for osmotic adjustment [19]. ABA can induce the accumulation of intracellular osmoprotectants, such as the LEA post-embryonic protein, chaperone proteins, carbohydrates, and Pro, which may be critical for survival under drought stress [20]. The relationship between water transmembrane transporter activity (GO: 0005372) and water channel activity (GO:0015250) was positively correlated with area (Fig. 4D, lightgreen), the black module was significantly negatively correlated ( Fig. 4D), and the expression of DHNs was generally regulated and induced by ABA, which can reduce root water conductivity [21][22][23]. LEA proteins bind to a large number of water molecules and maintain normal metabolism in cells [24]. In a previous study, rice OsANN3 was found to mediate Ca2 + influx by binding to phospholipids, and overexpression fortunei roots is repressed during the rehydration phase after extreme drought.

Conclusion
Under drought stress, T. fortunei reduces its light-trapping ability and the composition of the photosynthetic apparatus, thereby reducing photosynthesis and increasing drought resistance by leaf folding to prevent light-induced chloroplast ROS damage to dehydration.
As drought conditions worsen, phenylpropanoid biosynthesis in the roots is suppressed and ABA induces the accumulation of intracellular osmoprotectants, as well as DHNs, LEA, Annexin D2, NAC, and other genes that may protect the cellular membrane's permeability in T. fortunei root tissues. Fatty acid biosynthesis in T. fortunei roots is also repressed after rehydration.

Test materials
After the investigation and screening of excellent T. fortunei provenances in the early stage, the excellent paternal pollen collected in Guiding, Guizhou province, China, was institutional, national or international guidelines). The whole sibling progeny consisting of 1.5-year-old seedlings was selected for transplanting (four months) and placed in a greenhouse. The soil type in the pot was humus:yellow soil (1:3). Natural drought stress was simulated and plants were normally irrigated for 3 days before drought treatment. At the beginning of the stress experiment, physiological growth, roots, sample retention, and other indicators were measured. The aluminum box soil drying method was used to determine the soil absolute moisture content (AMC). A leaf area scanner was used to sequence the T. fortunei leaf area. Sampling was conducted between 8:00-9:00 AM. After experiment was conducted. After 12 hours (Re0.5d), 1 day (Re1d), 3 days (Re3d), and 6 days (Re6d), rehydration and related indicators were measured.

Determination of physiological indicators
First, mature leaves were collected using plant physiological and biochemical experimental principles and techniques [32]. The proline (Pro) content was determined by acid ninhydrin colorimetry. Malondialdehyde activity (MDA) was determined using the thiobarbituric acid method. Total superoxide dismutase (SOD) was determined using the nitrogen blue tetrazolium photoreduction method. Peroxidase (POD) activity was determined using the phenol method. T. fortunei roots were collected to determine the MDA content and root vitality (Rv) using the TTC method. The experiment included 3 biological and technical replicates.

Transcriptome material collection and library construction
The material selection nodes included 0 d, 9 d, 15 d, Re0.5d, and Re1d, each with three biological replicates. The leaf collection site was unified into mature healthy leaves. Root tips were collected from the roots, washed with distilled water, wrapped in tin foil, and stored in liquid nitrogen. Materials were extracted and purified using the CTAB method.

Data processing and analysis
Raw data were quality filtered and reads with adapters or poor-quality sequences were removed. Assembly of the reads was performed using the Trimmomatic tool [33]. Then, sequence alignment of the clean data and assemblage of the transcripts and unigene libraries were completed. Diamond was used to compare the unigenes to the NCBI NR, KOG (eukaryotic ortholog groups), Gene Ontology (GO), Swiss-Prot, eggNOG, and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases [34]. Functional analyses were conducted using the HMMER comparison of unigenes to the Pfam database [35]. The PlantTFDB database [36]was used to identify transcription factors (TFs).
Using Bowtie2 [37], the number of unigene reads in each sample was obtained. The expression of unigene Fragments Per Kilobase of transcript per Million Mapped reads (FPKM) was calculated using eXpress [38]. A principal components analysis (PCA) was performed using R (https://www.r-project.org/). The differentially expressed gene (DEG) screening threshold was set to p < 0.05 and |foldchange| > 2. The clusterProfiler package was used to perform GO and KEGG enrichment of the DEGs [39]. GOSemSim was used to calculate the similarity between GO terms [40]. GO term similarity clustering was conducted using ggtree [41]. Calculation results were selected if q < 0.01 in order to avoid the very general sets and limit the annotation of > 300 genes and GO terms of < 100 genes. According to the p-value screening, the top 10 KEGG pathways were used to map multiple sets of enriched analysis bubbles. Short Time-series Expression Miner (STEM) software was used to perform the trend cluster analysis [42], which was divided into 5 stages according to the time of occurrence; the combination of different stages of combined genes was input and normalized. The number of models was set to 50.
Screening for the top 5,000 genes of the leaf and root gene expression matrix (background gene set for all genes) was conducted according to the median absolute deviation (MAD) method. A matrix of the relationship between gene expression and sample traits was established using the weighted co-expression network analysis (WGCNA) package [43], which was subsequently transformed into a leader matrix to construct a joint analysis of the modules and traits.
Nine genes were selected to validate the transcriptome data using real-time quantitative PCR (RT-qPCR). SuperMix was removed and cDNA was synthesized using TransScript One- Step gDNA. The first strand of the cDNA fragment was synthesized from total RNA. RT-qPCR was performed on a real-time CFX96 Touch PCR instrument. The PCR reaction conditions were as follows: preheating at 95°C for 30 s, 40 cycles of heat denaturation at 95°C for 5 s, and annealing at 60°C for 34 s. The T. fortunei actin gene was used as the reference gene for data standardization. Each sample was repeated three times and the relative expression levels were calculated using the 2 -ΔΔCt method [44]. The primer sequences used in this study are provided (Table 1).

ABA: abscisic acid
ACCase subunit alpha: acetyl-coenzyme A carboxylase carboxyl transferase subunit alpha, chloroplastic All other data generated or analyzed during this study are included in this published article.

Competing interests
The authors declare that they have no competing interests.

Funding
The work was supported by Major Science and Technology Projects in Guizhou Province X.F wrote the manuscript. All authors have read and approved the manuscript.