Tissue-specific transcriptome analysis of drought stress and rehydration in Trachycarpus fortunei at seedling

Background Trachycarpus fortunei has broad economic benefits and excellent drought resistance; however, its drought response, adaptation, and recovery processes remain unclear. Methodology 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 suppressed, 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. Conclusion 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.


INTRODUCTION
Trachycarpus fortunei (Hook.) H. Wendl. (Fam.: Palmae; Gen.: Trachycarpus) is an evergreen tree that is widely distributed. It is an important economic and landscaping plant. The global climate is warming, the land area affected by water shortage is expanding, the frequency of droughts is increasing, and the severity and duration of droughts are increasing (Stocker et al., 2014;Yao et al., 2015;Lamaoui et al., 2018). The drought stress experiment of 2-year-old T. fortunei from 19 provenances in Guizhou and Yunnan of China found that soluble sugar and Superoxide dismutase (SOD) reached the maximum value at 14 date, due to the difference in the ability of the provenances to resist drought, the lethality of seedlings reached 40% after rehydration of some provenances (Lan et al., 2017). However, there is still a lack of research on the response and adaptation of T. fortunei to drought response.
Typical mechanisms for plants to resist drought include the development of vigorous root systems, the formation of epidermal wax, a thicker epidermal layer, regulation of stomatal conductance to reduce transpiration, elimination of reactive oxygen species, and activation of stress-related hormones. Under the pressure of drought environment, the leaves will close the stomata (Hummel et al., 2010), and the cells will adjust the osmosis to increase the amount of free sugar and free amino acids (especially proline) (Showler, 2013). When it persists, the protein may lose its activity or become denatured, producing excess reactive oxygen species (ROS), stimulating the activity of oxidative stress kinases in leaf cells (Anjum et al., 2011), photosynthesis is inhibited, leading to metabolic dysfunction and destruction of cell structure (Krasensky & Jonak, 2012;Li et al., 2016). ROS are rapidly cleared due to various drought-induced signaling pathways regulated by abscisic acid (ABA) (Contreras-Porcia et al., 2011). Productive ROS signaling requires both rapid increases in ROS, as well as the ability of cells to prevent ROS from reaching damaging levels through synthesis of antioxidants, including flavonols (Postiglione & Muday, 2020). ABA is a predominant hormone that regulates stomatal closure under drought stress. Plants recognize water deficit conditions in roots and that several molecular signals then move from roots to shoots (Takahashi et al., 2020). CLAVATA3/EMBRYO-SURROUNDING REGION-RELATED 25 (CLE25) peptide transmits water-deficiency signals through vascular tissues, and affects ABA biosynthesis and stomatal control of transpiration in association with BARELY ANY MERISTEM (BAM) receptors in leaves (Takahashi et al., 2018). ABA accumulation is enhanced mainly in the vasculature of the leaves, accumulated ABA is supposed to spread from the vasculature to all tissues to mediate stomatal movements and gene expression related to drought stress resistance (Kuromori, Seo & Shinozaki, 2018;Takahashi et al., 2020), drought-responsive genes mainly involved photosynthesis, signal transduction, lipid metabolism, sugar metabolism, wax synthesis, cell wall regulation, osmotic adjustment (Chen et al., 2020). When the damage exceeds a certain level, it will cause plant plastic damage, growth retardation, or even death. So far, little is known about the regulation during the stress response process of T. fortunei, and the dynamic changes of genes when T. fortunei reaches severe drought are still unknown.
The acquisition of full-sibling progeny is to retain the excellent traits of the family through sexual crossing, with the same genetic background, eliminating the influence of different genotype differences within the species. The dynamic changes of physiological indicators are the products of differential gene expression, and the tolerance to drought is a complex phenotypic trait controlled by multiple genes. Transcriptomics analysis is used to identify genes that play an important role in drought tolerance and infer the main mechanisms involved . At present, there is a lack of T. fortunei genome and transcriptome data, and relevant information is not available. In order to understand the molecular resistance and recovery mechanisms of T. fortunei, through the simulation of natural drought stress process at seedling stage, the expression profile changes of related resistance genes were obtained, combined with phenotypes and physiological indicators to evaluate the molecular response of T. fortunei.

