The three well-developed D. odorifera trees (about the age of 7 years) were grown at an artificial nursery of Hainan Province in China (19°38′56″N, 110°14′29″E). These three individual trees, which have formed heartwood, were selected as biological replicates without any treatment. These three trees were used for all the experiments. Samples of sapwood (Dosw1, Dosw2, Dosw3) and transition zone (Dotz1, Dotz2, Dotz3) were taken from xylem tissues near the cambium and near the heartwood of each tree, respectively, by following the previously described protocols (Yeh et al., 2020). All the samples were isolated from the trees with the help of a sharp chisel after removing the bark, and the tissues were immediately frozen in liquid nitrogen and stored at −80°C until RNA isolation.

Total RNA from Dotz and Dosw groups of samples was extracted by the CTAB method. Agarose gel electrophoresis (2% agarose gel) was used to evaluate RNA degradation and contamination, and a micro-spectrophotometer was used to assess the purity of RNA (260/280 ratio: 1.9–2.2, 260/230 ratio: ≥ 2.0). RNA integrity was evaluated with the help of RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2,100 system (Agilent Technologies, CA, United States), and samples with an RNA integrity number (RIN) greater than 7.0 were used for further experiments.

Magnetic beads with oligo (dT) were used to enrich mRNA by binding to its polyadenylated tail through A-T complementary pairing. A fragmentation buffer was used for fragmenting the mRNA. cDNA was synthesized using random hexamer primers, buffer, dNTPs, and DNA polymerase I. Double-stranded cDNA was purified by AMPure XP beads, and subjected to end repair, the addition of the poly-A tail, ligation of the sequencing linker, and fragment size selection. Finally, the cDNA libraries were subjected to PCR enrichment and sequenced with the Illumina HiSeq 2,500. After the original data was filtered, the redundancy was removed to obtain clean reads. The HISAT2 software was used to compare the clean reads of RNA-Seq to the reference genome sequence of D. odorifera.

The NEBNext^® Multiplex small RNA library prep set was used to generate small RNA-seq libraries for Illumina^® (NEB, Ipswich, MA, United States.) by following the manufacturer’s instructions. Then, the TruSeq SR Cluster Kit v3-cBot-HS (Illumina) was used to generate a cluster on cBot cluster generation system. Finally, all small RNA-seq libraries were sequenced on an Illumina HiSeq 2,500 platform.

After sequencing, clean reads were obtained by removing the reads containing poly-N, ploy A/T/G/C, with 5′ adapter contaminants, without 3′adapter or the insert tag, and the low-quality reads. In addition, Q20, Q30, and GC-content of the raw reads also were calculated. At the end, all the downstream analyses were performed on sequences ranging of 18–30 bp in length. Bowtie was used to locate/align the small RNAs to the reference genome of D. odorifera and analyze the distribution of small RNAs over the reference genome. The reads that mapped onto the reference genome of D. odorifera were compared with sequences in miRBase (20.0) to obtain the known miRNA; whereas novel miRNAs were predicted with the help of miREvo (Wen et al., 2012) and mirDeep2 (Friedlander et al., 2012). The expression level of miRNA was estimated by the TPM (transcript per million) value. The DESeq R package (3.0.3) was used to perform differential expression analysis of Dotz and Dosw with the |log₂ (fold change)| ≥ 1 and p-value < 0.05 as the threshold.

Clean reads of RNA-seq were obtained by removing low-quality sequencing fragments from the raw reads. Then, the read count of each gene was obtained by mapping the clean reads to the reference genome of D. odorifera with the help of HISAT2. The expression levels of mRNAs were estimated by the calculating FPKM (fragments per kilobase of exon model per million) values. A |log₂ (fold change)| ≥ 1 and q-value < 0.05 were used as criteria to screen differentially expressed genes (DEGs) between Dotz and Dosw samples.

