Expression Patterns and Regulation of Non-Coding RNAs during Synthesis of Cellulose in Eucalyptus grandis Hill

: Cellulose, an essential structural component in the plant cell wall and a renewable biomass resource, plays a signiﬁcant role in nature. Eucalyptus’s excellent timber tree species (including Eucalyptus grandis Hill) provide many raw materials for the paper and wood industries. The synthesis of cellulose is a very complex process involving multiple genes and regulated by various biological networks. However, research on regulating associated genes and non-coding RNAs during cellulose synthesis in E. grandis remains lacking. In this study, the wood anatomical characteristics and chemical indexes of E. grandis were analyzed by taking three different parts (diameter at breast height (DBH), middle and upper part of the trunk) from the main stem of E. grandis as raw materials. The role of non-coding RNAs (Long non-coding RNA, lncRNA; Micro RNA, miRNA; Circle RNA, circRNA) on regulating candidate genes was presented, and the network map of ceRNA (Competing endogenous RNA) regulation during wood cellulose biosynthesis of E. grandis was constructed. The transcriptome sequencing of nine samples obtained from the trunk of the immature xylem in E. grandis at DBH, middle and upper parts had a 95.81 G clean reading, 57,480 transcripts, 7365 lncRNAs, and 5180 circRNAs. Each sample had 172–306 known miRNAs and 1644–3508 new miRNAs. A total of 190 DE-lncRNAs (Differentially expressed long non-coding RNAs), 174 DE-miRNAs (Differentially expressed micro RNAs), and 270 DE-circRNAs (Differentially expressed circle RNAs) were obtained by comparing transcript expression levels. Four lncRNAs and nine miRNAs were screened out, and the ceRNA regulatory network was constructed. LncRNA1 and lncRNA4 regulated the genes responsible for cellulose synthesis in E. grandis , which were overexpressed in 84K ( Populus Alba × Populus glandulosa ) poplar. The cellulose and lignin content in lncRNA4-oe were signiﬁcantly higher than wild type 84K poplar and lncRNA1-oe . The average plant height, middle and basal part of the stem diameter in lncRNA4-oe were signiﬁcantly higher than the wild type. However, there was no signiﬁcant difference between the growth of lncRNA1-oe and the wild type. Further studies are warranted to explore the molecular regulatory mechanism of cellulose biosynthesis in Eucalyptus species.


Determination of Cellulose with Acid Detergent
1 g sample (20-30 tons of wood powder) (W 1 ) was weighed and heated to boiling in 100 mL of acid detergent. The reflux device was turned on, kept boiling for 60 min, and transferred to the crucible (W 2 ). A paste was formed by adding 72% sulfuric acid solution (Guangzhou Deli Chemical Co., Ltd., Guangdong, China) to the material in the crucible and stirred. After 3 h, 72% sulfuric acid solution was added to half of the crucible, and the temperature was kept between 20-23 • C for continuous stirring. It was followed by rinsing with hot water until neutral, drying in an oven (103 • C, 8 h), and weighing at the constant weight (W 3 ). Then, it was kept in a muffle furnace at 550 • C for 3 h of ashing and cooling and was finally weighed (W 4 ). ADF = (W 2 − W 4 )/W 1 , Cellulose = (W 2 − W 3 )/W 1 , Lignin = (W 3 − W 4 )/W 1 , Hemicellulose = NDF − ADF, ADF: Acid washing fiber content (%), W 1 : Dry sample quality (g), W 2 : Crucible mass (g), W3: Drying quality after acid (g), W 4 : Quality of crucible and wood powder after ashing (g).

Total RNA Extraction and Detection
Total RNA was extracted according to the instructions of the TRIzol kit (Thermo Fisher Scientific, Waltham, MA, USA), and Nanodrop detected RNA concentration, 28S/18S. Then, RIN (RNA integrity number) values were detected by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and the quality and integrity of the extracted RNA were detected by 1% gel electrophoresis.

