Comparative transcriptomics analysis of contrasting varieties of Eucalyptus camaldulensis reveals wind resistance genes

Plant material

For the trial, E. camaldulensis was planted in August 2012 in the South China Experiment Nursery, which is located in Suixi County, Guangdong Province. It was built for seed breeding from provenance forests of Australia and India, including 12 populations (provenances), 115 families, and more than 1700 individuals. The trial site and establishment have been described previously by Luo et al. (2014). In this study, the experimental materials were 7-year-old standing trees of two E. camaldulensis genotypes, CA5 (strong wind resistance) and C037 (weak wind resistance), whose wind resistance was measured and estimated based on wind damage caused by several typhoons. In September 2019, wind simulation experiments using tree-pulling tests were conducted in the field. The specific method involved applying the pulling force to the trunk at a tree height of 4 m by using a winch and cable to tilt the standing tree to 30° vertical to the trunk; three trees were selected and tested for each genotype. The immature xylem of the standing tree at breast height was taken at each of two time points, 0 h (before treatment) and 24 h (after treatment), and the wind-resistant CA5 had three biological replicates. Due to a trunk breakage of C037 in the test process, two biological replicates were considered. For simplicity, the samples of two genotypes extracted before and after treatment were named as C037_0 h and C037_24 h, and CA5_0 h and CA5_24 h, respectively, and then used for transcriptome sequencing.

Determination of wood traits

Fibre length and fibre width (FW) were determined by LDA 02 Hi-Res Fiber Quality Analyzer(OPTEST, Canada). Bending strength referred to a method of testing in bending strength of wood (GB/T1936.1-2009 of PRC national standards), bending elastic modulus followed by the method for determination of the modulus of elasticity in static bending of wood (GB/T19362-2009 of PRC national standards). Strength of structural timber parallel to grain followed by the method of testing in shearing strength parallel to the grain of the wood(GB/T 1937-2009 of PRC national standards). Compressive strength parallel to grain was a method used to test the compressive strength parallel to the grain of the wood (GB/T 1935-2009 of PRC national standards). Browning (1967) method was used for the determination of α-cellulose content and holocellulose content, and Huntley et al. (2003) method for the determination of lignin content.

Total RNA extraction and cDNA library construction

Total RNA was extracted using the RNAprep Pure Plant Kit (TIANGEN, Beijing, China), and the concentration and purity of the total RNA samples were detected using a NanoDrop 2000 (Thermo Scientific, Wilmington, DE, USA). The total RNA integrity of 28s/18s and RNA integrity number (RIN) values were detected by Agilent 2100 Bioanalyzer. Next, the mRNA was purified from total RNA with Poly-T oligomagnetic beads. The first-strand cDNA was synthesized with random hexamer primers and M-MuLV reverse transcriptase. The second cDNA strand was synthesised by DNA Polymerase I and RNase H. AMPure XP system (Beckman Coulter, Beverly, MA, USA) was used to purify the library fragment. Phusion High-Fidelity DNA polymerase, Universal PCR primer and Index (X) Primer were used for PCR. Finally, the PCR products were purified by the AMPure XP system, and the library quality was evaluated by Agilent Bioanalyzer 2100 system.

Transcriptome sequencing and functional annotation

High-throughput sequencing was performed on an Illumina HiSeq6000 platform with three technical replicates. The raw reads obtained by transcriptome sequencing were filtered to obtain clean reads by removing reads containing adapters, reads containing poly-N, and low-quality reads from raw data. Meanwhile, Q20, Q30, GC content and sequence repetition level of clean data were calculated. The TopHAT software (Trapnell, Pachter & Salzberg, 2005) was used to compare the sequence of clean reads with the transcriptome sequence of the third generation of E. camaldulensis, and accurate position information of the reference genome was obtained. There were 7 databases for gene annotation: NR (NCBI non-redundant protein sequences), NT (NCBI non-redundant nucleotide), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KO (KEGG homologous database), and GO (Gene Ontology). The personalised analysis platform and tools of the BMK cloud platform were used for in house analysis of transcriptome data and annotation of TFs.

Single-nucleotide polymorphism (SNP) mining