Test
In April 2017, we conducted a T. fortunei germplasm survey and selected excellent male and female plants for hybridization, the excellent paternal pollen collected in Guiding, Guizhou province, China, was sexually crossed with the excellent female plants in Huaxi, Guizhou province, China. The full-sibling progeny were obtained. The full-sibling progeny of 1.5-year-old seedlings was selected for transplanting and placed in a greenhouse (Normal growth 4 months after transplanting). The soil type in the pot was humus: yellow soil (1:3). Natural drought stress was simulated and all plants were normally irrigated for 3 days before drought treatment. After the irrigation was stopped, sampling was conducted between 8:00-9:00 in the morning, leaves and root tips were collected every three days from new plants (some leaves were immediately frozen in liquid nitrogen for transcriptome sequencing and RT-qPCR, while others were used for biochemical parameters measurements), phenotypic changes were also recorded. The leaf collection site was unified into mature healthy leaves, root tips were collected from the roots. Samples were collected every 3 days until the fresh leaves of the seedlings were extremely atrophic. A total of 15 days between the cessation of irrigation and the end of the severe drought. At this point, drought stress concluded and the rehydration experiment was conducted, rehydration was performed on the same day to the control (0d) level. After 12 h (Re0.5d), 1 day (Re1d), 3 days (Re3d), and 6 days (Re6d), rehydration and related 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 measure the T. fortunei leaf area.

Determination of biochemical parameters
First, mature leaves and root tips were collected. Using the methods described in The Principle and Technology of Plant Physiology and Biochemistry Experiment Book (Xuekui, 2006) to determine the following indicators. The proline (Pro) content in the leafs were determined according to the method of acid ninhydrin colorimetry (Demirel et al., 2020). The level of Malondialdehyde (MDA) content in the leaves (Ml) and roots (Mr) were determined by Thiobarbituric acid (TBA) method (Liu, Wei & Li, 2014). Leaf superoxide dismutase (SOD) activity was measured spectrophotometrically at 560 nm using the inhibition of photochemical reduction of nitro-blue tetrazolium (NBT) in the presence of riboflavin under light according to Giannopolitis & Ries (1977). Root activity (Rv) in the roots was determined usingtriphenyltetrazolium chloride (TTC) method. Peroxidase (POD) activity was determined using the phenol method. The experiment included 3 biological and technical replicates. The statistics and analysis of physiological related data are perfomed using R software (https://www.r-project.org/), Cluster analysis of the samples was performed using ggtree software package (Yu et al., 2017). A principal components analysis (PCA) was performed using FactoMineR software package (Le, Josse & Husson, 2008).

Transcriptome material collection and library construction
The material selection time nodes included 0 d (start drought), 9 d (ninth day of drought), 15 d (15th day of drought), Re0.5d (twelfth hour after rehydration), and Re1d (the first day/24 h after rehydration), collect the roots tips and mature leaves of the above time nodes, 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. The isolation of total RNA from above samples were isolated according to the instruction manual of the Trizol Reagent (Invitrogen, Carlsbad, CA, USA). Nandrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) was used to determine the concentration and purity of the RNA. The RNA integrity was assessed by agarose gel electrophoresis while its integrity number (RIN) value was measured by Agilent 2100 (Agilent Technologies, Santa Clara, CA, USA). After the quality of the total RNA sample met the requirements (with at least 1µg sample of concentration ≥ 50 ng/mL, OD260/280 = 1.8-2.2), the mRNA was enriched with Oligo (dT) magnetic beads. Further, the mRNA was added with fragmentation buffer and cut into short fragments. Using mRNA as templates, cDNA was reverse transcribed using six-base random primers. The double-stranded cDNA samples were purified, endrepaired, added with poly(A) tails and then ligated to the sequencing adapters to create cDNA libraries. After the libraries passed quality test, they were sent to Shanghai Ouyi Co. (Shanghai, China) for sequencing from both ends with an Illumina HiSeq X Ten machine.

Data processing and analysis
Raw data were quality filtered and reads with adapters or poor-quality sequences were removed using the Trimmomatic tool (Bolger, Lohse & Usadel, 2014). Assembly of the reads was performed using the Trinity (Grabherr et al., 2011) software paired-end to get the assemblage of the transcripts and unigene libraries. Diamond (Buchfink, Xie & Huson, 2015) was used to compare the unigenes to the NCBI NR, EuKaryotic Orthologous Groups (KOG), Gene Ontology (GO), Swiss-Prot, eggNOG, and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Functional analyses were conducted using the HMMER comparison of unigenes to the Pfam database (Mistry et al., 2013). The PlantTFDB database (Jin et al., 2013) was used to identify transcription factors (TFs).
Using Bowtie2 (Langmead & Salzberg, 2012), 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 (Roberts, 2013). First, the correlation heatmap and PCA analysis were performed using the expression levels of all samples, followed by PCA analysis for different tissue samples, roots and leaves. Use DESeq (Anders & Huber, 2012) software to standardize the number of Unigene counts in each sample (basemean value is used to estimate expression), calculate the difference multiple, and use NB (negative binomial distribution test) to test the difference reads significance, and finally screen differentially expressed gene (DEG) according to the difference multiple and the difference significance test results. 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 (Yu et al., 2012). GOSemSim was used to calculate the similarity between GO terms (Yu et al., 2010). GO term similarity clustering was conducted using ggtree. 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 (Ernst & Bar-Joseph, 2006), which was divided into five 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 (Langfelder & Horvath, 2008), 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 (Schmittgen & Livak, 2008). The primer sequences used in this study are provided (Table 1).