LncRNA Extract
After total RNA extraction, ribosomal RNA was removed with ribo-zero™ Magnetic Kit (Epicentre, Madison, WI, USA) to maximize the retention of all coding and ncRNAs. The RNA obtained was randomly broken into short fragments, then used as a template to synthesize the first strand of cDNA with random hexamers. Then buffer, dNTPs (dUTP instead of dTTP), DNA polymerase I, and RNase H were added to synthesize the second strand of cDNA. The second strand was purified by QiaQuick PCR kit (QIAGEN, Hilden, Germany), eluted with EB buffer solution, terminally repaired, base A added, and junction sequenced, then degraded by uracil-n-glycosylase (UNG) enzyme. Agarose gel electrophoresis was used for fragment size selection and PCR amplification. The Illumina HiSeqTM 4000 platform provided by the Gene Denovo Company (Guangdong, China) was used for sequencing.

Identification and Classification of lncRNAs
Clean outputs were obtained by filtering the machine's data and were compared with the reference genome of E. grandis (The NCBI gcf_000612305.1) using TopHat (v2.0.10) for analysis [35]. After the reads were matched with the genome, cufflinks (v2.2.1 parameter-f0.3-j 1-a 0-p 6-library-type fr-first strand) was used for assembly generation using reference annotation-based transcript (RABT) [35]. Novel transcripts were screened on the reference genome (screening criteria: transcript length ≥ 200 bp, exon ≥ 2, plant sample ≥ 1) based on their location, predicted as new lncRNAs. The coding ability of the new transcript was predicted by using the Coding protein calculator (CPC) (version CPC-0.9-r2) and Codingnon-coding Index (CNCI) [36]. Meanwhile, the new transcripts were also compared to the protein database UniProt. The intersection of these transcripts without coding potential and protein annotation information was considered a reliable prediction.

Screening of Differentially Expressed lncRNAs
Cuffdiff software [35] was used to analyze the expressions of lncRNAs and genes in Eg1, Eg2, and Eg3 (FPKM value). EdgeR [37] was used to evaluate the differences in expressions of lncRNAs. False discovery rate (FDR) and log2FC (Fold change) were used for differential screening. The screening conditions were set at FDR < 0.05 and |log2FC| > 1.

MiRNA Extraction and Identification Analysis
TRIzol was used to extract total RNA from the samples, followed by selecting 18-30nt fragments using agarose gel electrophoresis. Next, the 3 and 5 junctions were connected, and the miRNAs connected to the bilateral junctions were then subjected to reverse transcription and PCR amplification. Finally, about 140 bp strips were recovered and purified using agarose gel electrophoresis. The constructed libraries were controlled using Agilent 2100 and qPCR and sequenced by computer.
The blastall 2.2.25 (blastn) software (NCBI, Rockville, MD, USA) annotated the miRNA sequences obtained by sequencing using rRNA, scRNA, snoRNA, snRNA, and tRNA in GenBank samples and were removed as much as possible [38]. The bowtie (version 1.1.2) software was used to determine the position of the sequenced miRNA on the genome [38]. Tag sequences were used to compare the genomic location, and the mireap_v0.2 software was used to predict the unique secondary structure of miRNA and identify the new miRNA [38].

MiRNA Expression Analysis
The miRNAs identified from each sample were summarized. The expression level for each miRNA was calculated as tags per million (TPM) = T × 10 6 /N (T stands for miRNA, N stands for total miRNA (existing + existing edit + known + new predicted miRNA counts)), followed by the determination of expression profiles of all miRNAs.
EdgeR software was used to analyze the differences between various miRNAs using default parameters [39]. The screening criteria for differentially expressed miRNAs were more than two times, and the p-value was <0.05. We also analyzed the differential expressions of all miRNAs, existing miRNAs, known miRNAs, new miRNAs, and category I miRNAs. Patmatch_v1.2 software conducted the complementary pairing of miRNA and target genes through program screening and prediction [39].

Identification of circRNAs
After extracting the sequencing data, a TopHat comparison was conducted with the reference genome to obtain the comparative results for each sample. First, the unmapped reads were extracted from the alignment results, and then both ends of each unmapped read were intercepted (default 20 bp) to get the anchors. Finally, the genome was aligned with anchors, and the resulting alignment was submitted to Find_circ software to identify circRNAs [40].

Expression Analysis and Target Relationship Prediction of circRNAs
The expression of circRNAs was calculated using the back-spliced Reads Per Million mapped Reads (RPM) method, with RPM = 10 6 C/N. RPM (A) was the expression of gene A, C was the number of back-spliced Reads aligned to gene A, and N was the number of total back-spliced Reads aligned to the reference gene, p-value < 0.05, and |log2FC| > 1, for identification of differentially expressed circRNAs. GO annotation and KEGG pathway analyses were performed for differentially expressed genes. Patmatch_v1.2 software was used to analyze all the circRNAs and the miRNAs to predict the targeted relationship [39].