Using software SOAPsnp (http://soap.genomics.org.cn/soapsnp.html) to get all of the unigene sequences SNP detection. The genes of different SNPs in CO37 and CA5 were mapped to the GO annotation file of E. camaldulensis genome, and the Go functional classification statistics and KEGG metabolic pathway analysis were performed for all genes.

Analysis of differentially expressed genes (DEGs)

The CuffQuant and CuffNorm components of the CuffLinks software were used to quantify the transcription and expression levels of genes, and the fragment number per kilobase transcript (FPKM) was used as the index to measure the transcription or gene expression levels (Trapnell et al., 2010). For detecting DEGs, fold change ≥2 and FDR<0.01 were selected as screening criteria (Benjamini & Hochberg, 1995). After the DEGs were screened, BLAST analysis was performed between the DEG and NCBI databases (NR), Swiss-Prot, KEGG, COG, and KOG databases (E- value < 110-5). The DEGs were annotated with reference to the protein amino acid sequence information with high homology to the target gene. Goseq R packages based on Wallenius non- central hypergeometric distribution were used to classify the functions of DEGs in the GO database (Young et al., 2010). The KEGG signal pathway was enriched and analysed using the Kyoto Orthologics-based Annotation System2.0 (KEGG) signal pathway (Mao et al., 2005).

Quantitative Real-time PCR (qRT-PCR) verification

The cDNA of the transcriptome backup sample was used as a template for gene amplification. Twelve differentially expressed genes related to cellulose and lignin were selected to verify the transcriptome data. Primers for qRT-PCR of related genes were designed using Primer 3.0 software, and 18S rRNA was set as an internal reference gene, with three biological and three technical replicates. The primer information is included in Table S1. The cDNA of each sample was diluted and used as a template for the verification of differential genes by qRT-PCR. The qRT-PCR system is presented in Table S2. The qRT-PCR procedure was set according to the instructions corresponding to the fluorescent dye SYBR Green. qRT-PCR was performed using the analytikjena-qTOWER2.2 fluorescent quantitative PCR instrument. The procedure was as follows: 95 °C for 3 min, 95 °C for 10 s, and 58 °C for 30 s, for 39 cycles. The relative expression of genes was calculated using the 2^−ΔΔct method (Livak & Schmittgen, 2001).

Differential analysis of wood quality index between C037 and CA5

The differences in the bending strength, bending elastic modulus, structural timber parallel to the grain, compressive strength parallel to the grain, lignin content cellulose and holocellulose content between CA5 and C037 were very significant (P ≤ 0.01) (Table 1). The differences in the fibre length, fibre width, α-cellulose content between CA5 and C037 were significant (P ≤ 0.05). According to the fibre and physical mechanics analysis, it is found that the fibre length is between 0.53–0.60 mm, the fibre width means was 25.40 µm in the range of 24.00–26.60 mm. Therefore, Nevertheless, there were significant differences in fibre length and fibre width among different families. Moreover, the cellulose and lignin content of CA5 was higher than that of C037. Further analysis results showed that lignin content varied from 19.40% to 26.10%, the average content of Lignin content was 22.32%, the minimum of HC was 76.91%, the maximum value was 82.84%, and the average HC was 79.83%. The range of holocellulose content was 77.00%-82.40%, with an average of 79.90%. The results showed that bending strength ranged from 58.00 MPa to 89.00 MPa with an average of 73.45 MPa, and bending elastic modulus ranged from 5056.00 MPa to 8470.00 MPa with an average of 6962.17MPa. The minimum value of compressive strength parallel to grain was 37.10 MPa, the maximum value was 53.10 MPa and the average value was 45.63 MPa. The minimum value of structural timber parallel to grain was 16.20 MPa, the maximum value was 31.30 MPa, the average value was 24.35 MPa. There are obvious differences between CA5 and C037. The heritability of structural timber parallel to the grain, Compressive strength parallel to the grain, Holocellulose content and α-cellulose content was more than 0.8, indicating that these characters were less affected by environmental factors. The heritability of fibre length and fibre width was less than 0.5, indicating that these two traits are susceptible to environmental factors.

RNA sequencing data quality assessment

In total, 67.84 GB of clean data were obtained from 10 samples based on transcriptional analysis. For each sample, 5.89 GB of clean data were obtained, and the clean read number of the 10 samples was between 19697418 and 28879406. The variation in the number of clean bases was from 5893116176 to 8639188918. The GC content of the obtained data after filtration was between 51.01% and 52.05%. The Q20 and Q30 values were over 97% and over 94%, respectively. The average number of N (in the base) after filtration (Table S3). It indicates that the sequencing quality was high and accurate, and the raw RNA-seq data could be used for subsequent assembly. After a sequence similarity search against eight databases, we annotated 26588 different assembled unigenes (Table S4).

SNP detection

SNP_S were screened from C037 and CA5 samples using SOAPsnp and their locations were determined. Statistics show (Table S5) that the number of SNP_S detected in the samples of non-wind-resistant strain C037 were 87018, 64805 and 57496, and the number of SNPS detected in the samples of wind-resistant strain CA5 was 67361, 63584 and 61089 respectively. Further analysis showed that C/T and A/G had the highest frequency of occurrence, which was more than 25%, while the frequency of other four single nucleotide variations A/C, G/T, C/G and A/T were less than 10%. Among the six variation types, the frequency of C/T was the highest, which may be because the methylated cytosine residues on CpG dinucleotide are easily spontaneously deaminated to form thymine. In order to further understand the information of SNP locus, the genotype of this locus and the mutated genotype of each individual were analyzed. According to the number of reads supporting this locus and the genotype of this locus obtained by GATK3 software, the genotype was different in the samples of wind resistance. According to the statistical results, the heterozygosity of C037 was 40.53%, and that of CA5 was 43.16%.

GO classification of genes where SNP loci reside

In order to understand the function of the screened genes containing SNP, the GO annotation results were further classified. These genes were annotated into 29 functional regions in 3 categories (Fig. 1), including 90 genes in 13 functional regions of Biological process, 93 genes in 8 functional regions of Cellular Component and 57 genes in 8 Molecular functions. In biological processes, the most genes were involved in metabolic process (26), cellular process (25), single organization process (13) and biological regulation (10). Among the cell components, the genes involved in membrane (21), membrane part (20), cell and cell part (15) were the most. Among the molecular functions, the most genes were involved in binding (25) and catalytic activity (21). SNP containing genes was mainly related to the metabolic process of Eucalytus globulus.

KEGG pathway analysis of SNP locus gene

After processing the KEGG annotation results of 80 SNP containing differences unigene in the transcriptome data, it was found that 11 genes with known KEGG function have been annotated into 10 pathways (Fig. 2). Among them, cystaine and methionine metabolism had two unigenes, and the other pathways had one unigene. The results showed that both KEGG metabolic pathway analysis and GO classification were related to metabolism.

Screening for DEGs

The DEGs were screened with a standard fold change (FC) and false discovery rate (FDR) with significant differences across the tested samples (Fig. 3). In total, 7061 DEGs were detected across the four samples, including 4110 upregulated genes and 2951 downregulated genes. We found 1828 DEGs in CA5 under wind stress (CA5_0 h vs. CA5_24h), including 571 upregulated genes and 1257 downregulated genes. However, only 271 DEGs were detected in C037 under wind stress (C037_0 h vs. C037_24h). More DEGs were identified in CA5 under wind stress compared to that in C037. There were 3687 DEGs in C037_0 h vs. CA5_0 h before wind stress, of which 2570 were upregulated and 1117 were downregulated. There were 1275 DEGs in C037_0 h vs. CA5_0 h after wind stress, of which 841 were upregulated and 434 were downregulated. A total of 697 genes overlapped between C037_0 h and CA5_0 h and C037_24 h vs. CA5_24 h, and only 73 genes overlapped between C037_0 h vs. C037_24 h vs. CA5_0 h vs CA5_24 h. The annotated DEGs of the different groups are shown in Table S6. Based on the generated volcano map, we can further visualize the significant difference in the expression levels between the CA5 and C037 samples (Table S7).

In response to wind stress, 762 and 122 DEGs (Table S8) were exclusively detected in CA5 and C037, respectively; this indicates that wind resistance in Eucalyptus depends largely on the differential gene expression. After the clustering analysis of 7601 DEGs detected across these samples, the Clusters of Orthologous Groups of proteins (COG) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) analyses were conducted for the 726 DEGs specific to CA5 in response to the wind stress. We found that cell wall and phenylpropanoid biosynthesis pathway genes potentially play a key role in response to wind stress (Figs. 4; 5). The genes related to cell wall synthesis, especially those involved in the lignin synthesis pathway, could be the key factors for variation in wind resistance across the different Eucalyptus varieties.

