Bract samples used in this study were obtained from cuttings (3–5 years old) of B. glabra cultivars or accessions cultivated in the nursery of the Institute of Bougainvillea in Yuanshan, Zhangzhou, China. For each examined sample, mature bracts on several branches from a single accession were used. The information on the sampled accessions was listed in Supplementary Table S1 .

In order to quantify the variation of bract coloration, the CIE L-a-b value (i.e., L, representing the brightness dimension; a, value representing the red and green dimension; b, representing the yellow and blue dimension) of standard color space model was used to measure the bract surface color using the AvaSpec-ULS2048CL-EVO high-speed CMOS spectrometer [Avantes (Shanghai, China) Co., Ltd.] (Ogutcen et al., 2020). The measurement of L-a-b color values for each accession was repeated four times. The average value (AL, Aa, Ab) of each L-a-b value was calculated as color trait data.

The RNA sequencing (RNA-seq) data and de novo transcriptome assembly were downloaded from the Genome Sequence Archive (GSA), National Genomics Data Center (NGDC), under Bioproject number PRJCA011746 (Supplementary Materials).

The transcriptome reads of the 18 B. glabra accessions were mapped onto the de novo assembled pan-transcriptome assembled transcripts using Bowtie2 (Langmead and Salzberg, 2012). The variants [including single-nucleotide polymorphisms (SNPs) and Indels] were called from mapped BAM files using the GATK pipeline with HaplotypeCaller model (McKenna et al., 2010). The variants were filtered using the following parameters: 2 < depth (DP) < 100, minQ >20, maximum missing rate <40%, and minor allele frequency (MAF) >2%. The filtered set of variants was then used for the downstream population genomic analysis.

A maximum likelihood (ML) tree of the 18 B. glabra accessions was constructed with filtered variant sets using the software IQ-TREE 2.1.3 (Lam-Tung et al., 2015). Iterative tree construction was performed using the best model, PMB+F+R3 protein model, according to the Bayesian Information Criterion and an ultrafast bootstrap approximation algorithm. The final ML tree was visualized using the software iTOL v4 (Letunic and Bork, 2019).

To further analyze the genetic relationship among the 18 accessions, the population admixture was constructed using admixture v1.3.0 (Alexander et al., 2009). The optimal number of clusters (K) was identified based on the cross-validation (CV) error value, which was tested from K = 1 to 10. After choosing the optimal K value, the admixture matrix was visualized using iTOL v4.

The genetic relationship matrix among the 18 accessions was calculated using the software GCTA (Yang et al., 2011). The eigenvalues and eigenvectors were generated for PCA. The top 2 principal components were used to generate PCA scatter plot by ggplot2 in R (RAM and ZJ, 2019).

The expression levels of unigenes for bracts among the 18 accessions were quantified using the software RSEM (Dewey and Li, 2011) using BAM files generated from read mapping and normalized to the transcript per million (TPM) matrices. The statistical distribution of expression was constructed in box plot. The pheatmap package in R (Kolde and Kolde, 2018) was used to construct a heatmap with clustering tree for the correlation of expression levels among the 18 accessions. The PCA dimensionality reduction analysis was also performed on the expression levels of bracts to identify the clusters of the 18 accessions with high expression similarity.

In DEG analysis, expression data of three accessions with the highest AL values vs. three accessions with the lowest AL values were used for comparison in DESeq2 (Love et al., 2014) with multiple test correction method “BH.” Significantly DEGs with P-adjust values <0.05 and two times up/down fold changes were screened. Expression data of accessions with the highest and the lowest Aa and Ab values were also processed in a similar way.

To determine the correlation between bract phenotypical color variation and their gene expression levels, a gene expression network based on color variant trait data and gene TPM expression data of bracts from the 18 accessions was constructed using the WGCNA pipeline (Langfelder and Horvath, 2008). After the invalid insignificant expressed genes were filtered out from the TPM data, soft threshold of unsigned network was calculated to match the subsequent correlation analysis. Module detection was performed by Topological Overlap Matrix (TOM)-based similarity measurements, and the minimum number of genes per module was set at 50 genes. The module genes were then divided, while the similar modules were merged. For the module–trait correlation analysis, the eigengene value was calculated for each module to test the association with each average L-a-b color trait. The results of module–trait correlation were then constructed into heatmap to pick essential modules for subsequent analysis.