RNA Reverse Transcription and q RT-PCR
Two µg of RNA was taken, and reverse transcription was performed using the Novizan reverse transcription kit (R223, QIAGEN, Hilden, Germany). Fifteen ceRNA genes were selected as validation genes and expressed in three different parts of E. grandis viz., trunk, leaves, and roots. The total q RT-PCR reaction system was 20 µL containing 2 X ChamQ SYBR qPCR Master Mix 10 µL, forward primer (10 µmol·L −1 ) 0.4 µL, reverse primer (10 µmol·L −1 ) 0.4 µL, cDNA 4 µL, and RNase Free ddH 2 O 5.2 µL. The PCR amplification condition was 95 • C predenaturation for 90 s, 40 cycles of 5 s denaturation at 95 • C, 15 s of annealing at 60 • C, and 20 s of extension at 72 • C (the fluorescence signal was collected at the end of each cycle). Each analysis was repeated three times. Primers were designed by the Primer Express 2.0 Software (PE Applied Biosystems, Waltham, MA, USA) using default parameters. The primer sequences are given in Table S1. The actin gene was used as the internal reference gene, which was expressed stably during the whole growth process of E. grandis. The 2 −∆∆Ct method [41] and Microsoft Excel software (Albuquerque, NM, USA)were used to analyze the data.

Co-Expression of ceRNA Regulatory Network
The ceRNA co-expression regulatory network recognizes the interactions among the differentially expressed lncRNAs, miRNAs, and circRNAs. Pearson's correlations were used to calculate significant associations by constructing the co-expression network. Cytoscape 3.8.0 software was used to draw the networks [42].
2.13. Extraction of Total RNA from Immature Xylem of E. grandis and Reverse Transcription of cDNA The immature xylem of E. grandis preserved at −80 • C was grounded into a finely powdered mixture with liquid nitrogen. Genomic DNA was extracted from plant powder (about 50 mg) using a tree extract-DNA reagent (Wuhan Keyue Qilin Biotechnology Co., LTD., cat. 11-02, Wuhan, China). The reverse transcription reaction was performed using the TransScript II one-step gDNA Removal and cDNA Synthesis SuperMix (Beijing Chuanjin Biotechnology Co., LTD., cat. AH311-02, Beijing, China), and the cDNA finally obtained was stored at −20 • C.

Assembly of Plasmid Constructs for Plant Genetic Transformation
Positive clones were detected using colony PCR, and the plasmids were extracted and sequenced from the selected clones. The positive plasmids were transformed into Agrobacterium GV3101. The Agrobacterium colonies carrying 35S: lncRNA1 and 35S: lncRNA4, respectively, were used for poplar transformation according to the protocol by Liu et al. [44] and Lu et al. [45]. In addition, positive regeneration seedlings were identified using PCR according to Etchells [46].

Positive Transgenic Plant Screening and Phenotype Analysis
Genomic DNA was isolated from transformant 84 K poplar leaves using the MiniBEST Plant Genomic DNA Extraction Kit (TaKaRa, Bao Biological Engineering (Dalian) Co., Ltd. Liaoning, China). Positive transgenic plants were identified by PCR analysis using positive test primers (Table S3). The contents of cellulose, lignin, and hemicellulose were measured from 3-week-old plants using an ELISA (enzyme-linked immunosorbent assay) kit (Qiy, Shanghai, China). Phenotype analysis for heights, root lengths, and root numbers of 84 K poplar was carried out with 3-week-old plants. In contrast, phenotype analysis for poplar heights, internode diameters, and ground diameters was carried out with 10-week-old plants. The base and the middle of the stem were observed under a microscope (LNB701C, Shanghai LNB Instrument Co., Ltd., Shanghai, China). Four independent lines were used for the phenotype analysis of the transformants in each construct. Significant differences relative to the wild-type control plants were determined by Student's t-test (SPSS 20.0).