GO enrichment analysis

GO enrichment analysis of the candidate wind-resistant DEGs in the four groups revealed that CA5 and C037 were significantly different in their biological processes, cell composition, and molecular function under wind resistance and control conditions. A total of 5,377 GO annotation entries were obtained corresponding to 28,859 genes. Upon further classification and enrichment of the three functional categories, the top five GO classifications with the most significant differences in gene expression at different time points and categories were obtained (Table 2). In terms of cellular components, the integral component of membrane (GO: 0006952) and extracellular region (GO: 0005576) were significantly enriched in all groups, and the integral component of the plasma membrane (GO: 0005887) was significantly enriched in C037_0 h vs. C037_24 h, CA5_0 h vs CA5_24 h, and C037_24 h vs. CA5_24 h. The number of biogenesis genes in CA5 was significantly higher than that in C037 cells. In terms of molecular function, the GO terms corresponding to ADP binding, ATP binding, and iron ion binding were significantly enriched across all the samples, and the protein serine/threonine kinase activity was found significantly enriched in C037_0 h vs. CA 50 h, C037_0 h vs. C037_24 h, and CA5_0 h vs CA5_24 h. In the biological processes category, terms for defence response and protein phosphorylation were significantly enriched across all the samples, and the processes of lignin biosynthesis were significantly enriched in C037_0 h vs. CA5_0 h, C037_24 h vs. CA5_24 h, and CA5_0 h vs. CA5_24 h after wind stress treatment. Compared with C037, CA5 exhibited two more processes of lignin biosynthesis and cell wall biogenesis under wind stress. These findings indicate that the wind resistance trait in wind-resistant genotypes may be related to cell wall biogenesis.

