Chinese fir was grown in the greenhouse of Fujian Agriculture and Forestry University. The growth condition was maintained at ∼25°C with 16 h light/8 h dark photoperiod. We selected Chinese fir with relatively uniform height (1.2–1.3 m) as experimental material. The NW group without bending treatment was used as control (Figure 1). We also performed artificial bending treatment following previous method (Paux et al., 2005) to induce the generation of CW in Chinese fir. The angle between the 0.5–0.7 m height of the trunk and the vertical direction is 45° (Figure 1). At 10:00 AM on 21 October 2020, SDX tissue samples without bending treatment were collected as NW. The bending treatment of all groups started at 8:00 AM on 22 October 2020. The collection time for 2, 26, and 74 h treatment groups is 10:00 AM on October 22, 23, and 25, respectively (Figure 1).

After the bending treatment, the SDX samples were collected immediately. About 20 cm of the trunk at the bent part of the plant was cut. Then we peeled off the bark and scraped the SDX on the upper (the side far from the ground, marked as OW) and lower sides (the side near the ground, marked as CW) of the bent stem. The materials of NW were collected at the equivalent position. All the samples were immersed in liquid nitrogen, and stored at −80°C until RNA isolation.

We collected wood stems from Chinese fir plants that grow straight and lean, respectively. The bark was peeled off, and then the stems were cut into 0.8–1 cm. One-half number of the small fresh stems were cut into 100-micron slices using Leica rotary microtome. The 100-micron slices were dehydrated using alcohol and tert-butanol. Then the water-removed slices were used for scanning electron microscope observation. At the same time, the other small stems were dehydrated and transparentized using alcohol and tert-butanol, and then embedded using paraffin wax. The tissue wax block was cut into 12-micron paraffin slices using Leica rotary microtome. The paraffin slices were dewaxed by xylene and stained using safranin O and fast green. Then the slices were used for upright microscope observation.

The SDX tissue materials were grinded and crushed using the motorized tissue grinders. We extracted the total RNA by using the RNAprep Pure Plant Kit (Polysaccharides & Polyphenolics-rich) (Tiangen, Beijing, China) according to manual. Residual DNA was removed using RNase-free dnase I (Tiangen, Beijing, China). The quantity and purity of RNA samples were evaluated using the NanoDrop 2000 spectrophotometer to meet the requirements of library construction (1.8 < OD260/280 < 2.2). The RNA integrity and precise concentration were detected (using Agilent 2,100 RNA 6000 Nano kit) in Agilent 2,100 Bioanalyzer System. RNA samples with RNA Integrity Number (RIN) ≥ 7.5, RNA concentration ≥300 ng/μL, and total RNA amount ≥2 μg were used in the downstream experiment. In total, 14 RNAs from 7 experimental groups (NW, 2hCW, 2hOW, 26hCW, 26hOW, 74hCW, and 74hOW) with two biological repeats were extracted for RNA sequencing and qPCR validation.