The Content of Cellulose, Lignin, and Hemicellulose
The cellulose content in E. grandis was significantly higher in the middle (52.9%) than in the DBH and the upper part of the trunk. Lignin content in DBH and the upper trunk was significantly higher than the middle trunk. Hemicellulose content in DBH and the middle part of the trunk were significantly higher than the upper part ( Figure 1). The cellulose content differed in different trunk positions in E. grandis, followed by non-coding transcriptome sequencing analyses.

Assembly of Plasmid Constructs for Plant Genetic Transformation
Positive clones were detected using colony PCR, and the plasmids were extracted and sequenced from the selected clones. The positive plasmids were transformed into Agrobacterium GV3101. The Agrobacterium colonies carrying 35S: lncRNA1 and 35S: lncRNA4, respectively, were used for poplar transformation according to the protocol by Liu et al. [44] and Lu et al. [45]. In addition, positive regeneration seedlings were identified using PCR according to Etchells [46].

Positive Transgenic Plant Screening and Phenotype Analysis
Genomic DNA was isolated from transformant 84 K poplar leaves using the Mini-BEST Plant Genomic DNA Extraction Kit (TaKaRa, Bao Biological Engineering (Dalian) Co., Ltd. Liaoning, China). Positive transgenic plants were identified by PCR analysis using positive test primers (Table S3). The contents of cellulose, lignin, and hemicellulose were measured from 3-week-old plants using an ELISA (enzyme-linked immunosorbent assay) kit (Qiy, Shanghai, China). Phenotype analysis for heights, root lengths, and root numbers of 84 K poplar was carried out with 3-week-old plants. In contrast, phenotype analysis for poplar heights, internode diameters, and ground diameters was carried out with 10-week-old plants. The base and the middle of the stem were observed under a microscope (LNB701C, Shanghai LNB Instrument Co., Ltd, Shanghai, China). Four independent lines were used for the phenotype analysis of the transformants in each construct. Significant differences relative to the wild-type control plants were determined by Student's t-test (SPSS 20.0).

The Content of Cellulose, Lignin, and Hemicellulose
The cellulose content in E. grandis was significantly higher in the middle (52.9%) than in the DBH and the upper part of the trunk. Lignin content in DBH and the upper trunk was significantly higher than the middle trunk. Hemicellulose content in DBH and the middle part of the trunk were significantly higher than the upper part ( Figure 1). The cellulose content differed in different trunk positions in E. grandis, followed by non-coding transcriptome sequencing analyses.

Identification and Analysis of lncRNAs in E. grandis
A large number of clean reads were obtained using high-throughput sequencing (Table S4). Reads containing rRNAs were removed, and transcriptome assembly was performed using the E. grandis genome as a reference. 5375 new transcripts were predicted by using the Coding protein calculator (CPC), 4464 transcripts were predicted by using the Coding-non-coding Index (CNCI) and 7694 new transcripts were compared to the protein database UniProt. Over 3445 new transcripts were found (Table S5) Figure 2a. New lncRNAs can be divided into five categories based on their position relative to protein-coding genes on the genome: intergenic lncRNAs, bidirectional lncRNAs, intronic lncRNAs, antisense lncRNAs, and sense overlapping lncRNAs. In our study, the percentage of intergenic lncRNAs was the highest, and no intron lncRNA was found (Figure 2b).