KEGG-based pathways enrichment analysis of the related genes

The KEGG-based analysis further revealed that the number of DEGs in CA5 and C037 varied greatly in different physiological metabolic pathways under control and wind stress conditions (Fig. 6; Table S9). In C037_0 h vs C037_24 h, only pathways associated with protein processing in the endoplasmic reticulum (ko04141) and glycolysis/gluconeogenesis (ko00010) were significantly enriched under wind stress. In CA5_0 h vs CA5_24 h, eight KEGG pathways were significantly enriched under wind resistance. In C037_0 h vs. CA5_0 h, seven KEGG pathways were significantly enriched, and five metabolic pathways were significantly enriched in C037_24 h vs. CA5_24 h. Among them, stilbenoid, diallylheptanoid, and gingerol biosynthesis (ko00945), phenylpropanoid biosynthesis (ko00940), phenylalanine metabolism (ko00360), and flavonoid biosynthesis (ko00941) were involved in C037_0 h vs. CA5_0 h, CA5_0 h vs. CA5_24 h, C037_24 h vs. CA5_24 h. The four pathways may be related to the strong resistance of CA5 to wind stress. Moreover, phenylalanine biosynthesis (ko00940), phenylalanine metabolism (ko00360), and flavonoid biosynthesis (ko00941) pathways are all related to lignin biosynthesis. Therefore, it can be concluded that lignin content and various genes related to lignin synthesis are important factors, which cause variation in the levels of wind resistance in Eucalyptus.

Differential expression analysis of cellulose synthesis genes