The module genes with high correlation with color trait data were exported to calculate the weight nodes and screen the hub genes in the cytoHubba plug-in of Cytoscape software (Shannon et al., 2003). Maximal Clique Centrality (MCC) algorithm was used to calculate the top 50 weighted genes in picked module and saved to map former annotation. After essential modules were selected and hub genes were calculated, a co-expression network was constructed using software Gephi (Bastian, 2009). We uploaded each picked module gene and hub gene with their annotation. Edge degree and module classification was calculated in Gephi, and small nodes without annotation and high edge degree were filtered. Finally, the co-expression network for each module was constructed.

The selected genes from the DEG analysis and the genes in the essential modules from WGCNA were enriched based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database using clusterProfiler package in R (Yu et al., 2012). The top 20 terms of GO and KEGG enrichment were visualized in dotplot. The significance levels of GO and KEGG enrichment were tested with standard of P < 0.05.

The candidate genes from DEG analysis and the hub genes in the top WGCNA modules that correlated with L-a-b values were selected for expression pattern analysis. The expression patterns of selected genes were visualized using pheatmap package in R (Kolde and Kolde, 2018) according to the order of each bract L-a-b value from high to low among the 18 B. glabra accessions.

Real-time qPCR (RT-qPCR) was designed to verify the expression trend of candidate genes, and they were compared with transcriptome expression data (Bustin et al., 2009). Candidate genes screened out from DEGs and WGCNA were subjected to RT-qPCR experiment. Total RNA from bract samples with high, medium, and low L-a-b values was extracted for RT-qPCR using a Tiangen polysaccharide and polyphenol plant total RNA extraction kit [Tiangen Biotech (Beijing) Co., Ltd.]. The primers were designed by the Primer Premier 5.0 (Lalitha, 2000) ( Supplementary Table S2 ). cDNA was synthesized using a TaKaRa PrimeScript RT Reagent Kit with a gDNA Eraser (Perfect Real Time) kit. RT-qPCR was performed on the CFX 96™ Real-Time PCR Detection System. The relative expression data were quantified by the 2^−△△CT method and normalized by 18S rRNA (Schmittgen and Livak, 2008). The reaction mixture (20 μl) contained 10 μl of TB Green Premix Ex Taq II (Tli RNaseH Plus) (2×), 0.5 μl of PCR Forward Primer (10 μM), 0.5 μl of PCR Reverse Primer (10 μM), 1 μl of DNA template, and 8 μl of ddH₂O. Subsequently, the PCR programs were conducted with a pre-denaturation at 95°C for 3 min; 40 cycles of 95°C for 10 s, 60°C for 30 s; and a final 65°C–95°C slow increase for melting curve analysis. Finally, T-tests were performed to detect the significance between differential expressions between samples.

The color variation of bract surface among the 18 accessions was quantified and represented by L-a-b values ( Supplementary Tables S1, S3 ) and visualized as a three-dimensional (3D) scatter plot ( Figure 1 ). In each dimension of L-a-b, the color data presented a linear variation. On the red-green dimension (trait Aa), F_06_SZ had the highest red value with 44.81 and C_12_DW had the highest green value with -1.983. On the yellow-blue dimension (trait Ab), O_07_DP had the highest yellow value with 15.38 and F_06_SZ had the highest blue value with -24.88. On the brightness dimension (trait AL), C_12_DW had the brightest value with 67.90 and F_06_SZ had the darkest value with 43.04. According to the measurement result, as |Aa+| (red) and |Ab-| (blue) values increased, the bract color tended to change from light purple to dark purple. When |Aa+| and |Ab+| (yellow) values increased, the bract color tended to change from pink to orange. When |Aa-| (green) and |AL+| (brightness) values increased, the bract color changed from white to bright green, giving it a leaflike appearance. Appropriately, F_06_SZ is the darkest and deepest purple accession and C_12_DW is the brightest whitest accession without any red pigmentation. Accessions O_07_DP and L_21_SP had top |Ab+| value, and the |Aa| value was not less than 10 that gave them orange coloration with a little red blended pigmentation.

