Bougainvillea glabra is grown (1 age of these trees) in Nanning, Guangxi Province in China. Three stages in the development of the thorns were identified viz., as given in sample Table 1 below; please refer to the table for a detailed description of the correct stages.

In different stages of thorn formation on the stems, the thorns and stems of B. glabra were collected for transcriptome analysis. The three stages for analysis include the early stage of thorn formation, the stage in which thorns were hardened, and the stage that is between them. Triplicate samples for each stage were included.

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina), according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform, and paired-end reads were generated.

Raw data (raw reads) of fastq format were first processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing the adapter, reads containing ploy-N, and reads whose quality was low from raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. All the downstream analyses were based on clean data with high quality.

The left ends (read1 files) of fastq format were pooled into one big left.fq file, and the right end (read2 files), into one big right.fq file. The transcriptome assembly was accomplished based on left.fq and right.fq using Trinity (Grabherr et al., 2011) with min kmer cov set to 2 by default and all other parameters set to default. The assembled transcripts were hierarchically clustered to unigenes using shared reads and expressions by Corset (Davidson and Oshlack, 2014).

The correlation of gene expression levels between samples is a crucial indicator to test the reliability of an experiment and to check whether the sample selection is reasonable. It is essential to check the correlation of gene expression levels between samples before performing differential expression analysis. The Pearson correlation coefficient represents the correlation of gene expression levels between samples. The closer the correlation coefficient is to 1, the higher the similarity of expression patterns between samples. A correlation coefficient between 0.8 and 1 indicates a very strong correlation. If the correlation coefficient between samples of biological repeats is lower than 0.8, it means that the samples’ correlation is low. The correlation between pairwise comparisons of the three biological replicates at the same site in the same period is above 0.8, indicating that their respective gene expression levels are similar.

Gene function was annotated based on the NCBI non-redundant protein and nucleotide sequences. Gene function was also annotated by the Pfam (protein family), KOG/COG (clusters of orthologous groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KEGG Ortholog database (Kanehisa and Goto, 2000), and Gene Ontology (GO).

Picard-tools v1.41 and samtools v0.1.18 were used to sort, remove duplicated reads, and merge the bam alignment results of each sample. GATK3 software was used to perform SNP calling. Raw vcf files were filtered using the GATK standard filter method and other parameters (cluster: 3; window size: 35; QD < 2.0 or FS > 60.0 or MQ < 40.0 or SOR >4.0 or MQRankSum < −12.5 or ReadPosRankSum < −8.0 or DP < 10).

SSR of the transcriptome was identified using MISA (http://pgrc.ipk-gatersleben.de/misa/misa.html), and the primer for each SSR was designed using Primer3 (http://primer3.sourceforge.net/releases.php). Quantification of gene expression levels was estimated by RSEM (Li and Dewey, 2011) for each sample: 1. clean data were mapped back onto the assembled transcriptome; 2. read count for each gene was obtained from the mapping results.

Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the Differential Expression Sequence (DESeq) R package (1.18.0). DESeq provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted p-value <0.05 found by DESeq were assigned as differentially expressed.

Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted using the edgeR program package through one scaling normalized factor. Differential expression analysis of two conditions was performed using the DEGSeq R package (1.20.0). The p-values were adjusted using the Benjamini and Hochberg method. A corrected p-value of 0.005 and log2 (fold change) of 1 were set as the threshold for significantly differential expression.

GO enrichment analysis of differentially expressed genes (DEGs) was implemented using the GOseq R package, in which the gene length bias was corrected. GO terms with a corrected p-value of less than 0.05 were considered significantly enriched by DEGs. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism, and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/) (Kanehisa, 2019). We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways.

Protein–protein interaction (PPI) analysis of DEGs was based on the STRING database, with known and predicted PPIs. For the species existing in the database, we constructed the networks by extracting the target gene list from the database. Otherwise, Blastx (v2.2.28) was used to align the target gene sequences to the selected reference protein sequences, and then, the networks were built according to the known interaction of selected reference species.

Transcription factor (TF) prediction is done using iTAK software, the basic principle of which is to identify TFs by hmmscan using the categorized and defined TF (transcription factor) family and rules in the database.

Expression profiling of six genes of the phenylpropanoid pathway was carried out using the reverse transcription fluorogenic quantitative polymerase chain reaction (RT-qPCR). The total RNA from thorns and stems of B. glabra at three stages was extracted using the HiPure HP Plant RNA Mini Kit (Magen Biotechnology Co., Ltd, R4165-02). cDNA was reverse-transcribed from mRNA using HiScript II Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech Co., Ltd, R223-01). cDNA and ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd, Q711-02) were used for qPCR. Real-time quantitative PCR analysis was conducted using the qTOWER3 system (AnalytikJena, German). The primers used in RT-qPCR are provided in Supplementary Table S1.