GO classification and enrichment analysis further highlighted the potential involvement of DEGs in cell biological process (GO: 0030244), cell microfibril organization (GO: 0010215), cell metabolic process (GO: 0030243), cell catalytic process (GO: 0030245), cell synthesis (UDP forming) activity (GO: 0016760) and cell synthesis activity (GO: 0016759), that is, six pathways related to cellulose synthesis upon wind stress. To better understand the involvement of cellulose synthesis genes in wind resistance in E. camaldulensis, the related DEGs were further analysed. Upon wind stress treatment, 26 genes related to cellulose biosynthesis were differentially expressed (Table 3, In the table, EC_newGene is the abbreviation for Eucalyptus_camaldulensis_newGene.), which were mainly assigned the following GO terms: cellulase synthase (UDP-forming) activity (GO: 0016760), cellulose biosynthetic process (GO: 0030244), and cellular microfibril organisation (GO: 0010215). No significant accumulation of cellulose synthesis genes was observed in C037 before and after wind stress treatment. Sixteen DEGs that were significantly enriched in CA5 were downregulated upon wind stress treatment. In total, 39 DEGs were significantly enriched in C037_0 h vs. CA5_0 h, of which 37 were upregulated and 2 were downregulated. Further, all the 11 DEGs in C037_24 h vs. CA5_24 h were upregulated. EC_newgene_41710, Gene003075, Gene008307, and Gene011799 were upregulated or downregulated in all groups except C037_0 h and C037_24 h. EC_newgene_68496, Gene008307, and Gene011799 were upregulated more than 10 times between C037 and CA5 before and after wind stress treatment. Therefore, Gene008307 and Gene011799 can be used as key candidate genes for wind resistance research in the future.

Differential expression analysis of lignin synthesis Genes

KEGG-based pathways enrichment analysis revealed that some DEGs could be involved in lignin synthesis pathways. For instance, phenylpropanoid biosynthesis (ko00940) was observed Table S10) (Humphreys & Chapple, 2002). A total of 37 DEGs were annotated as phenylalanine ammonia-lyase (PAL), 4-Hydroxycinnamoyl-CoA ligase (4CL), cinnamate 3-hydroxylase(C3H), caffeic acid/5- hydroxyferulic acid O-methyltransferase (COMT), cinnamoyl-CoA reductase (CCR), cinnamyl-alcohol dehydrogenase (CAD), ferulate 5-hydroxylase (F5H), caffeoyl-CoA 3-O-methyltransferase (CCoAOMT), and peroxidase (POX). The expression levels of DEGs across the four samples were analysed (Fig. 7). The results showed that 8 out of 37 genes were present in the 4 groups. In C037_0 h vs. CA5_0 h, the downstream gene006451 (CAD) was down-regulated; however, the other 30 genes were upregulated, including EC_newGene_72280, gene002172, gene007373, gene009919, and gene012455, which were upregulated more than 10 times. In C037_24 h vs. CA5_24 h, all the other 18 key genes in the lignin monomer synthesis pathway were upregulated, except EC_newgene_70770 (POX) and Gene006451 (CAD). In CA5_0 h vs. CA5_24 h, the downstream regulatory genes EC_newgene_11514 (POX), EC_newgene_13932 (POX), and EC_NewGene_14226 (POX) were upregulated, while the other 20 key genes in the lignin monomer synthesis pathway were downregulated. EC_newGene_38808, EC_newGene_42611, EC_newGene_23446, EC_newGene_14226, and Gene012455 were upregulated more than 30 times across C037 and CA5 before and after wind stress treatment. EC_newGene_38808 and EC_newGene_14226, which were upregulated more than 45 times, could be considered key candidate genes involved in response to wind stress.

TFs related to wind stress

Based on differential expression analysis and the available transcriptomic sequence information, TFs were predicted across the two genotypes of E. maldulensis. A total of 78 different TFs, which may be related to the regulation of cellulose and lignin synthesis were obtained (Fig. 8; Table S11). These differentially expressed TFs were primarily classified as 20 MYB families, 15 NAC families, 12 WRKY families, 12 bHLH families, and a small number of LIM, MADS, SBP, C3H, and other TFs. Four TFs were downregulated and 37 TFs were upregulated in C037_0 h vs. CA5_0 h; 17 TFs were upregulated and 4 TFs were downregulated in C037_24 h vs. CA5_24 h, and 11 TFs were up-regulated before and after C037 vs. CA5 wind stress treatment. Sixteen TFs were upregulated and 31 were downregulated after CA5 wind stress treatment, while only three TFs were upregulated after C037 wind stress treatment. This indicated that there was little difference in the expression levels of C037 before and after wind damage. These TFs are involved in the regulation of plant growth and development, morphogenesis, stress resistance, and other biological metabolic pathways.