Identification and Analysis of miRNA in E. grandis
Clear reads were filtered to produce clear tags (Table S6). High-throughput sequencing was performed on mir-Eg-1, mir-Eg-2, and mir-Eg-3, and the clean tags obtained are shown in Table S7. MiRNA sequences were compared to the reference genome of E. grandis and plant miRNAs in miRBase to identify known miRNAs, as shown in Table S8. The hairpin structure of miRNA was predicted, and the number of new miRNAs identified in each sample was counted. More than 3000 new miRNAs were detected (a maximum of 3508) in all samples except mir-Eg-2-2 (1644) ( Table S9). A large number of clean reads were obtained using high-throughput sequencing (Table S4). Reads containing rRNAs were removed, and transcriptome assembly was performed using the E. grandis genome as a reference. 5375 new transcripts were predicted by using the Coding protein calculator (CPC), 4464 transcripts were predicted by using the Coding-non-coding Index (CNCI) and 7694 new transcripts were compared to the protein database UniProt. Over 3445 new transcripts were found (Table S5) after removing FPKM (Fragments Per Kilobase per Million) and other low coverage sequences, and the lncRNA prediction results from new transcripts are shown in Figure 2a. New lncRNAs can be divided into five categories based on their position relative to protein-coding genes on the genome: intergenic lncRNAs, bidirectional lncRNAs, intronic lncRNAs, antisense lncRNAs, and sense overlapping lncRNAs. In our study, the percentage of intergenic lncRNAs was the highest, and no intron lncRNA was found (Figure 2b).
(a) Heat maps can be clustered on two dimensions (samples, miRNA), and the topological relationship of clustering can be used to observe the relationship between samples or between miRNAs. In order to eliminate data noise, miRNAs with TPM < 1 are filtered out for heat map analysis. Based on miRNA expression, the relationship between samples and miRNA was hierarchically clustered, and the results of clustering were presented by heat map.
(e) Statistical graph of miRNA differential expression between groups.
The differential expression of lncRNAs was analyzed by comparing Eg-1, Eg-2, and Eg-3 (Figure 2c). It was observed that Eg-1 vs. Eg-2 identified one differentially expressed lncRNA (up-regulated). On the other hand, Eg-1 vs. Eg-3 identified 16 differentially expressed lncRNAs (six up-regulated, ten down-regulated). Likewise, six differentially expressed lncRNAs (one up-regulated, five down-regulated) were identified in Eg-2 vs. Eg-3. , and the topological relationship of clustering can be used to observe the relationship between samples or between miRNAs. In order to eliminate data noise, miRNAs with TPM < 1 are filtered out for heat map analysis. Based on miRNA expression, the relationship between samples and miRNA was hierarchically clustered, and the results of clustering were presented by heat map. (e) Statistical graph of miRNA differential expression between groups.
The predicted number of miRNA target genes and binding sites indicate that the relationship between miRNA and target genes was one-one, one-to-many, and many-toone (Table S10). Target genes of differentially expressed miRNAs in mir-Eg-1 vs. mir-Eg-2 samples were mainly enriched in cellular, metabolic, and single-organism processes, along with the binding and catalytic activities. Target genes of differentially expressed miRNAs in mir-Eg-1 vs. mir-Eg-3 samples were mainly associated with cellular and metabolic processes. Target genes of differentially expressed miRNAs in mir-Eg-2 vs. mir-Eg-3 samples were mainly involved in cellular and metabolic processes and binding and catalytic activities (Figure 3).

Circular RNA Statistical Analysis
CircBase database annotation was performed on circRNAs. First, annotated circRNAs were defined as existing circRNAs, while the unannotated ones were defined as new predicted circRNAs. Then, the existing circRNAs were annotated with the starBase miRNA targeting relationship and were defined as the existing targeting relationship. MiRNA targeting was predicted for all circRNAs, providing a new predictive targeting relationship. A total of 4347 miRNAs and 3871 circRNAs displayed 250,104 targeting relationships (Table S11).

Identification and Analysis of miRNA in E. Grandis
Clear reads were filtered to produce clear tags (Table S6). High-throughput sequencing was performed on mir-Eg-1, mir-Eg-2, and mir-Eg-3, and the clean tags obtained are shown in table S7. MiRNA sequences were compared to the reference genome of E. grandis and plant miRNAs in miRBase to identify known miRNAs, as shown in table S8. The hairpin structure of miRNA was predicted, and the number of new miRNAs identified in each sample was counted. More than 3000 new miRNAs were detected (a maximum of 3508) in all samples except mir-Eg-2-2 (1644) (Table S9).
The predicted number of miRNA target genes and binding sites indicate that the relationship between miRNA and target genes was one-one, one-to-many, and many-to-one (Table S10). Target genes of differentially expressed miRNAs in mir-Eg-1 vs. mir-Eg-2 samples were mainly enriched in cellular, metabolic, and single-organism processes, along with the binding and catalytic activities. Target genes of differentially expressed miRNAs in mir-Eg-1 vs. mir-Eg-3 samples were mainly associated with cellular and metabolic processes. Target genes of differentially expressed miRNAs in mir-Eg-2 vs. mir-Eg-Forests 2021, 12, x FOR PEER REVIEW 11 of 23