The phylogenetic tree clustered the 18 accessions into three major subgroups: Group I, Group II, and Group III ( Figure 2A ). Group I (C_02_SZ, L_19_DP, M_02_SP, and M_05_DP) clustered dispersedly. The accessions in Group I have a similar bract shape, which was narrowly elongated and sharply pointed at the front, while their bract color varied from white to pink-purple. Group II (Out_03_DZ, F_06_SZ, L_21_SP, N_07_SP, and O_07_DP) clustered apparently. Group II accessions also have a similar bract shape, which was mostly bigger and round overall, slowly pointed at the front, and slowly convergent at the base. Group III (C_12_DW, D_31_SZ, L_02_SZ, L_22_SP, M_04_SP, N_20_SP, O_08_SP, Out02_DP, and R_01_SW) clustered densely. The accessions in Group III mostly have round overall bracts with color varying from pink-purple to white.

Admixture analysis showed that the optimal clusters of the 18 accessions are K = 2 and K = 3 as revealed by the lowest CV error ( Figure 2B ). When K = 2, Group III of the phylogenetic tree mostly consisted of S1 stratification in admixture groups. Group I and Group II of the tree consisted of two stratifications, S1 and S2. When K = 3, most of the accessions in Group III only consisted of S1 stratification. Group I and Group II of the tree consisted of one, two, or three stratifications. The accessions N_07_SP and L_21_SP of Group II consisted of only S2 stratification, while C_02_SZ and M_05_DP of Group I consisted of only S3 stratification.

The clustering results of PCA were consistent with that of the phylogenetic tree and population structure. According to the PCA scatter plot of the 18 accessions ( Figure 2C ), the first two principal components in the PCA explained 20.20% and 13.50% of variation, respectively, with a combined 33.70% of the total variation.

To detect the expression variation and their relationships among B. glabra accessions, we analyzed the expression distribution and correlation of all unigenes based on the pan-transcriptome assembly of the 18 accessions. The expression distribution was revealed by the log₁₀(TPM) values in box plot ( Figure 3A ). The expression TPM level of each accession varied between -2 and 6. A total of 23,857 unigenes were expressed in all accessions, 11,100 unigenes were expressed in the most orange bract accessions O_07_ DP and L_21_SP, and 9,188 unigenes were expressed only in C_12_DW, which has the whitest bracts.

The correlation analysis of expression level showed that the 18 accessions were clustered into four subgroups ( Figure 3B ). Subgroup I contained accessions F_06_SZ, M_04_SP, Out_03_DZ, and L_22_SP, among which three of them had the deepest purple bracts. Subgroup II contained accessions M_02_SP, N_20_SP, C_02_SZ, M_05_DP, and O_07_DP, among which three of them had light purple bracts. Subgroup III consisted of accessions C_12_DW, L_02_SZ, Out02_DP, and R_01_SW, which had two white bract accessions. Subgroup IV contained accessions L_21_SP, N_07_SP, L_19_DP, D_31_SZ, and O_08_SP, among which three of them had orange bracts.

According to the subgroup classification in the correlation analysis, the PCA clustering of expression levels among the 18 accessions ( Figure 3C ) showed that Subgroup I separated dispersedly and most accessions are purple bracts. Several accessions of Subgroup II clustered together on PC1 dimension but dispersedly along PC2 dimension whose bract color varied from purple to orange. Subgroup III clustered densely that were mostly purple bract and white bract accessions. Subgroup IV contained three accessions with orange bract relatively clustered together on PC2 but dispersedly on PC1 dimension.

To screen the gene response to the bract color variation, accessions with the highest and lowest L-a-b values were selected to conduct the DEG analysis. By comparing the gene expressions of the three accessions (N_20_SP, C_12_DW, and R_01_SW) with the lowest Aa values and the three accessions (F_06_SZ, M_04_SP, and Out_03_DZ) with the highest Aa values, 101 significantly upregulated/downregulated genes (76 upregulated/25 downregulated) were identified ( Figures 4A, D ; Supplementary Table S4 ). These genes were mostly involved in the extracellular region (gene number/P value, 8/0.00), cell wall (6/0.00), external encapsulating structure (6/0.00) according to GO enrichment ( Supplementary Figure S1A left) and were enriched in plant–pathogen interaction (3/0.02), cyanoamino acid metabolism (2/0.00), and pentose and glucuronate interconversions (2/0.01) according to KEGG enrichment ( Supplementary Figure S1A right).