In total, 14 PacBio Iso-Seq libraries were constructed using the SMARTer PCR cDNA Synthesis Kit (Clontech, #634925). Total RNA samples were reverse-transcribed into first-strand cDNA. The distribution characteristics of amplified products were confirmed by electrophoresis and Agilent 2,100. The cycle number of PCR amplification of cDNA was determined by the pre-amplification to obtain a sufficient amount of cDNA. The purified PCR amplicons were used to construct the SMRT sequencing library using the Template Prep Kit (PacBio, #100–259–100). The sequencing library was purified with PB beads and quantified using Qubit, for which quality control was performed using Agilent 2,100. In brief, 14 libraries with 1–10 kb insert sizes were constructed and sequenced by PacBio Sequel II.

In this study, 14 libraries were also constructed for second-generation transcriptome sequencing. SDX mRNA was enriched by Oligo (dT)25 magnetic beads. In first-strand synthesis reaction buffer, with elevated temperature, the rRNA-depleted RNA fragments are cut into smaller fragments using divalent cations. By using random primer and reverse transcriptase, the first-strand cDNA was synthesized. By using buffer, dNTPs (containing dUTP instead of dTTP), RNaseH, and DNA polymeraseⅠ, the double-strand cDNA was synthesized. Then, double-strand cDNA was purified. Next, the purified double-stranded cDNA was performed terminal repair, poly-A tail addition, sequencing adapter connection, and the fragment size selection. The sequencing library was obtained by PCR enrichment. The sequencing library was initially quantified using Qubit, and the effective concentration was accurately quantified by qPCR. In total, 14 strand-specific cDNA libraries were constructed for Paired-End RNA Sequencing with reads length of 150 nt.

We used SMRT Link 9.0 command-line tool isoseq3 program (version 3.3.0) to process subreads (https://github.com/PacificBiosciences/IsoSeq). First, we used the ccs module (using the polish parameter) of the isoseq3 program to generate the circular consensus sequencing (CCS read). Next, we used the lima module of the isoseq3 program to remove the primers of the CCS sequence. We used the refine module to trim poly-A tails and remove concatemer sequences. After this step, we obtained full-length non-concatemer reads (FLNC reads). Finally, we used the cluster module to cluster the FLNC sequences and obtain the redundant isoforms sequence. After the combination of the FASTA files from above 14 processed Iso-Seq libraries, CD-HIT-EST software (Huang et al., 2010) was used to remove redundant sequences using following parameter: -c 0.99 -M 0 -T 0 -G 0 -aL 0.90 -AL 100 -aS 0.99 -AS 30. To improve the quality of isoforms sequences (removed redundant), we used Illumina reads to correct isoforms sequences by LoRDEC software (Salmela and Rivals, 2014) using following parameter: -t 5 -b 200 -e 0.4 -T 40 -k 21 -s 3. Finally, we again used CD-HIT-EST software (-c 0.95 -M 0 -T 0) to cluster isoforms sequences after error correction for redundancy removal. In brief, we obtained non-redundant, non-chimeric, full-length unigenes, which makes subsequent analyses more accurate. The flowchart for the bioinformatics analysis is shown in Figure 1.

For exhaustive annotation of above full-length unigenes, we used BLAST software (version 2.5.0) (Altschul et al., 1997) to search against NR (NCBI Non-redundant Protein Sequence) (Wheeler et al., 2007), NT (Nucleotide Sequence database), and Swiss-Prot (Apweiler et al., 2004), respectively. The NT, NR, and Swiss-Prot annotations correspond to unigenes using blastn, blastx, and blastp, with E-value <1*10^–6. Moreover, the unigenes were also annotated based on the Pfam (protein family) database (Finn et al., 2014) using HMMSCAN module of HMMER software (version 3.3.2) (Potter et al., 2018). Finally, Gene Ontology (GO) (Ashburner et al., 2000) terms were assigned to unigenes by Blast2GO software (Conesa and Gotz, 2008).

ORF of unigenes with a base number >200 nt was predicted using getorf module from EMBOSS software (version6.3.1) (Rice et al., 2000). The unigenes whose ORF amino acids were less than 120 were extracted. We used CNCI (Sun et al., 2013), CPC2 (Kang et al., 2017), and PLEK (Li et al., 2014) to identify lncRNA. The sequences predicted by the three software were intersected to obtain high confidence lncRNAs.

For TF identification, we removed lncRNAs from unigenes firstly. Then the CDS region of unigenes was identified and translated into a protein sequence using ESTScan 3.0, which was provided by the Plant Transcription Factor database (PlantTFDB 5.0) from PlantRegMap (Tian et al., 2020). The TF families were identified by mapping the protein sequence to PlantTFDB 5.0.

Clean data from 14 Illumina RNA-seq libraries were mapped to full-length unigenes using the RSEM software (--bowtie2, version 1.3.0) (Li and Dewey, 2011). We obtained genes read count using rsem-calculate-expression module of RSEM software from the mapping results. The Spearman’s correlation coefficient (SCC) between replicates in each experimental group was calculated as a measure of repeatability. Then, the differentially expressed unigenes among different treatment groups were calculated using DESeq2 (version 1.26.0) (Love et al., 2014) using p-value < 0.05 and fold change>=2.0 as cutoff. According to the classification of lncRNA and TF, the differentially expressed lncRNAs and differentially expressed TFs were obtained.

We further classified differentially expressed unigenes into lncRNAs and coding genes. Pearson correlation coefficient (PCC) between lncRNA and coding gene was calculated and selected to construct the gene co-expression network by Cytoscape software (version 3.8.0) (Kohl et al., 2011) using PCC >0.9 as cutoff.

Using the clusterProfiler package (version 4.0.2) (Wu et al., 2021) with p-value < 0.05 as cutoff, we performed GO functions enrichment analysis on differentially expressed coding genes from the 6 treatment groups and from the co-expression network in 20 GO terms, respectively.

We synthetized the first-strand cDNA using Hifair Ⅲ SuperMix plus, and 1 μg of total RNA was used. The cDNA was used as template, Hieff qPCR SYBR Green Master Mix was used for quantitative reverse transcription-polymerase chain reaction on Mx3005P QPCR system (Stratagene, Santa Clara, CA). Primers for the selected 12 genes (genes related to the lignin and cellulose synthesis pathway, TFs related to the regulation of secondary growth, and predicted lncRNAs) were designed using NCBI Primer-BLAST tool. The primers are listed in Supplementary Table S1. We chose the EF1α (elongation factor1-alpha) gene as an internal reference for normalization. This gene is stably expressed in Chinese fir across different tissues and organs. The final 20 μL reaction contained: primers with final concentration of 0.2 μM, 2 μL template cDNA(1:20 dilution), and 1X Power SYBR Green dye. PCR amplification conditions: 95°C for 5 min, 40 cycles of 95°C for 10 s, 60°C for 30 s. After amplification, a thermal denaturation cycle was performed to determine the dissociation curve for verification of amplification specificity: 95°C for 15 s, 60°C for 15 s, and 95°C for 15 s. All qRT-PCR reactions were repeated 4 times.

In order to observe cell structure of the NW, OW, and CW in xylem of Chinese fir, the paraffin sections from straight and leaning trunks were observed by upright microscope, and the fresh slices after dehydrate were observed by scanning electron microscopy (Figure 2). The width of the entire xylem of the NW is relatively uniform (Figure 2A). The tracheids in the NW were quadrilateral or polygonal and arranged tightly according to the paraffin section under the X200 field of view (Figure 2B), which was further confirmed under X500 scanning electron microscope (Figure 2C). Compared with the straight Chinese fir, the leaning Chinese fir showed obvious eccentric growth, and the xylem of the CW was wider than that of the OW (Figure 2D). Safranin and fast green staining revealed that the color of CW is redder than OW. The lignin content of the xylem in CW is higher than that of OW, which was consistent with the results of earlier studies in Pinus and Chinese fir (Nanayakkara et al., 2009; Li et al., 2021). The tracheid morphology in the OW was consistent with the tracheid morphology in NW (Figure 2E), which was further confirmed under the X500 scanning electron microscope (Figure 2F). The tracheids in the CW were obviously elliptical or round, with intercellular spaces, and the tracheid cell wall also showed a certain degree of thickening (Figure 2G), which were consistent with the results of earlier studies in European yew, red cedar, and Norway spruce (Burgert et al., 2004). The special tracheid morphology in the CW was further confirmed under X500 scanning electron microscope (Figure 2H).

The full-length transcriptome libraries of 14 SDX samples from Chinese fir were sequenced using PacBio Sequel II platform. A total of 382.78 G subreads bases were generated, which included a total of 236,935,007 reads of insert (ROI). The statistics of PacBio raw data are shown in Supplementary Table S2. The average number of ROI generated by each library was 16,923,929 and the average read quality was >99.9%. After using the SMRTlink9.0 command-line tool to process subreads from the ZMW, we obtained a total of 4,136,486 circular consensus sequencing (CCS) reads, and the average number of CCS reads generated by each library was 295,463. The CCS sequence was processed with 5′, 3′ sequencing primer removal and chimeric sequence removal. Finally, a total of 3,682,735 FLNC sequences were obtained. The average FLNC reads number in each library was 263,052. The FLNC reads were clustered by using the cluster module of isoseq3 software resulted in a total of 415,423 high-quality isoform sequences.

We merged 415,423 high-quality isoform sequences from 14 PacBio full-length transcriptome libraries to perform the first de-redundancy for removing redundant sequences with sequence similarity ≥99% using CD-HIT-EST software. In total, there are 139,318 high-quality sequences remaining. Subsequently, the filtered high-quality sequences were further corrected using LoRDEC software and the short-read sequence generated by Illumina RNA-Seq. The corrected sequence was again removed redundant sequences by using the CD-HIT-EST software with a similarity ≥95%. In total, 41,253 high-quality sequences with an average length of 2,245 bp were obtained as full-length unigenes for consequent analysis.

The 41,253 unigenes were searched for NT, NR, and Swiss-Prot databases using BLAST for functional annotation. There were 25,453, 34,823, and 27,207 genes which obtained functional annotations from NT, NR, and Swiss-Prot, respectively. Moreover, 32,134 and 26,767 unigenes were annotated using HMMSCAN and BLAST2GO, respectively. The statistics of the full-length unigenes annotations were listed in Supplementary Table S3. GO classifications of the full-length unigenes were shown in Supplementary Table S4.

In total, there were 41,225 unigenes remaining after removing the unigenes sequences whose base number was less than 200. Subsequently, the ORF of the remaining sequences was predicted using EMBOSS-getorf to exclude the sequences whose ORF is greater than 120. There were 5,527 remaining unigenes for lncRNA identification. In total, there were 4,949, 5,363, and 4,049 lncRNAs based on the result from CNCI, CPC2, and PLEK, respectively. Finally, a total of 3,602 unigenes were annotated as lncRNAs using the common intersection from three software (Supplementary Table S5).

After removing 3,602 lncRNAs from 41,253 unigenes, we annotated TF using 37,651 unigenes according to PlantTFDB database. A total of 1,185 unigenes were annotated as TFs. The most abundant group of TF family was the bHLH family(93), followed by ERF(75), bZIP(74), Trihelix(74), C2H2(73), and C3H(73) family, respectively (Supplementary Table S6).

NGS technology has a greater sequencing depth than third-generation sequencing. Thus, we used the Illumina RNA-seq libraries to quantitatively analyze the expression of the full-length unigenes. The 14 SDX strand-specific mRNA transcriptome libraries were sequenced and generated a total of 324,215,605 paired-end reads. On average, each library produced 23,158,258 reads of 150 nt. The statistics of Illumina RNA-Seq data are shown in Supplementary Table S7. The SCC showed good repeatability and reliability between replicates in each experimental group (Supplementary Figure S1).

According to the alignment of Illumina RNA-Seq, we obtained all the expression levels of unigenes in 14 samples using RSEM software (Supplementary Table S8). The differentially expressed unigenes between the control group (NW) and the six treatment group (2 h_CW, 2 h_OW, 26 h_CW, 26 h_OW, 74 h_CW, 74 h_OW) were identified using DESeq2 R package (p-value < 0.05, log₂FoldChange> 1). From six pairwise comparisons, we obtained the differentially expressed coding genes (DECGs), which included 838 DECGs in all comparisons (Figure 3A). Compared with the NW group, the number of the DECGs detected in the 2 h_CW, 2 h_OW, 26 h_CW, 26 h_OW, 74 h_CW, and 74 h_OW was 4,829, 4,470, 3,211, 3,555, 4,219, 4,447, respectively (Supplementary Table S9). Overall, the number of DECGs in the 2 h treatment group was the largest, followed by the 74 h treatment group, and the least number of DECGs was in the 26 h treatment group. The number of DECGs unique to the six treatment groups was 770, 499, 209, 251, 285, 390, respectively (Figure 3A). The number of unique DECGs in the 2 h treatment group was the largest, and the up-regulated DECGs accounted for about 53% of the total DECGs (Figure 3B). The number of up-regulated DECGs in the 26 h treatment group and 74 h treatment group was significantly reduced, accounting for 24 and 17% of the total DECGs, respectively. For down-regulated DECGs, the number of down-regulated DECGs in the 26 and 74 h treatment groups was higher than that in the 2 h treatment group, and the 74 h treatment group had the largest number of down-regulated DECGs (Figure 3B). In addition, the number of DECGs obtained by CW compared with OW gradually increased in a time-dependend manner upon compression stress (Supplementary Figure S2, Supplementary Table S10).

Among the DECGs, we found a total of 371 differentially expressed TFs (DETFs) (Supplementary Table S11). The percentage of DETFs in DECGs among the six treatment groups was 4.47, 4.70, 2.96, 3.15, 2.75, and 2.74%, respectively (Figure 3C). The proportion of DETFs in CW and OW samples gradually decreases as time increases upon bending treatment (Figure 3C). The 2 h treatment group had the highest proportion of DETFs. Among the DETFs, the number of ERF, bHLH, C2H2, bZIP, MYB, and NAC families was 47, 29, 29, 27, 21, and 20, respectively (Figure 3D).

GO enrichment analysis was performed on the DECGs (p-value < 0.05). The top 20 GO terms with the most significant p-value were retained (Supplementary Table S12). The DECGs of 2 h_CW and 2 h_OW were enriched in GO terms, including cell wall organization, cell wall, xyloglucan metabolic process, xyloglucosyl transferase activity, cell wall biogenesis, etc. (Figure 3E). In addition to the above terms, the 26 h_CW and 26 h_OW treatment groups were enriched in plant-type cell wall, xylan 1,4-β-xylosidase activity, arabinose catabolic process, xylan catabolic process, etc. (Figure 3E). For the 74 h treatment group, only the 74 h_OW treatment group had DECGs enriched in cell wall, xyloglucan metabolic process, cell wall biogenesis, xyloglucosyl transferase activity, etc. (Figure 3E). Xyloglucan is a component of plant hemicellulose (Fry, 1989). Arabinose is often combined with other monosaccharides, and exists in the form of heteropolysaccharides in hemicellulose, pectic acid, the heartwood of coniferous trees, and in some glycosides (Maki-Arvela et al., 2011). According to the GO enrichment analysis of the DECGs in the six treatment groups, we revealed that the DECGs related to the SDX of Chinese fir are involved in the process of plant cell wall production and plant hemicellulose catabolism after the artificial bending treatment. This might be one of the reasons that Chinese fir produces reaction wood tissue to resist bending stress and restore straight growth.

Chinese fir generates CW in response to stress, and the lignin and cellulose content of its xylem have changed. We analyzed the expression of the genes involved in the lignin and cellulose synthesis pathways and the TFs related to the plant secondary growth based on homologous genes from Populus trichocarpa. According to RNA-Seq result, we revealed the differential expression levels of lignin, cellulose-related genes, and secondary growth-related TFs (Figure 4).

Previous studies have found that the phenylalanine biosynthesizes three kinds of lignin-subunits (H- subunit, G-subunit, S-subunit) through a metabolic grid consisting of 11 enzyme families and 24 metabolites (Chen et al., 2015a; Wang et al., 2018; Liu et al., 2021). In gymnosperms, the G lignin monomer accounts for most of its lignin content. In this study, we found that genes of 13 enzymes (involving 9 gene families) were differentially expressed in treatment groups (Figure 4A, Supplementary Table S13). Differentially expressed genes involved in lignin synthesis included four phenylalanine ammonia-lyases genes (PAL1, PAL2, PAL4, PAL5), one cinnamic acid 4-hydroxylase gene (C4H1), one 4-coumaric acid 3-hydroxylase gene (C3H3), one 4-coumaric acid: CoA ligase gene (4CL5), one cinnamyl alcohol dehydrogenase gene (CAD1), one hydroxycinnamoyl-CoA shikimate hydroxycinnamoyl transferase gene (HCT1), two caffeoyl-CoA O -methyltransferases genes (CCoAOMT1, CCoAOMT2), one coniferaldehyde 5-hydroxylase gene (CAld5H1), and one peroxidase gene (PO2), respectively. Most of the genes involved in the lignin monomer synthesis pathway were generally down-regulated first, and then restored to normal expression levels or up-regulated (Figure 4A).

For genes related to the cellulose synthesis pathway (Chen et al., 2015a; Liu et al., 2021), we found 5 cellulose synthase genes (CESA.2, CESA2.1, CESA2.2, CESA4, and CESA8) and two glycosyltransferase family mannan synthase genes (MANS.1 and MANS.2) are differentially expressed upon bending treatment (Figure 4B, Supplementary Table S14). Except for the down-regulation of CESA8 in 26 h_OW, the other genes showed an obvious down-regulation at 2 h stage first, and then recovered to a certain extent at 26 h groups. Finally, their expression levels were up-regulated at 74 h groups.

Two key TFs PtrHSFB3-1 and PtrMYB092 regulate cell wall biosynthesis during the formation of TW in poplar (Liu et al., 2021). PtrMYB092 can directly activate 11 monolignol genes. In the study of the poplar wood formation process, PtrGATA12, as a nuclear localization transcription activator, coordinately regulates the biosynthetic pathway of secondary cell wall components (Ren et al., 2021). In our study, we also detected 11 TFs genes (Figure 4C, Supplementary Table S15), including GATA12, WND6A, WND6B, MYB2, MYB21, MYB092, etc., which are differentially expressed at different treatment groups. The qRT-PCR analysis of these genes confirmed the expression changes detected by RNA-Seq (Figure 4D). The expression analysis of genes related to secondary growth provides valuable clues for explicating the molecular mechanism of CW formation.

Among 9,785 differentially expressed unigenes during the artificial bending process, of which 434 are predicted lncRNAs, and the remaining 9,351 are DECGs (including 371 TFs). The differentially expressed lncRNAs (DElncRNAs) in CW and OW samples at 2, 26, and 74 h treatment groups were shown in Figure 5A (Supplementary Table S16). The expression changes of these lncRNAs were confirmed by qRT-PCR (Figure 5B). And the results of qRT-PCR showed that the expression trends of 3 lncRNAs were consistent with those obtained from RNA-Seq, indicating the accuracy of the RNA-Seq data.

According to the DECGs and DElncRNAs expression, we calculated the PCC (Supplementary Table S17), and the correlation between the expression of DElncRNAs and DECGs was shown through the network diagram (Figure 5C). We revealed a significant correlation (r > 0.9) between 372 DElncRNAs and 6,533 DECGs (6,263 mRNA, 270 TF). GO enrichment analysis was performed on these 6,533 coding genes (p-value < 0.05, keep the 20 GO terms with the most significant p-value, Figure 5D, Supplementary Table S18). Coding genes presenting correlation with lncRNA were enriched in GO terms such as cell wall, xyloglucan metabolism, cell wall organization, cell wall biogenesis, xyloglucosyl transferase activity, plant-type cell wall, etc. (Figure 5D). This may imply that lncRNA is involved in regulating the process of compression stress response in Chinese fir.