Circular RNA Statistical Analysis
CircBase database annotation was performed on circRNAs. First, annotated circR-NAs were defined as existing circRNAs, while the unannotated ones were defined as new predicted circRNAs. Then, the existing circRNAs were annotated with the starBase miRNA targeting relationship and were defined as the existing targeting relationship. MiRNA targeting was predicted for all circRNAs, providing a new predictive targeting relationship. A total of 4347 miRNAs and 3871 circRNAs displayed 250,104 targeting relationships (Table S11).
Twelve non-coding RNAs were expressed in three different tissues in the trunk, leaves, and roots of E. grandis. The results displayed that XR_001985124.1 and XR_727233.2 were highly expressed at the top and the root of the tree. XR_720796.2 was expressed at DBH, middle, and top of the trunk but significantly in the root. Most genes were highly expressed within the roots and leaves (Figure 4). Q RT-PCR verified high throughput sequencing results, and the expression patterns of genes were consistent with transcriptome data, depicting reliability.

Construction of ceRNA Co-Expression Regulatory Network
Based on relatedness to cellulose synthesis, ceRNA and mRNA were screened out to construct a co-expression regulatory network, with miRNA as the center, connecting different lncRNAs, mRNAs, and circRNAs. As shown in Figure 5A, four lncRNAs, 19 mRNAs, and 17 circRNAs constituted a co-expression network with novel-m0277-5p as the center. With mir6476-x as the central node, ten lncRNAs, 18 mRNAs, and two circRNAs constituted the co-expression network ( Figure 5B). A novel-m1920-3p, novel-m1950-3p, novel-m2147-5p and novel-m3931-5p constructed the central framework ( Figure 5C). A novel-m1950-3p as the center, one lncRNA, one mRNA, and six circRNAs also constituted a co-expression network ( Figure 5D), while another co-expression network was found with both novel-m2401-3p and novel-m4146-3p as the center ( Figure 5E).

Construction of ceRNA co-Expression Regulatory Network
Based on relatedness to cellulose synthesis, ceRNA and mRNA were screened out to construct a co-expression regulatory network, with miRNA as the center, connecting different lncRNAs, mRNAs, and circRNAs. As shown in Figure 5A, four lncRNAs, 19 mRNAs, and 17 circRNAs constituted a co-expression network with novel-m0277-5p as the center. With mir6476-x as the central node, ten lncRNAs, 18 mRNAs, and two circRNAs constituted the co-expression network ( Figure 5B). A novel-m1920-3p, novel-m1950-3p, novel-m2147-5p and novel-m3931-5p constructed the central framework ( Figure 5C). A novel-m1950-3p as the center, one lncRNA, one mRNA, and six circRNAs also constituted a co-expression network ( Figure 5D), while another co-expression network was found with both novel-m2401-3p and novel-m4146-3p as the center ( Figure  5E).

Gene Cloning and Vector Construction
The extracted genomic DNA is shown in Figure S2. Using a Q5 high-fidelity DNA polymerase (New England Biolabs, M0491), XR_001980078.1 (lncRNA1), XR_720796.2 (lncRNA2), XR_727233.2 (lncRNA3), and XR_001985124.1 (lncRNA4) sequences with relatedness to cellulose biosynthesis pathway were successfully amplified from the prepared genomic DNA of E. grandis ( Figure S3). pCAMBIA vector was successfully constructed. Plasmids were then cloned from positive bacteria and identified by double enzyme digestion ( Figure S4). Finally, the clones were selected for sequencing and were successfully transformed to Agrobacterium GV3101.

Gene Cloning and Vector Construction
The extracted genomic DNA is shown in Figure S2. Using a Q5 high-fidelity DNA polymerase (New England Biolabs, M0491), XR_001980078.1 (lncRNA1), XR_720796.2 (lncRNA2), XR_727233.2 (lncRNA3), and XR_001985124.1 (lncRNA4) sequences with relatedness to cellulose biosynthesis pathway were successfully amplified from the prepared genomic DNA of E. grandis ( Figure S3). pCAMBIA vector was successfully constructed. Plasmids were then cloned from positive bacteria and identified by double enzyme digestion ( Figure S4). Finally, the clones were selected for sequencing and were successfully transformed to Agrobacterium GV3101.