By comparing the gene expressions of the three accessions (F_06_SZ, M_04_SP, and Out_03_DZ) with the lowest Ab values and the three accessions (O_07_DP, L_21_SP, and L_19_DP) with the highest Ab values, 196 significantly upregulated/downregulated genes (75 upregulated/121 downregulated) were identified ( Figures 4B, E ; Supplementary Table S4 ). These genes were mostly involved in the extracellular region (6/0.01), monooxygenase activity (5/0.00), and cell wall organization or biogenesis (4/0.01) according to GO enrichment ( Supplementary Figure S1B left) and were enriched in plant–pathogen interaction (5/0.00) and phenylpropanoid biosynthesis (4/0.00) according to KEGG enrichment ( Supplementary Figure S1B right).

According to the comparisons of the gene expressions for the three accessions (F_06_SZ, M_04_SP, and Out_03_DZ) with the lowest AL value and the three accessions (C_12_DW, R_01_SW, and M_02_SP) with the highest L value, 119 significantly upregulated/downregulated genes (23 upregulated/96 downregulated) were identified ( Figures 4C, F ; Supplementary Table S4 ). These genes were mostly involved in the extracellular region (11/0.00), catabolic process (8/0.01), and defense response (6/0.00) according to GO enrichment ( Supplementary Figure S1C left) and were enriched in sesquiterpenoid and triterpenoid biosynthesis (2/0.01) and pentose and glucuronate interconversions (2/0.02) according to KEGG enrichment ( Supplementary Figure S1C right).

The L-a-b values as the trait data combined with a total of 19,484 gene expression data of bracts for the 18 accessions were subjected to WGCNA. We picked the soft threshold (power) at 7, which passed the scale independence and mean connectivity standard for the network analysis. Gene expression was detected as an outlier based on the “ward.D2” module with all accessions passed. A total of 40 modules were constructed in the WGCNA ( Figure 5A ). Modules MEdarkgrey, MEred, and MEcyan had more than 60% correlation (Pearson correlation coefficient) with trait Aa but low to -1% for traits Ab and AL. Modules MEyellowgreen, MEsaddlebrown, and MEorange had more than 50% correlation with trait Ab. Only module MEdarkolivegreen had more than 50% correlation with trait AL. The genes in the modules that had more than 60% correlation with bract color Aa value were constructed into a co-expression network by Gephi ( Figure 5B ) and calculated by the cytoHubba plug-in in Cytoscape based on the MCC algorithm.

The top 50 genes in the co-expression network were selected as the hub genes of each module according to the connect weight values calculated by Cytoscape. A total of 96 notable annotated hub genes from co-expression networks of the top 3 modules were screened ( Figure 5B ; Supplementary Table S5 ). For trait Ab, the genes in the modules that had more than 50% correlation with bract color Ab value were constructed into the co-expression network. A total of 98 annotated hub genes from co-expression networks of the top 3 modules were selected ( Figure 5C ; Supplementary Table S5 ). For trait AL, the genes in the modules that had more than 30% correlation with L value were constructed into the network. A total of 62 annotated hub genes from co-expression networks of the top 2 modules were selected ( Figure 5D ; Supplementary Table S5 ).

Genes in the modules that had more than 50% Pearson correlation coefficient with each color trait value and P value <0.05 were extracted for GO and KEGG enrichment ( Figure 6 ). For trait Aa, the genes in modules MEdarkgrey (Pearson correlation coefficient/P value, 0.65/0.003), MEcyan (0.61/0.007), and MEred (0.62/0.006) were selected for enrichment analysis. GO enrichment indicated that 1,621 genes in these top modules were mainly involved in oxidoreductase activity (gene number/P value, 9/0.00), terpene synthase activity (6/0.00), GTPase activity (5/0.00), and response to UV-B (3/0.00) ( Figure 6A left). KEGG enrichment indicated that 1,621 genes were mainly enriched in phenylpropanoid biosynthesis (7/0.01), plant–pathogen interaction (6/0.01), and sesquiterpenoid and triterpenoid biosynthesis (4/0.00) ( Figure 6A right).