Homologous evolution analysis was performed between the differentially expressed MYB and NAC TFs as well as the TFs involved in secondary wall synthesis in Arabidopsis thaliana (Figs. 9; 10). EC_newGene_8075 exhibited 93% similarity with ATMYB46, ATMYB6, ATMYB83, EC_newGene_2841, while EC_newGene_35442 showed 99% similarity with ATMYB58 and ATMYB63 (Fig. 7). Gene004569 and gene010621 exhibited a 100% similarity to ATSND1, likewise, EC_newGene_19196 was 100% similar to ATVND1, ATVND2, ATVND3. EC_newGene_5820 and Gene004903 exhibited 100% similarity with ATVND4 and ATVND5, while EC_newGene_27222 and EC_newGene_27223 were 100% similar to ATSND2 and ATSND3.

Verification of sequencing results by qRT-PCR analysis

To verify the accuracy of the transcriptome sequencing results, the expression levels of 12 different genes related to cellulose and lignin synthesis were analysed by conducting qRT-PCR (Fig. 11). The qRT-PCR-based expression trends of these genes were consistent with the results of transcriptome sequencing, and the correlation coefficients ranged from 0.9050 to 0.9997, which are significant and indicate that the transcriptome sequencing data had high repeatability and quasi-certainty.

SNP sites based on RNA-sep

In the natural environment, the mutation frequency of plants is very low. The substances that can induce mutation are usually substances that have certain damage or stimulating effects on plants. Based on the transcriptome data of wind-resistant strain CA5 and non-wind resistant strain C037, SNP loci were screened, which confirmed that there were abundant SNP loci in E. camaldulensis. SNP loci located on genes were annotated. Based on these annotations, the selected wind resistance related genes or SNPs in the pathway may be directly related to wind resistance. However, not all SNPs related to traits are located on genes, and some SNPs located in non-coding regions also play an important role in biological phenotypes.

Relationship between cellulose anabolism and wind stress

The cell wall of plants has a strong filamentous network structure, which can provide mechanical support for cells and the whole plant body; therefore, it can promote the maintenance of stem mechanical strength. The content of the cell wall can reflect the plant anti-lodging ability to a certain extent (Xiang et al., 2010; Hagiwara et al., 1999). In this study, 42 functional annotation genes were differentially expressed and were closely related to cellulose synthesis. The variation trend was consistent with the determination of cellulose and lignin content in the secondary wall, that is, the expression levels in CA5 were much higher than that in C037. This is consistent with the previous results corresponding to the lodging resistance of crops. The increase in cellulose content in wheat straw significantly improves the mechanical strength of the stem and enhances the compressive capacity of the wheat stem (Fan et al., 2012). A significant positive correlation was observed between the breaking resistance of the rice basal stem and the contents of lignin and cellulose (Yang et al., 2009). High cellulose content is beneficial for improving the stem strength of soybeans and enhancing their lodging resistance (Deng et al., 2016). Therefore, the cellulose content in the cell wall and its mechanical properties may have a significant impact on the wind resistance of trees. In this study, EC_newgene_68496, Gene008307, and Gene011799 were upregulated more than 10 times between C037 and CA5 genotypes before and after wind stress treatment, which can be considered as key candidate genes for wind resistance research in the future.

Relationship between lignin anabolism and wind resistance