The samples were dried at 80°C until they reached a constant weight, pulverized, passed through a 40-mesh sieve, and weighed for a certain amount (denoted as W). Acetylation of phenolic hydroxyl groups of lignin in samples was performed. Following the acetylation of phenolic hydroxyl groups within lignin, a discernible absorption peak at 280 nm became evident. Notably, the absorbance measurement at 280 nm exhibited a direct positive correlation with the lignin content present. In this study, the lignin content was characterized by the absorbance value at 280 nm. The formula for converting 280 nm absorbance to lignin content is as follows: Lignin (mg/g) = (ΔA-0.0068)÷0.0347×V×10–3÷W×T = 0.0294×(ΔA-0.0068) ÷0.002 × 50. (V, total volume of the reaction; W, sample quality (dry weight); T, dilution factor; A, measurement of absorbance at 280nm; ΔA, measurement of absorbance at 280 nm of the sample and absorbance value at 280 nm for blank control).

For research purposes, the plants of B. glabra were selected at three different stages, ranging from thorn formation (stage 1) to thorn hardening (stage 2 to stage 3), as shown in Figure 1. We extracted the total RNA from both thorns and stems and constructed transcriptome libraries for each stage using three biological triplicates of samples. We sequenced the libraries using the Illumina platform and obtained raw reads ranging from 41,266,290 to 49987188 from 18 libraries. The base calling accuracy of more than 99.9% (Q30(%)) for each treatment was over 90%, as shown in Table 2.

The transcript sequence assembled by Trinity was used as the reference sequence for subsequent analysis. After hierarchical clustering by Corset, the longest cluster sequence was obtained for subsequent analysis. The lengths of transcripts and clustered sequences were counted separately, and the results are shown in Figure 2.

Gene function was annotated using various databases including NCBI, Pfam, KOG/COG, Swiss-Prot, KEGG Ortholog, and GO. A total of 8,937 genes were co-annotated using these databases (as shown in Figure 3A). The FPKM density distribution reflects the gene expression pattern of each sample as a whole. The graph shows a non-standard normal distribution with an area of 1, indicating that the sum of probabilities is 1. The peak of the density distribution curve represents the highest number of genes at the expression level (Figure 3B). In this experiment, hierarchical clustering analysis was performed using the FPKM values of differential genes under different experimental conditions as expression levels. Genes in the same cluster exhibit similar expression level changes under various treatment conditions. Genes within the same group have similar expression patterns and might participate in the same biological processes or have similar functions. As shown in Figure 3C, the gene expression patterns in stages 2 and 3 are significantly different from those in stage 1, both in thorns and stems. Furthermore, the gene expression patterns are notably different between thorns and stems.

Compared to the thorn formation stage (stage 1), at the stage of thorns turning hard (stage 2), a total of 1,045 genes were upregulated, and 391 genes were downregulated significantly in thorns (Figure 4A), while 918 genes were upregulated and 784 genes were downregulated markedly in stems (Figure 4B). As the thorns become hard and turn brown (stage 3), a total of 63 genes exhibited a noteworthy increase in expression and 98 genes showed a significant decrease in expression in the thorns (Figure 4A), while 46 genes showed a significant upregulation in expression and 29 genes showed downregulation within the stems, compared with the stage 2 (Figure 4B). The results indicate that stage 2 is the most active period of gene expression during the hardening process of B. glabra thorns.

