One-year-old P. deltoides “Danhong” (Pd) and P. simonii “Tongliao1” (Ps) were cultivated in the experimental field of the Chinese Academy of Forestry, Beijing, China (116.256°E, 40.007°N). We collected shoot apical meristem (SAM, Pd_S and Ps_S) and scraped phloem (inside of the bark, Pd_P and Ps_P) and developing xylem (newly formed xylem cells about 2–3 mm, Pd_X and Ps_X) from Pd and Ps, respectively, at diameter breast height (DBH) during the fast-growing period (July 20, 2019). Each tissue had three biological replicates. The samples (2 genotypes × 3 tissues × 3 biological replicates) used for RNA extraction were frozen immediately in liquid nitrogen and stored at −80°C. Shoot tips and cuneiform blocks (phloem, cambium, and xylem) at DBH for histologic analysis were fixed in a mixture of formalin, glacial acetic acid, and 70% ethanol in the ratio 5:5:90 vol.; FAA under vacuum for at least 24 h.

Stem pieces were embedded with Spurr resin as described by Zhang et al. (2020d). A cross section of 4 µm thick was obtained from the stem by Leica M205FA, while the SAM sections of 40 μm were obtained using a rotary microtome (Leica VT1200S, Wetzlar, Germany). The sections were stained using 0.05% toluidine blue O (TBO) and were examined with a microscope (Zeiss). The number and diameter of vessel cells in the same area (1,260 × 980 µm) were counted by ImageJ (version 1.8.0).

In order to understand the difference of wood properties between Pd and Ps, we measured the plant height and ground diameter and collected the stems to measure the wood properties including basic density, fiber length, fiber width, microfibril angle, cellulose, holocellulose, and lignin content of 1-year-old trees in December 2019. The basic density was determined by using the drainage method. A 10-cm-high wood segment was cut from the base of the trunk without bark and pith. It was softened by heating in 30% nitric acid and a small amount of potassium chlorate and converted into wood pulp by forced oscillation. The length and width of the fiber were measured 50 times using Shyygx Measure 2.0. Wood flour (40–60 mesh) from a 5-cm basal stem segment was used to determine the chemical composition. The content of holocellulose and lignin was calculated according to Chinese standards GB/T 2677.10-995 and GB/T 2677.8-1994, respectively. To evaluate the content of cellulose, the specimens were extracted with a mixed solvent of nitric acid and ethanol (v/v = 1/1) (Zhan et al., 2015). Three replicates were performed for each variety.

Total RNA was isolated from the 18 samples (SAM, phloem, and developing xylem) using the RNAprep Pure Plant Plus Kit (TIANGEN, China). An index of the reference genome (P. trichocarpa v3.0) was built using HISAT2 (Kim et al., 2015). StringTie was used to calculate FPKMs of both lncRNAs and coding genes in each sample (Pertea et al., 2015). Sequencing data are available in NCBI SRA database (SRA number: SRP2343030 to SRR13961247).

We used four filtration steps to identify lncRNAs from the transcriptome assembly: 1) Transcripts with an exon number ≥2 and length ≥200 bp were selected. 2) CuffCompare software was used to screen out transcripts that overlap with the database annotation exon field. 3) Evaluation of Coding Potential Calculator (CPC) (Kong et al., 2007), Coding Potential Assessment Tool (CPAT) (Wang et al., 2013), and Coding-Non-Coding Index (CNCI) (Sun et al., 2013) was carried out to screen whether there is coding potential. 4) They were referred to the the Hugo Gene Nomenclature Committee (HGNC) to name the novel_lncRNA of this analysis (Wright, 2014).

Two methods were used to predict lncRNA target genes. Cis target genes were predicted according to the location relationship between lncRNA and mRNA, and the screening range was within 10 kb (Jia et al., 2010). Co-expression–related target genes were predicted according to the expression correlation between lncRNA and mRNA, and the screening condition is that the correlation coefficient is greater than 0.95 and p-value< 1.68E-09. The mRNA–lncRNA regulatory network was further modeled and visualized using Cytoscape 3.8 (Otasek et al., 2019).

In order to identify lncRNAs that may be used as precursors of miRNAs, we compared the published miRNAs of P. trichocarpa in miRBase (http://www.mirbase.org/search.shtml) with lncRNAs. lncRNAs as targets of miRNAs were predicted by Novomagic, a free online platform for data analysis (https://magic.novogene.com).

To identify the differential expression of lncRNAs and mRNAs between Pd and Ps, we performed the read count of pair-wise comparisons (Pd_S vs Ps_S, Pd_P vs Ps_P, Pd_X vs Ps_X, Pd_S vs Pd_P, Pd_S vs Pd_X, Pd_P vs Pd_X, Ps_S vs Ps_P, Ps_S vs Ps_X, Ps_P vs Ps_X) by DESeq R package with a q-value < 0.05 (Love et al., 2014). Finally, those putative cis- and coexpression-targets of lncRNAs were analysed using Gene Ontology (GO) analysis (Ashburner et al., 2000; Conesa et al., 2005), and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment of DE genes was performed based on a corrected p-value < 0.05.

We selected three DELs and three DEGs from the results of the transcriptional analysis and confirmed them through qRT-PCR. PtrActin and PtrUBQ were used as internal reference genes (Wang et al., 2020a) (Supplementary Table 1). Real-time PCR was conducted on a LightCycler 480 (Roche, Basel, Switzerland) using the SYBR Premix Ex Taq™ Kit (Takara, Dalian, China). The relative expression levels of the genes were calculated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001), and the data are presented as mean ± SD from three independent biological replicates.