Lignin can improve the strength of the cell wall and the mechanical strength of the stem. The existence of lignin not only enhances the defence ability of plants against biological stress, but also has great significance in waterlogging resistance, cold resistance, lodging resistance, and other abiotic stresses (Whetten & Sederoff, 1995). The results of a study on lodging resistance of crops have revealed that lignin can improve the cell wall strength and enhance the mechanical strength of the stem, and the content of lignin is closely related to stem rigidity (Turner & Somerville, 1997). Changing the composition of lignin does not affect the growth and development of plants (Vanholme et al., 2012); however, increasing its content can significantly enhance the compressive and lodging resistance in stems. Moreover, the strength of the lodging resistance is proportional to the extent of mechanical strength (Jones, Ennos & Turner, 2001; Baucher et al., 2003). For example, the presence of lignin improves the lodging resistance of wheat, rice, rape, corn, and other important food crops. The lignin content of varieties with strong lodging resistance is significantly higher than that of varieties without lodging resistance (Li et al., 2015; Buranov & Mazza, 2008; Kim & Dale, 2004; Ookawa & Ishihara, 1992). Some studies have shown that the expression levels of PAL, 4CL, C4H, and C3H in the phenylalanine secondary metabolism pathway not only affect the synthesis and content of lignin but also affect the biosynthesis of other secondary metabolites (Baucher et al., 2003). The expression levels of F5H, COMT, and CCOAOMT have a greater impact on the composition of lignin monomers, while CAD and CCR significantly affect the composition of lignin monomers (Boerjan, Ralph & Baucher, 2003). In this study, we identified several DEGs in the lignin monomer synthesis pathway in the phenylpropanoid biosynthetic pathway (ko00940); the expression levels of PAL, C3H, 4CL, COMT, CCoAOMT, CAD, F5H, and POX in CA5 were higher than those in C037, which further indicates that the expression of lignin synthesis genes in wind-resistant strains was higher than that in C037. These DEGs are likely to be involved in the regulation of cellulose and lignin content/component changes in E. camaldulensis, which provides a theoretical basis for screening the major genes affecting wind resistance in E. camaldulensis. Based on their expression profile, EC_newGene_38808 and EC_newGene_14226 can be used as key candidate genes for wind resistance research in the future, as these were found to be upregulated more than 45 times in CA5.

The relationship between TFs and wind resistance

Transcription factors play an important role in the regulatory network of plant stress resistance (Singh et al., 2002). Some transcription factors participate in the regulation of a variety of stress responses. The main transcription factor families involved in plant stress resistance include WRKY family, NAC family, AP2 / ERF family, bZIP family, MYB family, et. Many TFs are involved in the regulation of wind damage in E. camaldulensis. Among the 78 differential TFs that may be related to the regulation of cellulose and lignin synthesis, the MYB and NAC families were the most abundant ones. ATMYB83 and ATMYB46 are nonspecific synthesis activators, which directly target genes of SND1 and node genes regulating secondary wall formation in Arabidopsis thaliana. They regulate lignin synthesis in addition to the overall synthesis of the secondary cell wall including cellulose and hemicellulose. When overexpressed in A. thaliana, it can activate the biosynthetic pathways of cellulose, xylan, and lignin, and activate the expression of promoter genes related to secondary wall synthesis (Zhong, Richardson & Ye, 2007; Mccarthy, Zhong & Ye, 2009). EC_newGene_8075 has 93% similarity with ATMYB83 and ATMYB46; therefore, we can speculate that EC_newGene_8075 may have similar functions. The specific expression of ATMYB58 and ATMYB63 in cells with thickened secondary wall material deposition is regulated by the transcription switch factor ATSND1 and its downstream TF ATMYB46. They activate gene expression by binding to AC elements on promoters of key enzymes in the lignin monomer synthesis pathway, such as 4CL and F5H. Hence, they are important regulators of lignin biosynthesis (Zhou et al., 2009; Hartmann, Sagasser & Mehrtens, 2005). EC_newGene_2841 and EC_newGene_35442 have 99% similarity with ATMYB58 and ATMYB63; therefore, it can be speculated that these two TFs have similar functions. Arabidopsis NAC domain transcription factor SND1 is a key transcription switch that regulates the synthesis of a secondary wall in fibres, which can activate the expression of the ATMYB46 gene, induce the highly upregulated expression of secondary wall-related genes, and promote the deposition of cellulose, xylan, lignin, and secondary wall in cells (Hartmann, Sagasser & Mehrtens, 2005; Zhong et al., 2011). Studies have shown that ATSND2 and ATSND3 are secondary TFs that are directly regulated by ATSND1 (Zhong et al., 2008) and are involved in the regulation of secondary cell wall thickening. It can be inferred that Gene004569, Gene010621, EC_newGene_27222, and EC_newGene_27223 may have similar functions. Among the TFs which are highly similar to MYB and NAC in Arabidopsis, EC_newGene_2841(MYB) and EC_newGene_5820(NAC) were upregulated more than nine times in CA5. The two TFs may account for the significant differences in the wind stress resistance between the two genotypes; therefore, these can be used as candidate genes for further functional verification analysis through overexpression experiments.