Active gene expression means vigorous biosynthesis and metabolism in the organism. According to GO enrichment analysis, at the stage of thorns turning hard, the three most gene-enriched biological processes in thorns (C2) were single-organism metabolism, carbohydrate metabolism, and the oxidation–reduction process, compared with the early stage of thorn formation (C1). At the same time, the molecular function exhibiting the highest degree of gene enrichment was catalytic activity (Figure 5A). The gene enrichment in stem samples was generally consistent with that observed in thorn samples (Figure 5B). Comparing the gene enrichment analysis at stage 3 with stage 2 of thorns (C3 VS C2), the bioprocess of metabolism was the most enriched, while the molecular function of oxidoreductase activity was enhanced as the top activity. The number of genes enriched in these biological processes decreased significantly when the thorns entered 2–3 stage. However, the biological processes in which genes are mainly enriched are still the same three processes—the single organism metabolism process, carbohydrate metabolism process, and oxidation–reduction process (Figure 5C). The situation of gene enrichment in the thorn (J3 VS J2) samples was comparable to what was noted in the stem samples (Figure 5D).

A directed acyclic graph (DAG) was drawn based on the results of differential gene GO enrichment analysis. Further comprehensive analysis of the biological process from stage 1 to stage 3 shows that the synthesis of cellulose is significantly enriched in this process (Figure 6). It is found that genes related to cellulose synthesis are significantly enriched in thorns during this process, whereas genes related to cellulose degradation and synthesis are significantly enriched in the stems. We speculate that degraded cellulose in the stem may be transported to the thorns and participate in their hardening process.

Through KEGG enrichment analysis, we found that phenylpropanoid biosynthesis was a key step in the hardening process of the thorns of B. glabra. From stage 1 to stage 3, in thorns, the most significantly enriched pathway was the phenylpropanoid biosynthesis (Figure 7; Supplementary Tables S2–5). From stage 1 to stage 2, many genes in phenylpropanoid biosynthesis were upregulated remarkably in thorns (Figure 8). The above results were also verified by RT-qPCR of the important enzymes of the pathway. It indicates that PAL (EC:4.3.1.24, K10775, phenylalanine ammonia-lyase) (Figure 9A), CYP73A (EC:1.14.14.91, K00487, trans-cinnamate 4-monooxygenase) (Figure 9B), 4CL (EC:6.2.1.12, K01904, 4-coumarate-CoA ligase) (Figure 9C), CCR (EC:1.2.1.44, K09753, cinnamoyl-CoA reductase) (Figure 9D), CAD (EC:1.1.1.195, K00083, cinnamyl-alcohol dehydrogenase) (Figure 9E), and POX (EC:1.11.1.7, K00430) (Figure 9F) were the key genes for the phenylpropanoid biosynthesis pathway, leading to the structural component of plant secondary cell walls and lignin biosynthesis. The expression levels of these genes exhibited a notable upregulation during the second stage compared with the first stage. The phenylpropanoid biosynthesis is the key step of lignin biosynthesis. The determination of the lignin content in thorns and stems at three stages (Figure 9G) also agreed with the results of transcriptome analysis and RT-qPCR. To sum up, phenylpropanoid biosynthesis is a key step in the hardening process of the thorns of B. glabra as it was upregulated during maturation with maximum gene expression of the key enzymes during stage 2 of maturation.

From the first stage to the third stage of thorn formation, 2,395 TFs were completely identified in thorns. They included MYB, AP2-EREBP, orphans, C3H, HB, bHLH, WRKY, NAC, C2H2, and bZIP (Figure 10A). In stems, 2,444 TFs were completely identified. They included similar types of TFs as thorns, differing only in the percentage of TFs in each class from thorns (Figure 10B). From the first stage to the second stage of the thorn formation, 39 genes were significantly altered in the phenylpropanoid biosynthesis pathway in thorns (Supplementary Table S2). By correlation analysis, there were many TFs with correlation coefficients greater than 0.8 with these genes. These TFs MYB and HB are the most numerous (Figure 10C).