Effects of Overexpression of lncRNA1 and lncRNA4 on Plant Growth
To verify the functions of lncRNA1 and lncRNA4, they were transformed into 84 K poplar, viz., lncRNA1-oe, and lncRNA4-oe, respectively (Figures S5 and S6). The results showed that cellulose, lignin, and hemicellulose contents of lncRNA4-oe with three-weekold plants were significantly higher than that in the lncRNA1-oe and wild-types ( Figure 6). In addition, the average height, root lengths, and the number of roots of lncRNA4-oe of three-week-old plants were significantly higher than lncRNA1-oe and wild-types (Figure 7).
To verify the functions of lncRNA1 and lncRNA4, they were transformed into 84 K poplar, viz., lncRNA1-oe, and lncRNA4-oe, respectively (Figures S5 and S6). The results showed that cellulose, lignin, and hemicellulose contents of lncRNA4-oe with three-weekold plants were significantly higher than that in the lncRNA1-oe and wild-types ( Figure 6). In addition, the average height, root lengths, and the number of roots of lncRNA4-oe of three-week-old plants were significantly higher than lncRNA1-oe and wild-types ( Figure  7).
The average height of lncRNA4-oe of ten-week-old plants was 26.31% and 22.17% higher than the controls and lncRNA1-oe, showing significant differences (Figure 8a). At the middle of the stem, the diameters of lncRNA4-oe were significantly longer than lncRNA1-oe and controls. At the base of the stem, the diameters of lncRNA4-oe were significantly longer than lncRNA1-oe and controls (Figure 8b). In addition, the increase of xylon cells in the middle and the base part of the lncRNA4-oe was significantly higher than lncRNA1-oe and wild type 84 K poplar (Figure 9).    The average height of lncRNA4-oe of ten-week-old plants was 26.31% and 22.17% higher than the controls and lncRNA1-oe, showing significant differences (Figure 8a). At the middle of the stem, the diameters of lncRNA4-oe were significantly longer than lncRNA1-oe and controls. At the base of the stem, the diameters of lncRNA4-oe were significantly longer than lncRNA1-oe and controls (Figure 8b). In addition, the increase of xylon cells in the middle and the base part of the lncRNA4-oe was significantly higher than lncRNA1-oe and wild type 84 K poplar (Figure 9).

Discussion
The present study found 3000 new lncRNAs in each sample compared to a previous study that showed 551 new lncRNAs in the leaves and stems of E. grandis [19]. Annotation of lncRNA targets revealed their involvement in many physiological processes like photosynthesis, ubiquitin-mediated proteolysis, and protein output and processing from the endoplasmic reticulum. A total of 4,481 intergenic lncRNAs, 761 antisense lncRNAs, and 710 intron lncRNAs were identified in different stages of development in poplar [47]. In Cassava, 833 transcripts were classified into 652 intergenic and 181 antisense lncRNAs