For trait Ab, genes in the modules MEorange (0.52/0.03), MEsaddlebrown (0.58/0.01), and MEyellowgreen (0.66/0.003) were selected for enrichment analysis. GO enrichment indicated that 473 genes in these top modules were mainly involved in catalytic activity (5/0.00), ATPase activity, coupled to transmembrane movement of substances (3/0.04), endopeptidase inhibitor activity (2/0.01), and S-adenosylmethionine-dependent methyltransferase activity (2/0.02) ( Figure 6B left). KEGG enrichment indicated that 473 genes were mainly enriched in the ABC transporters (3/0.00), mitogen-activated protein kinase (MAPK) signaling pathway–plant (2/0.00), and starch and sucrose metabolism (2/0.01) ( Figure 6B right).

For trait AL, genes in the module MEdarkolivegreen (0.53/0.02) were selected for enrichment analysis. GO enrichment indicated that 129 genes in this module were mainly involved in cytochrome-c oxidase activity (3/0.00), response to auxin (2/0.01), and response to stimulus (1/0.01) ( Figure 6C left). KEGG enrichment indicated that 129 genes were mainly enriched in tropane, piperidine, and pyridine alkaloid biosynthesis (2/0.00), alanine, aspartate, and glutamate metabolism (1/0.03), and aminoacyl-tRNA biosynthesis (1/0.05) ( Figure 6C right).

Candidate DEGs between accessions with high vs. low L-a-b values might potentially be involved in the regulation of bract coloration in B. glabra. In DEGs of high vs. low Aa values, one peroxidase 42 (E1.11.1.7) gene, one (-)-germacrene D synthase-like (GERD) gene, one chalcone synthase (CHS) gene, and one S-adenosylmethionine synthase 2-like (metK) gene were selected to show the consistency of expression variation with variation in bract Aa value ( Figure 7A ; Supplementary Table S6 ). In the accessions F_06_SZ, M_04_SP, and Out03_DZ with the highest Aa value, these genes were generally upregulated, while in the accessions N_20_SP, R_01_SW, and C_12_DW with lower Aa value, these genes were generally expressed less ( Figure 7A ). Moreover, in DEGs of high vs. low Ab values, one primary-amine oxidase (AOC3) gene, one cytochrome P450 (CYP716A) gene, and one cytochrome P450 76T24-like (CYP76C) gene were selected. In accessions O_07_DP, L_21_SP, L_19_DP, and N_07_SP with the highest Ab value, these genes were upregulated considerably, while in accessions F_06_SZ, Out03_DZ, M_04_SP, and C_02_SZ with lower Ab value, genes were generally downregulated ( Figure 7B ). For DEGs of high vs. low AL values, only E1.11.1.7 gene was selected. In the accessions R_01_SW and C_12_DW with the highest AL value, E1.11.1.7 was upregulated, while expressed less in other accessions ( Figure 7C ).

In the top WGCNA module MEdarkgrey correlated with trait Aa, two stay-green protein (SGR) genes, RING finger and CHY zinc finger domain-containing protein 1 (RCHY1) gene, and one GERD were selected ( Figure 8A ; Supplementary Table S6 ). In module MEred, two xanthoxin dehydrogenase (ABA2) genes, one metK, ditrans,polycis-polyprenyl diphosphate synthase (DHDDS) gene, and chlorophyll b reductase (NOL) gene were selected to show the consistency of expression variation with variation in bract Aa value. In module MEcyan, one ubiquitin-conjugating enzyme E2 A (UBE2A) gene, one allene oxide cyclase (AOC) gene, one molecular chaperone HtpG (HSP90A) gene, and cinnamyl-alcohol dehydrogenase (K22395) gene were selected. In the accessions F_06_SZ and Out03_DZ with the highest Aa value, these genes were generally upregulated, while in the accessions R_01_SW and C_12_DW with lower Aa value, these genes were generally expressed less ( Figure 7A ).