The psRobot is focused on plant small RNA analysis, which has been widely used for target gene predictions (Wu et al., 2012; Zhao et al., 2020; Liao et al., 2021). Therefore, we used its online version (http://omicslab.genetics.ac.cn/psRobot/target_prediction_1.php/) to predict the potential target genes of miRNAs, using the default parameters (Wu et al., 2012). Then, correlation analysis of accumulation of target genes and miRNAs was performed. When the expression levels of target genes had a strong negative correlation of at least −0.8 (p < 0.05) with miRNAs, they were selected for further analysis. All DEGs and the target genes of miRNAs were mapped to individual terms in the GO database (http://www.geneontology.org/) and the number of genes per term was calculated. Then, the GOseq software was used for GO enrichment analysis of DEGs. Analysis of gene regulatory pathways in the KEGG pathway database (http://www.genome.jp/kegg/pathway.html/) was performed with the help of KOBAS (3.0) software. Finally, Cytoscape 3.9.0 was used to construct a co-expression network graph.

cDNAs were synthesized by reverse transcription of total RNA from six D. odorifera samples (Dosw1, Dosw2, Dosw3, Dotz1, Dotz2, and Dotz3). Gene-specific primers for the targets were designed with the help of Primer Premier v5 software (Supplementary Table S1). Four important DEGs (WRKY22, AP2, FPP7, and PAL2), five miRNA target genes (SPL12, ARF8, FH20, GCS1, and MMT1) of miRNA-mRNA correlation network, and six differentially expressed miRNAs (miR156g, miR167e, miR168a-5p, novel_9, novel_15, and novel_52) were chosen to verify expression levels. For the genes, RT-qPCR analysis was conducted with TB Green^® Premix Ex Taq™ (Tli RNaseH Plus; Takara, Beijing, China) following the manufacturer’s recommendations. The amplification was performed on a preheated (94°C) thermal cycler, and samples were incubated at 94°C for 2 min, followed by 40 cycles of 95 °C for 5 s and 60°C for 30 s. The actin gene served as an internal control for normalization (Meng et al., 2019). We used the miRNA RT-qPCR Detection Kit (Aidlab, Beijing, China) for RT-qPCR analysis of miRNAs, following the manufacturer’s recommendations. PCR amplification was performed at a preheated (94°C) thermal cycler and samples were incubated at 94 °C for 2 min, followed by 40 cycles of 94°C for 15 s and 60°C for 40 s. The U6 gene served as an internal control for normalization (Meng et al., 2022). The 2^–△△Ct method was used to calculate the expression levels of the miRNAs and genes against the internal controls (Schmittgen and Livak, 2008). Three technical replicates per sample were analyzed to ensure reproducibility and reliability.

To understand the molecular mechanism of xylem differentiation, six cDNA libraries from the Dotz and Dosw of D. odorifera were sequenced. RNA-seq generated 20,567,345 (Dotz1), 20,373,919 (Dotz2), 25,875,484 (Dotz3), 21,492,446 (Dosw1), 22,895,259 (Dosw2) and 23,017,502 (Dosw3) clean reads (Table 1). Of these, 19,643,871 (Dotz1), 19,430,607 (Dotz2), 24,690,387 (Dotz3), 20,041,706 (Dosw1), 21,867,262 (Dosw2) and 21,684,789 (Dosw3) clean reads, respectively, were mapped to the D. odorifera genome with mapping ratios of 93.25–95.51% (Table 1). These results show that RNA-seq captured a significant portion of the genes in the genome of D. ororifera. For differential analysis of gene expression between the Dotz and Dosw groups, the FPKM values were used to normalize the reads from RNA-seq, and a cutoff of |log₂ (fold change)| ≥ 1 and q-value <0.05 was used. A total of 395 mRNAs were differentially accumulated, of which 72 were up- and 323 were down-regulated, respectively, in Dotz compared to Dosw (Figure 1A; Supplementary Table S2). All 395 differentially expressed mRNAs (Supplementary Table S2) obtained from Dotz and Dosw were used for subsequent analysis.

To annotate and reveal the function of DEGs in different tissues of xylem, we used GO classification. Up-regulated mRNAs were associated with 11 molecular functions, 49 biological processes, and four cellular components (p-value < 0.05), including “anion transport,” “anion transmembrane transporter activity,” and “inorganic anion transport” associated with ion transport (Figure 2A). Down-regulated mRNAs were associated with 41 molecular functions, 54 biological processes, and three cellular components (p-value < 0.05), including “inorganic anion transmembrane transporter activity”, “transmembrane transport”, “sulfate transport” and “sulfate transmembrane transporter activity” terms associated with anion transport (Figure 2B). The DEGs were further referenced through the KEGG database. Eight KEGG pathways (p-value < 0.05) were significantly enriched, including “carotenoid biosynthesis,” “diterpenoid biosynthesis,” and “alpha-linolenic acid metabolism” (Figure 3).

Functional annotation of DEGs revealed the processes, such as the transcriptional regulation, lignin synthesis, flavonoids, and terpenoids synthesis, which could play important roles in the xylem differentiation of D. odorifera. We found that 11 DEGs might encode for transcription factors, of which eight belong to the WRKY transcription factor family (Figure 1B). Four genes (PER47 (Peroxidase 47), COMT (Caffeic acid 3-O-methyltransferase), CCR2 (Cinnamoyl-CoA reductase 2), and UXS2 (UDP-glucuronic acid decarboxylase 2)) are related to lignin synthesis and cell wall composition, two genes (KAO2 (Ent-kaurenoic acid oxidase 2) and GA20OX1 (Gibberellin 20 oxidase 1)) are part of gibberellin synthesis, two genes (OPR2 (12-oxophytodienoate reductase 2) and CYP74A (Allene oxide synthase)) participate in jasmonic acid synthesis, and two genes (PAL2 (Phenylalanine ammonia-lyase 2) and CYP71A1 (Cytochrome P450 71A1)) are related to flavonoids and terpenoids synthesis (Figure 1B).

Small RNA-seq was performed to unveil the possible role of miRNAs in regulating gene expression during xylem differentiation in D. odorifera. 15,809,909 (Dotz1), 15,412,631 (Dotz2), 15,469,161 (Dotz3), 15,422,891 (Dosw1), 15,394,233 (Dosw2) and 15,755,203 (Dosw3) raw reads were obtained from the small RNA libraries generated from the Dotz and Dosw groups. After removing connectors and low-quality reads, a total of 45,496,409 (Dotz) and 45,444,352 (Dosw) clean reads were obtained (Table 2). The length of most clean reads ranged between 21 and 24 nucleotides (Supplementary Figure S1). In comparison by Bowtie software, 90.87–95.30% of small RNA reads matched with the genome sequence of D. odorifera (Table 2).

The 18–30 nt clean reads that mapped to the reference genome were compared to sequences in miRBase to identify conserved/known miRNAs. Subsequently, novel miRNAs were predicted using miREvo and miRDeep2 tools on the basis of structural characteristics of miRNA. We identified 38 known miRNAs in Dotz and 36 in Dosw samples, respectively. Together, these accounted for a total of 40 known miRNAs that belonged to 22 miRNA-families (Figure 4A; Supplementary Table S3; Supplementary Figure S2). Further analysis of the distribution of the first base of these known miRNAs showed that “U” was the most dominant base (Supplementary Table S4; Supplementary Figure S3). Also, a total of 123 novel miRNAs were identified, where 122 were present in Dotz and all 123 in Dosw (Figure 4A; Supplementary Table S5; Supplementary Figure S4). The statistical analysis of the distribution of the first base of these novel miRNAs also found that “U” was the most dominant base (Supplementary Table S6; Supplementary Figure S5).

To understand the miRNA-driven mechanism of regulation of xylem differentiation, miRNAs that accumulated differentially between the two tissue types were identified by the TPM method. Compared with Dosw, four known (miR168a-5p, miR167f-5p, miR167h-5p, and miR167e) and three novel miRNAs (novel_15, novel_9, and novel_165) were up-regulated in Dotz, whereas two known (miR390a, miR156g) and five novel miRNAs (novel_45, novel_52, novel_110, novel_47, and novel_130) were down-regulated (Figure 4B). Among these differentially expressed miRNAs, the most up-regulated miRNAs were miR167e (52.21 fold), miR167h-5p (20.53 fold), and miR167f-5p (16.33 fold). The top down-regulated miRNAs were miR390a (52.75 fold), novel_45 (21.44 fold), and novel_52 (17.17 fold) (Supplementary Table S7).

To understand the potential function of differentially expressed miRNAs, 1,056 candidate target genes were predicted with the help of psRobot. Then, the target genes with a strong negative correlation of at least −0.8 (p < 0.05) with miRNAs were selected for further analysis (Supplementary Table S8). We, thus obtained 21 putative target genes for nine differentially expressed miRNAs. Subsequently, GO enrichment analysis indicated that these 21 target genes are associated with 42 GO terms (p-value < 0.05) from biological processes, molecular functions, and cellular components. Few examples of such processes included “amino acid binding,” “cellulase activity,” and “DNA beta-glucosyltransferase activity” (Supplementary Figure S6A). Further, KEGG enrichment indicated that the target genes were enriched in five pathways (p-value < 0.05), including “protein processing in endoplasmic reticulum,” “fatty acid elongation,” and “biosynthesis of unsaturated fatty acids” (Supplementary Figure S6B).

To explore the relationship between miRNAs and mRNAs during xylem differentiation of D. odorifera, a regulatory network diagram of the nine differential expression miRNAs and their 21 target genes was constructed (Figure 5; Supplementary Table S8). It was evident that one miRNA could target 1-7 mRNAs, among which miR156g and novel_52 could negatively regulate seven and six mRNAs, respectively, while miR167h-5p and miR167f-5p could simultaneously regulate one mRNA (pyrH, Uridylate kinase). Similarly, miR156g might target the transcripts of Squamosa promoter-binding-like protein 6 (SPL6), which is involved in programmed cell death (Wang Q.-L. et al., 2018), as well as the SPL12. A novel miRNA, novel_130, might target mRNA of mannosyl-oligosaccharide glucosidase (GCS1), involved in cellulose synthesis. The up-regulated miRNA novel_9 was predicted to target auxin response factor 8 (ARF8) that is involved in auxin dynamics balance (Tian et al., 2004) (Figure 5).

We determined the relative transcript levels of miRNAs and genes identified from DEGs and the interaction network by RT-qPCR, and observed that only miR168a-5p was different with RNA-seq. The expression patterns of all nine genes and other miRNAs estimated by RT-qPCR and RNA-seq followed a similar trend in Dotz and Dosw samples. Although the fold-change (FC) values calculated by sequencing did not exactly match the expression values obtained by RT-qPCR, the expression profiles were basically consistent for all tested miRNAs and genes (Figure 6). These analyses validate the reliability of the gene expression values generated from sequencing results.

By analyzing the patterns of mRNA accumulations in two tissue types related to wood formation in D. odorifera, we identified 395 differentially expressed mRNAs between Dotz and Dosw (Figure 1A). Further, the enriched GO terms associated with the DEGs were largely expected. Anion transport and inorganic anion transport were highly enriched, indicating the function of these genes in wood formation. KEGG analysis showed that DEGs were commonly enriched in some pathways, including “carotenoid biosynthesis,” “diterpenoid biosynthesis,” and “alpha-linolenic acid metabolism”. This result is also consistent with our expectations as these three pathways might play an important role in the differentiation of xylem, especially in the biosynthesis of terpenoids and plant hormones (Biesgen and Weiler, 1999; Regnault et al., 2014).

We have identified genes that are possibly involved xylem differentiation process. Eleven genes that may encode transcription factors were identified, eight of which belong to the WRKY family. Our findings on differential expression of many genes encoding transcription factors are supported by previous investigations on other trees, such as rubber tree (Meng et al., 2021) and Populus tomentosa (Chen et al., 2015). In P. tomentosa, some genes of WRKY family were also differentially accumulated amoung the xylem tissues of tension wood, opposite wood, and normal wood (Chen et al., 2015). Further, earlier studies have shown that WRKY transcription factors were widely involved in plant responses to biotic, abiotic, and hormonal stresses, and regulate the biosynthesis of phenols, terpenes, and alkaloids (Schluttenhofer and Yuan, 2015). For example, MdWRKY11 promotes the expression of F3H, FLS, DFR, ANS, and 10UFGT in Malus domestica callus, and then increases the accumulation of flavonoids and anthocyanins (Wang N et al., 2018). Further, we also found some genes (such as PER47, COMT, CCR2, UXS2, KAO2, GA20OX1, OPR2, CYP74A, PAL2, and CYP71A1) related to lignin synthesis and cell wall composition, gibberellin synthesis, jasmonic acid synthesis, and flavonoids and terpenoids synthesis. Similarly, in the study of Sun et al. (2020), seven PAL family genes involved in the phenylpropane pathway were also significantly differentially expressed between the discolored wood by pruning and the normal wood in D. odorifera. In addition, COMT silenced transgenic poplar trees have significantly reduced 17% lignin levels (Jouanin et al., 2000). The trunk dry mass of a three-month-old PdGA20ox1 overexpressing transgenic poplar was four times higher than that of untransformed control poplar; the contents of xylose and glucose were also significantly increased in these transgenic plants (Park et al., 2015).

MiRNAs are an important class of non-coding factors that regulate genes expression at the post-transcriptional level in nearly every aspect of plant development such as germination, growth, tissue differentiation, and flowering (Das et al., 2019; Glazińska et al., 2019; Qiu et al., 2019). The results of small RNA-seq demonstrate that the xylem of D. odorifera contains a large and diverse small RNA population, a finding similar to previous studies in rubber tree (Meng et al., 2022) and Chinese fir (Wan et al., 2012).

A total of 40 known miRNAs belonging to 22 miRNA families and 123 novel miRNAs were identified in these small RNA libraries, among which 14 miRNAs were differentially expressed (Figure 4). Similarly, the miR156, miR159, miR166, miR319, miR396, miR398, and miR408 families have also been identified in the xylem of rubber tree (Meng et al., 2022). The miR396a and miR156g are also differentially expressed in the primary stem, transition stem, and secondary stem of Populus trichocarpa (Wang R et al., 2021). Then, we found many miRNA families in D. odorifera, for example, miR156 (Wang et al., 2011), miR166 (Chen et al., 2018), and miR397 (Lu et al., 2013), are known to play a role in wood formation in other species. Subsequently, functional analysis of the target genes of the differentially expressed miRNAs revealed that these genes were commonly enriched in a few pathways including, “protein processing in endoplasmic reticulum,” “fatty acid elongation” and “biosynthesis of unsaturated fatty acids” These results further suggest that these differentially expressed miRNAs possibly regulate xylem differentiation of D. odorifera.

We predicted 1,267 target-miRNA pairs between the 14 differentially expressed miRNAs and 1,056 candidate genes with the help of psRobot in Dotz and Dosw. Among these miRNA target related gene pairs, 22 showed a significantly negative correlation. Some studies have found that miR156 regulates plant growth and development, morphogenesis, anthocyanin accumulation, gibberellin synthesis, and stress response by inhibiting the expression of SPL transcription factors in A. thaliana at the post-transcriptional level (Schwarz et al., 2008; Jung et al., 2011). In our study, miR156g appears to negatively regulate seven target genes, including SPL6 and SPL12. Novel miRNAs were equally differentially accumulated and significantly negatively correlated to their target genes, thus indicating an important regulatory role in xylem differentiation. In our study, the novel_130 is predicted to target the expression of up-regulated GCS1 in Dotz. In A. thaliana, α-Glucosidase I (encoded by GCS1) is required for cellulose biosynthesis and morphogenesis (Gillmor et al., 2002). These results indicate that novel_130 may affect cellulose synthesis in the Dotz of D. odorifera by regulating GCS1.