Wood is the secondary xylem of trees, mainly composed of cellulose, hemicellulose, and lignin. All xylem cell types first undergo secondary cell wall (SCW) thickening and programmed cell death (Zhong and Ye, 2015; Zhang et al., 2020d). The plant height and ground diameter of annual Pd were significantly higher than those of Ps (Figures 1A,B). The SAM is an important regulatory site of plant height growth, and phloem and developing xylem are important parts of plant radial growth. The slice results showed that the SAM was surrounded by young leaves, and the SAM of Pd was conical and convex, while the Ps was flat, which suggested that the growth point of Ps might not be obvious enough, which caused the high growth to be slower (Figure 1C). The radial section of the stem showed that the phloem of Pd had wider phloem fibers, and the cambium was more obvious (Figure 1D). The average diameter of the vessel cells in the xylem of Pd was 72.5 μm, which was significantly larger than that of Ps (Figures 1D,E). We further determined the quality of wood. The basic density of basic Pd was less than Ps (Figure 1F), and there was no significant difference in the microfibril angle (Figure 1G). Compared with Ps, the fiber of Pd is short and thick (Figures 1H,I). Further material property determination found that the content of the three major elements of Pd is higher (Figures 1J–L).

As an important fast-growing tree, it is very important to understand the molecular pathways of growth and development of poplar. After trimming adapters and removing low-quality and contaminated reads, in total, 246.24 Gb clean data were obtained from 18 libraries, with an average Q30 of 93.00% (Supplementary Table 2). Finally, we identified 6,355 lncRNAs, with protections of 2,454 sense_overlapping lncRNAs, 2,004 of lincRNAs, and 1,897 of antisense lncRNAs (Table S3).

In order to analysis the characteristics of these lncRNAs, we evaluated the distribution of chromosome location, transcript length, exon number, and expression level of lncRNAs. In general, lincRNA, sense_overlapping, and antisense lncRNAs were evenly distributed on 19 chromosomes, although they had different emphases (Figure 2A). The average length of lncRNAs was 990 bp, and about 63.4% contained two exons (Figures 2B,C). Antisense lncRNAs ranged in length from 201 to 9,830 bp, and the average was 940 bp. lincRNAs ranged between 201 and 4,734 bp (average = 783 bp), and average length of sense_overlapping lncRNAs was 1,197 bp. The GC content of antisense lncRNA was 41.36%, which was significantly higher than that of lincRNAs and sense_overlapping lncRNAs (Figure 2D). For expression levels, the lncRNA expression levels were different and showed fewer average counts (FPKM = 4.44) than the coding transcripts (FPKM = 18.96) (Figures 2E,F).

A principal component analysis (PCA) plot of the whole dataset revealed a sequential order of the different samples. The results showed that the SAM, xylem, and phloem of the two species were clustered into three groups, and the similarity of lncRNAs in tissues was greater than that between genotypes (Figure 3A).

In order to further analyze whether these genes were differentially expressed between the two genotypes and different tissues, nine comparative combinations were carried out. Finally, 3,572 differentially expressed (DE) lncRNAs and 27,582 DEmRNAs were obtained. Among them, the DElncRNAs of Pd_S vs Ps_S were the highest in number, including 1957 lncRNAs (Figures 3B,C; Supplementary Figure 1). DEmRNAs participated in molecular functions such as “ADP binding” and “catalytic activity” (Supplementary Figure 1C). We also compared the DElncRNAs in three tissue difference genes between the two genotypes, and there were 980 DElncRNAs in the three comparison combinations (Pd_ X vs Ps_ X, Pd_ P vs Ps_ P and Pd_ S vs Ps_ S) and 322, 333, and 420 specifically expressed lncRNAs in the SAM, xylem, and phloem, respectively (Figure 3D).

Since lncRNAs play important roles in regulating gene expression, identification and analysis of their target genes may help us explore their potential functions. We calculated and predicted 13,932 co-localization pairs consisting of 3,413 lncRNAs and 10,627 RNAs and identified 72,038 co-expression pairs consisting of 1,975 lncRNAs and 11,709 RNAs (Tables S4, S5). To further analyze the function of these lncRNAs, we performed GO and KEGG analyses on their target genes. The colocation target genes of DElncRNAs were mainly enriched in 60 GO terms such as “ADP binding” and “nucleoside binding” (Supplementary Figure 2). Some target genes were enriched in the photosynthesis pathway, including 39 lncRNAs and 68 mRNAs. For example, TCONS_00135489 showed the same trend as its target, and the expression level of related genes was high in the SAM (Supplementary Figure 3). And 3,489 genes were found by co-location and co-expression of lncRNAs. They participated in immune response, cell death, purine nucleotide binding, and ATP binding progress, and the KEGG enrichment analysis shows that they were enriched in cytochrome P450, chalcone synthase, and phenylalanine ammonia-lyase which were important parts of phenylpropanoid biosynthesis (Figure S4).