For trait Ab in WGCNA module MEyellowgreen, one signal peptide peptidase-like 2B (SPPL2B) gene, one transcription factor C subunit 3 (TFC3) gene, one phospholipase C (PLC) gene, and one acylaminoacyl-peptidase (APEH) gene were selected to show the consistency of expression variation with variation in bract Ab value ( Figure 8B ; Supplementary Table S6 ). In module MEorange, ATP-binding cassette, subfamily A (ABC1) and member 3 (ABCA3), vacuolar protein sorting 34 (PIK3C3), one protein-lysine N-methyltransferase EEF2KMT (EEF2KMT) gene, one UBX domain-containing protein 7 (UBXN7) gene, one type I protein arginine methyltransferase (CARM1) gene, and one UDP-glycosyltransferase (FG3) gene were selected. In module MEsaddlebrown, one small subunit ribosomal protein S29e (RP-S29e) gene, one glutathione S-transferase T1 (GST) gene, one urate oxidase (uaZ) gene, one sphingosine kinase (SPHK) gene, and one ATP-dependent RNA helicase DHX36 (DHX36) gene were selected. In accessions O_07_DP, L_21_SP, L_19_DP, and N_07_SP with the highest Ab value, these genes were upregulated considerably, while in the accessions F_06_SZ, Out03_DZ, M_04_SP, and C_02_SZ with lower Ab value, genes were generally downregulated ( Figure 7B ).

For trait AL in module MEdarkolivegreen, three COX1, COX2, and COX3, one cytochrome P450 82G1 monooxygenase (CYP82G1) gene, and one F-type H+-transporting ATPase subunit a (ATPeF0A) gene were selected to show the consistency of expression variation with variation in bract AL value ( Figure 7C ; Supplementary Table S6 ). In module MEsteelblue, four jasmonate ZIM domain-containing protein (JAZ) genes, one dCTP diphosphatase (DCTPP1) gene, and one tryptophan synthase beta chain (trpB) gene were selected. In the accessions R_01_SW and C_12_DW with the highest AL value, these genes were generally upregulated, while genes were generally expressed in low levels in other accessions ( Figure 7C ).

Among genes screened from DEG analysis, a GERD gene (TRINITY_DN19569_c0_g1) upregulated in accessions of high L-a-b values was chosen for RT-qPCR validation in C_12_DW, D_31_SZ, and F_06_SZ accessions following increased Aa value ( Figure 8A left). RT-qPCR result showed that the trend of GERD gene expression was similar as that in transcriptome expression analysis ( Figures 7A , 8A left). Compared with C_12_DW and D_31_SZ, the F_06_SZ (T-test, P < 0.01; P < 0.001) exhibited increased expression of GERD. An AOC3 gene (TRINITY_DN6496_c0_g2) upregulated in high Ab value accessions was also selected for RT-qPCR validation in F_06_SZ, L_22_SP, and O_07_DP accessions following increased Ab value ( Figure 8B left). Results showed a reduced expression in L_22_SP (T-test, P < 0.0001) compared with F_06_SZ and an increased expression in O_07_DP (T-test, P < 0.001) compared with F_06_SZ. A peroxidase (E1.11.1.7) gene (TRINITY_DN1078_c0_g1) upregulated in high L value accessions was chosen to perform RT-qPCR validation in F_06_SZ, O_07_DP, and C_12_DW accessions following increased AL value ( Figure 8C left). The results showed an increased expression in O_07_DP (T-test, P < 0.01) compared with F_06_SZ and a reduced expression in C_12_DW (T-test, P < 0.01) compared with F_06_SZ.

For the hub genes from each co-expression network, SGR gene (TRINITY_DN12271_c0_g2) was chosen for RT-qPCR validation in C_12_DW, D_31_SZ, and F_06_SZ accessions following increased Aa value ( Figure 8A right). The trend of relative expression was consistent with the trend in transcriptome expression analysis ( Figures 7A , 8A right). Compared with C_12_DW, SGR was significantly upregulated (T-test, P < 0.001) in F_06_SZ and D_31_SZ accessions. GST gene (TRINITY_DN671_c0_g2) was selected for RT-qPCR validation in F_06_SZ, L_22_SP, and O_07_DP accessions following increasing Ab value. The expression of GST in O_07_DP (T-test, P < 0.01) was increased compared with F_06_SZ but no significant difference between F_06_SZ and L_22_SP ( Figure 8B right). JAZ gene (TRINITY_DN16399_c0_g1) was also selected for RT-qPCR validation in F_06_SZ, O_07_DP, and R_01_SW accessions following increased AL value. The JAZ gene showed reduced expression in O_07_DP (T-test, P < 0.001) compared with F_06_SZ but significant upregulation in R_01_SW compared with F_06_SZ (T-test, P < 0.001) ( Figure 8C right).