## 3 Results
### Functional Annotation  
Functional annotation of the assembled unigenes was performed to classify their potential functions in slash pine. To this end, sequence alignment based on similarity was carried out with the several public databases, including NR, Swiss-Prot, GO, COG, KOG, KEGG, and Pfam using BLAST software with E-value < 1 \(\times\) 10\({}^{-5}\). As a result, a total of 30,330 unigenes were annotated, with 9782; 16,225; 11,170; 17,782; 22,178; 18,980; 27,598; and 27,405 genes distributed in the NR, Swiss-Prot, GO, COG, KOG, KEGG, and Pfam  
Figure 1: The overall DEGs between the high- and low-yielding slash pines. **(A)** Detailed sampling sites of high- and low-yielding slash pines for RNA-seq. White bar = 10 cm. **(B)** Principal component analysis (PCA) of 24 samples based on the DEGs (Log2 \(\geq\) 1 or \(\leq\)\(-\)1). **(C)** Venn diagram and histogram. The Venn diagram shows the number and distribution of DEGs between high- and low-yielding slash pines in April, July, and October, and the histogram indicates DEGs numbers. The blue stands for the total DEGs, red means the up-expressed DEGs, and green is the repressed DEGs. **(D)** Heat map exhibits the congregation of relative expression levels in different groups.  
databases, respectively. In addition, 11,417 unigenes were identified with 300 nt \(\leq\) lengths \(\leq\) 1000 nt, and 18,913 unigenes exhibited the length \(\geq\)1000 nt.  
To investigate the DEGs in relation to oleoresin yield of slash pine, the RNA-seq data of 24 samples were classified into six sequencing libraries based on high- or low-yielding slash pine. The changes in the transcriptome between high- and low-yielding slash pine at different oleoresin-yielding stages (April, July, and October) were analyzed. Three comparisons of P4-L vs. P4-H, P7-L vs. P7-H, and P10-L vs. P10-H (H or L refer to high or low oleoresin-yielding slash pine) were conducted (Supplementary Materials Tables S3-S5). In P4-L vs. P4-H, a total of 39 DEGs were identified with 23 upregulated and 16 downregulated genes. In P7-L vs. P7-H, there were 39 DEGs identified with 17 upregulated and 22 downregulated genes. Significantly, 291 DEGs were identified between P10-L and P10-H, with 102 upregulated and 189 downregulated genes. In Figure 1C, the Venn diagram shows that there are 30, 26, and 274 unique DEGs in P4-L vs. P4-H, P7-L vs. P7-H, and P10-L vs. P10-H groups, respectively. Only three DEGs are commonly expressed across all three compared groups. The list of all DEGs of P4-L vs. P4-H, P7-L vs. P7-H, and P10-L vs. P10-H is presented in Supplementary Materials Table S4.  
The functions of these DEGs were further annotated by aligning their sequences with public databases. According to GO enrichment analysis, all of the DEGs between high and low oleoresin-yielding slash pine were enriched to three classes, including biological process, cellular components, and molecular function. In the biological process, the significant two enriched GO terms are \"metabolic process\" and \"cellular process\". In cell component, the primary term is \"cell and cell part\", and the largest proportion in molecular function is \"catalytic activity and binding\".  
In KEGG enrichment analysis, a total of 20, 9, and 111 DEGs in P4-L vs. P4-H, P7-L vs. P7-H, and P10-L vs. P10-H were annotated, respectively. The DEGs of P4-L vs. P4-H and P7-L vs. P7-H were significantly enriched in the \"galactose metabolism\" pathway. In regard to the terpenoid biosynthesis (Supplementary Materials Table S6), DEGs of P4-L vs. P4-H were also enriched in \"ditrepenoid biosynthesis\" and \"carotenoid biosynthesis\" pathways, while DEGs of P10-L vs. P10-H were significantly enriched in \"ubiquinone and other terpenoid biosynthesis\" pathway. However, DEGs of P7-L vs. P7-H were not enriched in any known terpenoid biosynthesis-related pathways.