The co-expression genes of Pd_ X vs Ps_ X, Pd_ P vs Ps_ P, and Pd_ S vs Ps_ S were mainly related to “ADP binding”, “heme binding”, and other biological functions (Figure 4A). And they were significantly enriched in the “phenylpropanoid biosynthesis” and “cyanoamino acid metabolism” pathway (Figure 4B). These possible target genes provide new insights into the role of lncRNAs in poplar development.

Based on the predicted GO conditions of growth differential lncRNAs and the pathways associated with target genes, we speculated that lncRNAs might play an important role in phenylalanine biosynthesis in poplar. Phenylpropanoids are a group of plant secondary metabolites derived from phenylalanine, which has a variety of structural and signal molecular functions (Costa et al., 2003). It is the starting compound for biosynthesis of lignin, flavonoids, anthocyanins, etc., and a core mediator of crosstalk between development- and defense-related pathways (Alessandra et al., 2010). Further analysis co-expression network of these lncRNAs, and structure genes of lignins, and flavonoid biosynthesis found that TCONS_00128372, lincRNA, was located in Chr12 and interacted with MYB46, secondary wall–associated NAC domain2 (SND2), cinnamate-4-hydroxylase (C4H), caffeoyl-CoA 3-O-methyltransferase (CcoAMT), and laccase (LAC). And sense_overlapping lncRNA TCONS_00079190 co-expressed with MYB83, MYB46, NAC secondary wall thickening promoting factor1 (NST1), and LAC (Figure 5A). Similarly, the results of the study by Quan et al. (2019) and Zhou et al. (2017) showed that lncRNAs and miRNAs regulated lignin biosynthesis by regulating TFs such as bHLHs and bZIPs or directly coexpressing with structural genes in the process of lignin biosynthesis. And TCONS_00174042, TCONS_00101258, and TCONS_00136338 regulated the structural genes of flavonoid biosynthesis. TCONS_00174042 and TCONS_00101258 co-expressed with MYB3, chalcone and stilbene synthase (CHS), leucoanthocyanidin reductase (LAR), leucoanthocyanidin dioxygenase (LDOX), and naringenin 3- dioxygenase (F3H). In addition, TCONS_00136338 can coexpress with LDOX, CHS, UFGT, LAR, and dihydroflavonol 4-reductase (DFR) (Figure 5B). Lignin-related lncRNAs and structural genes were highly expressed in xylem of Pd and Ps, while flavonoid-related genes were highly expressed in the SAM (Figure S5).

Hormones such as auxin, cytokinin, and gibberellin play important roles in primary and secondary growth (Nieminen et al., 2008; Xu et al., 2019). In order to further analyze the relationship between lncRNAs and plant hormones, we constructed a co-expression network including the structural genes of auxin, cytokinin, and gibberellin biosynthesis and their possible regulatory lncRNAs (Figure 6A). For example, TCONS_00134627 could be coexpressed with GA2OX8 (Potri.011G134000), GASA10 (Potri.009G092600), and SAUR94 (Potri.009G127300) (Figure 6A). The expression levels of DEGs and DElncRNAs related to auxin, cytokinin, and gibberellin biosynthesis were shown in the heatmap (Figure 6B). The differential genes related to auxins and gibberellins were mainly concentrated in the high expression of Pd_S and Pd_X, while only a few genes are highly expressed in the phloem. And the genes related to cytokinins are highly expressed in Ps_P (Figure 6B).

To further analyze whether lncRNAs were used as ceRNAs to absorb miRNAs, inhibited the effect of miRNAs, and promoted mRNA expression, we conducted sequential analysis. MicroRNAs, a major class of small RNAs with 20–24 nucleotides, create various aspects of plant development and stress responses through posttranscriptionally regulated gene expression (Yu et al., 2019). A total of 658 miRNA–lncRNA pairs consisting of 188 miRNAs and 200 lncRNAs were identified, including 19 plant hormone–related pairs and 28 phenylalanine-related pairs (Supplementary Table 6). TCONS_00066905, hormone-related lncRNA, was predicted to be a target mimic of miR396a and miR396b. In addition, phenylalanine-related TCONS_00023606 and TCONS_00093325 were the target mimics of flavonoid-related regulatory genes miR156h and miR828a (Figure 7A). miR396-GRF was an important regulatory module of plant growth and development. We found that 12 GRFs were differentially expressed and highly expressed in the SAM of Pd, which may be an important reason for the rapid growth of Pd (Figure 7B).

In order to verify the accuracy of transcriptome data, differential genes related to xylem development and hormone signal transduction were selected for qRT-PCR verification. The results showed that the expression trends of three lncRNAs and three mRNAs in six tissues were consistent with those of qRT-PCR, which proved the reliability of transcriptome data (Figure S6).