Plant-related physiological and biochemical indicators
During the experiment, leaves entered the folded state after 9 d (Fig. 1A), AMC decreased continuously from 0 to 15 d (Fig. 1B). After 9 d, the POD activity reached its maximum value (Fig. 1C), while SOD activity was at its lowest value (Fig. 1D). After 12 d, the leaves were folded and curled, leaf area decreased (Fig. 1E), and the MDA content in the roots increased to its maximum value (Fig. 1F), and the MDA content in the leaves reached its maximum value in Re0.5d (Fig. 1G), but was delayed compared to the roots. The Pro content in the leaves reached its maximum value (Fig. 1H). After 15 d, SOD activity reached its peak value, and Rv (Fig. 1I) and area decreased to their lowest values. Each indicator showed a correlation with the positive/negative modes of stress, Area was positive correlation with Rv (Correlation coefficient = 0.70, P = 1.918869e−05( ÃÃÃ )), Area and Pro show negative correlation (Correlation coefficient = −0.83, P = 1.062141e−08 ( ÃÃÃ )) ( Fig. 1J). Samples of 12d and 15 d were clustered into member 2 (Fig. 1K). After rehydration, most T. fortunei gradually stretched at Re0.5d. PCA of physiological indexes showed that PC1 on the X axis captures the greatest variation and accounts for 50.1% of the total variance around the PCs whilst PC2 on the Y axis captures the second greatest variation and accounts for 25.4% of the total variance around the PCs (Fig. 1L), the PCA results revealed that Re6d and 0 d were closer and coercion was lifted.
Data quality control, splicing, comparison, expression quantification, and differential analysis The raw data of each sample (Q30) ranged from 96.88% to 97.57%. The effective data ranged from 6.98 to 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. 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. 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), this transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under the accession GIYF00000000.
The distance clustering heatmap between the samples revealed that the samples were reproducible and the leaf and root tissues were distinct ( Fig. 2A). PCA revealed that the repeatability was better among groups and the cumulative contribution rate was high, PC1 on the X axis captures the greatest variation and accounts for 29.69% of the total variance around the PCs whilst PC2 on the Y axis captures the second greatest variation  Figure 1 During the experiment, leaves entered the folded state after 9 d (Fig. 1A). and accounts for 9.56% of the total variance around the PCs, samples of leaves and roots are closely grouped together, which is related to tissue specificity (Fig. 2B). PCA analysis of the leaves (Fig. 2C) and roots (Fig. 2D) revealed that the samples of the leaves in each time group were more clearly separated, which may indicate that the leaves responded more to the drought than the roots. After screening for DEGs, compared to 0d_L (Fig. 2E), the DEGs of the corresponding parts of each stage included 4,791 (9 d), 9,108 (15 d), 5,480 (Re0.5d), and 3,211 (Re1d) DEGs. Compared to 0d_R (Fig. 2F), the DEGs of the corresponding parts of each stage included 7,995 (9 d), 8,556 (15 d), 4,968 (Re0.5d), and 6,204 (Re1d) DEGs. 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. 2G), 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 chlorophyl 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 chlorophyl 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.
In the root co-expression analysis network (Fig. 4B), the black module was significantly negatively correlated with area (Correlation coefficient = −0.99, P = 2e−13). Response to water (GO:0009415), response to cold (GO:0009409), and response to salt stress (GO:0009651) (BP category) were enriched (Fig. 4D_1). DHNs COR410-like (TRINITY_ DN51483_c1_g1_i2), LEA14 (TRINITY_DN41505_c0_g2_i2), NAC 47 (TRINITY_ DN48903_c0_g2_i1), and Annexin D2 (TRINITY_DN46177_c1_g1_i1) aggregated in response to water (GO:0009415). The light green module was significantly ; the left of the double axis is the tick mark ("# 6B9EC2") that represents the leaves; the right of the double axis is the tick mark ("# EE861A") that represents the roots; the text includes the number of gene modules, module names, and GO enrichment analysis results (displayed with the smallest p-value for each of BP, CC, and MF category); The gene trends in the module are processed by taking the logarithm of the gene expression in the module, and fitted in a linear simulation.