Discussion
The present study found 3000 new lncRNAs in each sample compared to a previous study that showed 551 new lncRNAs in the leaves and stems of E. grandis [19]. Annotation of lncRNA targets revealed their involvement in many physiological processes like photosynthesis, ubiquitin-mediated proteolysis, and protein output and processing from the endoplasmic reticulum. A total of 4481 intergenic lncRNAs, 761 antisense lncRNAs, and 710 intron lncRNAs were identified in different stages of development in poplar [47]. In Cassava, 833 transcripts were classified into 652 intergenic and 181 antisense lncRNAs strictly based on genomic locations [48]. This study detected 2122 intergenic lncRNAs, 239 bi-directional lncRNAs, no intron lncRNAs, and 341 antisense lncRNAs. Studies show that lncRNAs regulate gene expression in both the nucleus and cytoplasm via diversified pathways. In the nucleus, lncRNAs modulate gene expression by affecting chromatin remodeling, epigenetic modifications, and alternative splicing [49]. Our research screened four lncRNAs, viz., XR_001980078.1, XR_720796.2, XR_727233.2, and XR_001985124.1, which could be closely related to cellulose synthesis in E. grandis. In P. tomentosa, 1377 lncRNAs were identified in the xylem of tension wood, contrast wood, and normal wood, and 15,691 lncRNAs were identified in the vascular cambium, developed xylem and mature xylem. Both UGTRL and NERDL might be involved in wood formation through regulation [16,17]. In E. grandis, lncRNA might act as a potential target or mimic to interact with miRNA, consistent with poplar lncRNA's interaction [47]. In this study, two lncRNAs related to cellulose synthesis of E. grandis were cloned, and vectors were constructed. Further experiments will be carried out about genetic transformation to verify its function.
In this study, five miRNAs, viz., mir5298-y, mir3951-x, mir5198-y, mir5156-x, and mir5298-y, were closely related to cellulose synthesis of E. grandis. Mcnair [50] found expression patterns of 12 miRNAs in typical and fast-growing Eucalyptus trees. Pappas et al. [27] used high-throughput sequencing to excavate the xylem of two E. globulus genotypes, and obtained significant miRNA information. Lin et al. [51] showed that 41 miRNAs were associated with the wood formation of E. grandis, mainly regulating transcription factors. Although many miRNAs have been identified in plants, their functions are still not fully elucidated. Future studies should explore the combined functions of miRNA and its regulatory mechanism. Furthermore, the regulatory roles of miRNAs and their interrelationships and regulatory networks should be clarified to use these miRNAs better and apply them to crop breeding and production practice through corresponding biotechnological means.
This study selected ceRNA and mRNA related to cellulose synthesis to construct coexpression regulatory networks, with miRNA as the center, connecting different lncRNAs, mRNAs, and circRNAs. CircRNAs have been found to function primarily as miRNA sponges. CDR1as, also known as cIRS-7, has been found to contain more than 70 conserved miR7-binding sites. Due to the binding of ciRS-7 with miR7, the expression levels of the miR7 target gene significantly increased due to decreased activity. Once the cyclic NRAciRS-7 was degraded, miR7 was released [56].
Furthermore, circRNAs positively regulated Pol II transcription by interacting with Pol II [55]. Studies have shown that EIciRNAs could micro-regulate the expression of parent genes by interacting with U1 snRNA. This discovery revealed the role of circRNAs in transcriptional regulation and the regulatory mechanism of specific interaction between EIciRNA and U1 snRNA [57]. In summary, the interaction of circRNAs with transcription mechanisms provides new insights into the regulatory mechanisms of gene expression in cells and ceRNA network mechanisms. In this study, lncRNA1 and lncRNA4 were overexpressed in poplar. Moreover, cellulose and lignin contents in lncRNA4-oe at three weeks of seedling age were significantly higher than those in lncRNA1-oe and wild-type 84K poplar. The hemicellulosic content of lncRNA4oe was significantly higher than that of lncRNA1-oe and wild-type poplar 84K, indicating that the cis-acting target gene CesA7 of lncRNA4 promoted cellulose and lignin synthesis. The average plant height, root length, and root number of lncRNA4-oe at three weeks of seedling age were significantly higher than wild-type poplar 84K and lncRNA1-oe. Similar results were obtained even at ten weeks of seedling age. The increase of xylem cells in lncRNA4-oe was more extensive than that in wild-type plants and lncRNA1-oe. The mechanism of lncRNA regulating target genes needs to be further explored.

Conclusions
In this study, RNA-sequence transcriptome analysis was performed on different tissues of immature xylem in the trunk to identify ceRNA (lncRNA, miRNA, circRNA) in the process of cellulose formation in E. grandis. As a result, four lncRNAs, viz., XR_001980078.1, XR_001985124.1, XR_727233.2 and XR_720796.2, and five miRNAs, such as miR5298-y, miR3951-x, miR5198-y, miR5156-x, and miR5298-y that could be closely related to cellulose synthesis were screened out. In addition, the ceRNA co-expression network related to cellulose synthesis of E. grandis was constructed, and lncRNA4 showed potential. These results have laid a foundation for further research on the expression and regulation of genes related to cellulose synthesis in E. grandis.  Figure S2: Total RNA extraction from E. grandis immature xylem M2Kb, molecular marker2000 bp, Figure S3: Electrophoresis for PCR product of lncRNA1 and lncRNA4 cloning, Table S1: Primer sequence of target gene qPCR, Table S2: Primer sequence, Table S3: Primer sequences of positive control test, Table S4: Genome comparative statistics, Table S5: Statistics on the number of transcripts, Table S6: Filtering statistics of each sample data, Table S7: The tag statistics of samples, compared with the genome, Table S8: Statistics of known miRNA number and tag abundance identified in each sample, Table S9: The number of novel miRNA identified in each sample and the abundance statistics of tag, Table S10: Predictive statistics of all miRNA target gene loci, Table S11: Statistics on new predictive targeting relationships of circRNAs.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.