DISCUSSION
Drought is one of the main abiotic stressors affecting plant growth and development. Drought responses may vary in different plant organs or tissue and the analysis in both mature leaves and root tips is critical for plants to adapt to water deficits (Feng et al., 2017).
After entering the drought state, T. fortunei leaves fold, root activity continuously decreases, Pro activity continuously increases, and the genes related to photosynthesis (ko00195) and photosynthesis-antenna proteins (ko00196) are continuously downregulated (0, −1, −2, −1, 0). Under drought stress, the light absorption of antenna proteins decreases, electron transport rates of PSII and PSI decreases, ROS accumulation increases, photosynthetic pigments are destroyed, RAC activity decreases (Din et al., 2011;Sharma et al., 2012;Guo et al., 2017;Dalal & Tripathy, 2018;Luo et al., 2018;Michaletti et al., 2018). Resuscitation plants are 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 (Dace et al., 1998;Collett et al., 2003). In this study, under drought stress, the expression of chlorophyl 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 thylakoidassociated cellular components proliferated after rehydration (Fig. 4C_3, 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 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 due to dehydration. Plant roots may induce specific stress responses to cope with the early perception of soil water loss (Haas et al., 2019). MDA is a product of membrane lipid peroxidation that is produced when plants are under stress, and the content of MDA can reflect the degree of . The MDA content in the root reaches its peak at 12d, which is faster than the time node of the highest MDA value in the leaf (Figs. 1F and 1G) In the R2 group, phenylpropanoid biosynthesis (ko00940) was curbed, most of the DEGs (48/53) were downregulated in 15d_R, and some genes (27)  . ABA signaling is majorly composed of three core components: ABA receptor, protein kinases, and protein phosphatases (Saini, Singh & Pandey, 2020). PP2C belongs to a group of phosphatases involved in ABA signaling and is a negative regulator (Zhang & Gan, 2012). After ABA, polyethylene glycol and dehydration treatments, the accumulation rate of OsPP2C09 (Oryza sativa L.) transcripts in roots was more rapid and greater than that in shoots (Zhang et al., 2020). Both PP2C-related DEGs in L2 and R2 groups were upregulated after 15 d. In order to absorb water and survive under drought stress, various organic solutes accumulate in the cytoplasm and chloroplasts for osmotic adjustment (Rajasheker et al., 2019). 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 (Wang et al., 2018). The relationship between water transmembrane transporter activity (GO:0005372) and water channel activity (GO:0015250) was positively correlated with area ( Fig. 4D_6, lightgreen), the black module was significantly negatively correlated (Fig. 4B), and the expression of DHNs was generally regulated and induced by ABA, which can reduce root water conductivity (Choi, Zhu & Close, 1999;Hand et al., 2011;Li et al., 2019a). LEA proteins bind to a large number of water molecules and maintain normal metabolism in cells (Ingram & Bartels, 1996). In a previous study, rice OsANN3 was found to mediate Ca 2+ influx by binding to phospholipids, and overexpression significantly increased drought stress survival (Li et al., 2019b). DHNs COR410-like (TRINITY_ DN51483_c1_g1_i2), LEA14 (TRINITY_DN41505_c0_g2_i2), NAC 47 (TRINITY_ DN48903_c0_g2_i1), and Annexin D2 (TRINITY_DN46177_c1_g1_i1) aggregated in response to water (GO:0009415) term which is enriched in root black modules (Fig. 4D_1). qRT-PCR showed Annexin D2 (Fig. 5B), Dehydrin COR410 (Fig. 5D), NAC47 (Fig. 5F) upregulated in roots at 15d. As drought persisted, phenylpropanoid biosynthesis in the roots was suppressed, the DHNs, LEA, Annexin D2, NAC, and other genes were expressed, possibly to protect cell membrane permeability in T. fortunei root tissues.

CONCLUSIONS
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, 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 repressed after rehydration.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
The work was supported by Major Science and Technology Projects in Guizhou Province (Ministry of Science and Technology Major Project [2014] 6024). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.