QHP with bright purple leaves, L2025 with green leaves, and ZSY with bright purple and green leaves were planted in the experimental field of Nanjing Botanical Garden, Memorial Sun Yat-Sen (32°3′N, 118°49′E). The fourth and fifth fully expanded mature leaves were collected from two-year-old seedlings of QHP (P) and L2025 (G), respectively. Fully expanded purple leaves (F_P) and green leaves (F_G) from two-year-old ZSY seedlings were also harvested. The leaf samples were immediately frozen in liquid nitrogen and stored at −80°C until use for the determination of chlorophyll, carotenoid, and anthocyanin contents, metabolite detection, RNA-seq, and qRT-PCR analyses. Three independent biological replicates were used for each experiment.

The contents of chlorophyll, carotenoid, and anthocyanin of the QHP, ZSY, and L2025 leaves were determined based on a previously described method (Tian et al., 2021b). Fresh leaves were used to measure total chlorophyll content based on the method described by Wu et al. (2020). Leaf tissue (0.1 g) was ground into powder and extracted in 5 mL of 95% ethanol at 50°C for 2 h. The mixture was vortexed and centrifuged at 5,000 rpm for 5 min. The absorbance of the supernatant was measured at 470, 649, and 665 nm using an ultraviolet–visible spectrophotometer (Wu et al., 2020; Zhuang et al., 2021b). Fresh leaves were used to determine the carotenoid content in accordance with the method used by (Gao et al., 2021). Approximately 100 mg of fresh leaf tissue were cut into pieces with scissors, and extracted in 10 mL of 1% (v/v) HCl–ethanol at 60°C for 30 min. The mixture was vortexed and centrifuged at 13,000 × g for 5 min. The absorbance of the supernatant was measured at 530, 620, and 650 nm using an ultraviolet-visible spectrophotometer (Gao et al., 2021).

The total anthocyanin content was determined in accordance with the method of Xu et al. (2017) using functional leaves (the third to fifth leaves from the main branches) of two-year-old seedlings of QHP, ZSY, and L2025.

Metabolite profiling of the extracted leaf samples was performed using a liquid chromatography–tandem mass spectrometry (LC-MS/MS) (Ma et al., 2006). Chromatographic separation of anthocyanin metabolites was performed on an ACQUITY UPLC HSS T3 column (100 mm × 2.1 mm, 1.8 µm), using a gradient elution of 0.1% formic acid (Solution A) with acetonitrile (Solution B) at a flow rate of 0.2 mL/min and temperature of 40°C. The elution profile was set as follows: 80:20 VA/VB at 0 min, 80:20 VA/VB at 1 min, 25:75 VA/VB at 3 min, 25:75 VA/VB at 4 min, 80:20 VA/VB at 4.05 min, 80:20 VA/VB at 6.5 min, and 80:20 VA/VB at 7.5 min. The injection volume was 10 μL.

A quadrupole-time-of-flight (TOF) tandem mass spectrometer was operated in both positive and negative ion modes. The optimized parameters were as follows: the ion source temperature 500°C; ion spray voltage 5.5 kV (positive), −4.5 kV (negative); ion sources gas1, ion sources gas2, and curtain gas were set at 50, 50, and 35 psi, respectively. The collision gas and ion spray voltage were set to medium and 4.5 kV, respectively. Data acquisition and analysis were performed using Analyst Software (Sciex, Foster City, CA, USA).

Data preprocessing was performed using Progenesis QI software to obtain m/z values, retention time (RT), and peak intensity. Metabolites were identified by searching and comparing the m/z values, RT, and fragmentation patterns with standards in the Human Metabolome Database (HMDB; http://www.hmdb.ca/) and Metlin database (https://metlin.scripps.edu/). The variable importance of projection (VIP) score of the orthogonal projections to latent structures–discriminant analysis model was used to detect global metabolic changes among comparable groups. Metabolites with a t-test adjusted p-value (p _adj) ≤ 0.05 and VIP ≥ 1 were defined as differential metabolites. The Bioconductor R package from on the Majorbio Cloud Platform (https://cloud.majorbio.com) was used to perform multivariate statistical analysis. An unsupervised method was used to perform a principal component analysis (PCA).

Total RNA was extracted from the frozen leaf samples of each poplar cultivar (QHP, ZSY, and L2025) using the EASYspin Plus Plant RNA Kit (Aidlab, Beijing, China). The extracted RNA was treated with DNase I (TaKaRa) to remove the contaminant DNA. The Dynabeads™ mRNA Purification Kit (Invitrogen) was used for the purification of mRNA from total RNA. The purified mRNAs were cleaved into small fragments and reverse-transcribed into first-strand cDNAs with random hexamer primers. The double-stranded cDNAs were synthesized using the NEBNext Ultra RNA Library Prep Kit (NEB). Adaptors were added to the double-stranded cDNAs, which were further fragmented and enriched by PCR to construct the final cDNA libraries for Illumina paired-end sequencing using an Illumina NovaSeq 600 platform.

After trimming the adapter sequences and removing low quality sequences, the clean reads were mapped to the P. deltoides reference genome (https://www.ncbi.nlm.nih.gov/data-hub/genome/GCA_015852605.2/) and assembled using Hisat2 and StringTie, respectively. The obtained transcripts were annotated using eight databases (Nr, Nt, GO, Swiss-Prot, COG, KOG, Pfam, and KEGG). Gene expression levels were calculated and normalized using fragments per kilobase of transcript per million mapped fragments (FPKM). The DEGs between two groups were analyzed using the DESeq R package. A false discovery rate (p _adj) ≤ 0.05 and |log₂(fold change)| ≥ 1 were used as the criteria to identify significant DEGs (Ali and Bhaskar, 2016). In addition, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of DEGs were performed using the DAVID and KOBAS online tools, respectively. The RNA-seq dataset is available from the Gene Expression Omnibus (accession no. GSE216754).

Pearson’s correlation analysis was conducted by calculating the correlation coefficient between anthocyanin content and the expression of DEGs enriched in the flavonoid–anthocyanin biosynthesis pathway (ko00941 and ko00942). Furthermore, R2R3-MYB and bHLH TFs with differential expression levels were used to perform a correlation analysis with differentially accumulated anthocyanins. Interaction networks between DEGs and differentially accumulated anthocyanins were visualized using Cytoscape 2.8.2 (Cho et al., 2016).

The differential expression patterns of selected structural genes and TFs in the anthocyanin biosynthesis pathway, as detected from the RNA-seq data, were investigated using qRT-PCR (Kodama et al., 2018). The Actin gene was used to normalize the gene expression levels in all groups (Wang et al., 2015). The relative expression level for each gene was calculated using the 2^−ΔΔCt method. Three technical replicates and three biological replicates were analyzed for each sample. All primers used in this study are listed in Supplementary Table 1 .

Heatmap analysis, PCA, and multiple-testing of metabolome and transcriptome data were conducted using R software. Results are expressed as the mean ± standard deviation (SD). The calculated p-value was calibrated by false discovery rate (FDR) correction.

To detect the pigment constituents in purple leaves of QHP and ZSY, the contents of total chlorophyll (Chl-a and Chl-b), carotenoids, and anthocyanins in leaves of L2025, QHP, and ZSY were measured ( Figures 1A, B ). The Chl-a, Chl-b, and carotenoid contents of L2025 leaves (G) were 2.85, 1.54, and 1.50 times those of QHP leaves (P), respectively, and the anthocyanin content of QHP leaves (P) was 14.82 times that of the L2025 leaves (G) ( Figures 1C–E ). The Chl-a, Chl-b, and carotenoid contents of ZSY green leaves (F_G) were 1.70, 1.54, and 1.71 times those of ZSY purple leaves (F_P), whereas the anthocyanin content of ZSY purple leaves was 15.69 times that of ZSY green leaves ( Figures 1C–E and Supplementary Table 2 ). These results indicated that the purple color of the leaves may have been caused by anthocyanin accumulation.

To further explore the anthocyanin components and the corresponding molecular mechanisms of anthocyanin biosynthesis in purple leaves of QHP and ZSY, the same leaf samples, comprising the green leaves of L2025 and ZSY (G and F_G), and the purple leaves of QHP and ZSY (P and F_P), were used for metabolite and transcriptome analysis. Metabolite profiling of the leaves (G, P, F_G, and F_P) was performed using LC-MS/MS analysis. A total of 39 anthocyanin metabolites were identified: 12 cyanidins, 7 pelargonidins, 5 delphindins, 5 procyanidins, 4 malvidins, 2 petunidins, 2 peonidins, 1 pseudopurpurin, and 1 gentisin ( Table 1 ). The PCA showed clear separation of the sample groups. Among the leaf samples, 90.6% of the total variance was explained by PC1 (61.7%), PC2 (16%), and PC3 (12.9%) ( Figure 2A ). The PCA revealed lower variability among the biological replicates. To identify anthocyanin compounds that were differentially accumulated between the purple and green leaves, the anthocyanins from the leaf samples were subjected to comparative analysis according to their relative contents: P vs. G and F_P vs. F_G. The comparison between F_P and F_G showed that 22 anthocyanins increased, whereas 11 decreased ( Figure 2B ). Between P and G, 24 anthocyanins increased, whereas two decreased ( Figure 2C ). In particular, 22 differentially accumulated anthocyanins, including 20 upregulated anthocyanins (7 cyanidins, 5 procyanidins, 2 delphinidins, 2 pelargonidins, 2 peonidins, 1 malvidin, and 1 pseudopurpurin) and two downregulated anthocyanins (cyanidin and malvidin O-hexoside), were detected in both comparisons P vs. G and F_P vs. F_G ( Supplementary Table 3 ).

Correlation heat maps were constructed from the differential anthocyanin metabolites detected from the comparisons P vs. G and F_P vs. F_G ( Figure 3 ). Of 741 pairwise correlations from 39 anthocyanin metabolites, 367 positive (yellow to red) and 95 negative (yellow to blue) correlations (altogether approximately 62.3% of the total) were observed for P vs. G ( Figure 3A ), and 297 positive and 142 negative correlations (altogether approximately 69.2% of the total) were observed for F_P vs. F_G ( Figure 3B ). Significant changes in metabolite–metabolite correlations (both positive and negative) were observed between P vs. G and F_P vs. F_G. This in turn showed that adjustments in the branches of the anthocyanin metabolic pathway occurred during the development of purple leaves in QHP and ZSY. In P vs. G, 31.51% of 367 positive correlations and 47.91% of 95 negative correlations were for procyanidins, and cyanidins–other anthocyanins correlations, respectively. 33.67% of 297 positive correlations and 43.26% of 142 negative correlations were for procyanidins, and cyanidins–other anthocyanins correlations, respectively, that were observed in F_P vs. F_G ( Supplementary Table 4 ). These results indicated that procyanidins and cyanidins were the main contributors to the differences in accumulated anthocyanin metabolites between QHP and ZSY.

To explore the molecular regulatory mechanisms of anthocyanin biosynthesis during the development of purple leaves of QHP and ZSY, 12 cDNA libraries were constructed from the same samples used for metabolite analysis to enable high-throughput RNA-seq ( Supplementary Table 5 ). In total, 535,787,836 raw reads in the range of 42,021,546–46,112,070 and 510,395,564 clean reads in the range of 41,235,222–44,584,822 were obtained for each library ( Supplementary Table 6 ). The Q30 percentage and average GC content percentages were more than 97% and 43.93% for the libraries, respectively. Overall, 80.23%–92.71% of the clean reads were mapped to the genome of P. deltoides for each library. Thus, 32,819, 32,704, 34,047, 32,331, 32,321, 32,337, 34,290, 34,098, 33,214, 33,014, 33,300, and 32,998 genes for each library were identified. These results demonstrated that the RNA-seq data were of high quality and suitable for subsequent analyses.

To identify the DEGs between the purple and green leaves of QHP, L2025, and ZSY, the correlation coefficients (r) between the purple and green leaves for the repeated samples of QHP, L2025, and ZSY were computed based on the FPKM value of all transcripts in each sample. The r values between the biological replicates of the four poplar leaf samples were all greater than 0.9 ( Figure 4A ). The samples were scattered among the four different groups and clustered within the same group, which indicated that the biological replicates could be further used for the detection of DEGs ( Figures 4B, C ). Based on pairwise comparisons with fold change ≥ 1 and FDR ≤ 0.05, 8,217 (2,005 upregulated and 6,212 downregulated) and 11,484 (5,041 upregulated and 6,643 downregulated) DEGs were identified in F_G vs. F_P and P vs. G, respectively ( Figure 4D ).

A GO enrichment analysis was performed to explore the biological functions of DEGs using Blast2GO. A total of 2,345 (F_G vs. F_P) and 3,297 (G vs. P) DEGs were clustered into 83 and 69 functional groups, respectively ( Supplementary Table 7 ). The highest-ranked significant GO terms in the enrichment analysis for the DEGs in F_P vs. F_G and P vs. G were classified into three GO categories: 11 and 11 in the biological process category, 4 and 7 in the cellular component category, and 35 and 32 in the molecular function category, respectively ( Figures 4E, F ). Among biological processes, the most enriched terms were microtubule-based movement, microtubule-based process, and cellular carbohydrate metabolic process. In the cellular component category, the cell wall, external encapsulating structure, and cell periphery were the most represented GO terms. For the molecular function category, the most enriched terms were hydrolase activity acting on acid anhydrides, hydrolase activity hydrolyzing O-glycosyl compounds, and hydrolase activity acting on glycosyl bond catalytic, binding, and transporter activities.

To further explore the DEGs involved in the metabolic pathways, a KEGG pathway analysis was performed. A total of 1,158 (F_P vs. F_G) and 1,818 (P vs. G) DEGs were assigned to 121 KEGG pathways ( Supplementary Table 8 ). The 50 highest-ranked significant KEGG pathways were classified into five categories: metabolism, genetic information processing, cellular processes, environmental information processing, and organismal systems. Among these pathways, the five metabolic pathways with the largest number of DEGs were annotated using the KEGG database in two pairwise comparisons: F_P vs. F_G and P vs. G ( Figures 5A, B ). The five metabolic pathways with the largest number of DEGs from F_P vs. F_G were plant hormone signal transduction (pop04075, 89 genes), phenylpropanoid biosynthesis (pop00940, 78 genes), carbon metabolism (pop01200, 69 genes), starch and sucrose metabolism (pop00500, 61 genes), and plant–pathogen interaction (pop04626, 60 genes). For the comparison of P vs. G, plant–pathogen interaction (pop04626, 113 genes), plant hormone signal transduction (pop04075, 112 genes), phenylpropanoid biosynthesis (pop00940, 92 genes), starch and sucrose metabolism (pop00500, 72 genes), and amino sugar and nucleotide sugar metabolism (pop00520, 65 genes) were the most highly enriched pathways. Notably, plant hormone signal transduction (pop04075), phenylpropanoid biosynthesis (pop00940), starch and sucrose metabolism (pop00500), and plant–pathogen interaction (pop04626) were present among the significantly enriched pathways in both F_P vs. F_G and P vs. G ( Figures 5C, D ).

The DEGs encoding key enzymes involved in anthocyanin biosynthesis were identified based on gene functional annotation and KEGG pathway enrichment analysis ( Figure 6 and Supplementary Table 9 ). As a result, 28 DEGs (four CHS, three FLS, four F3H, two F3’H, two F3’5’H, one DFR, four ANS, and eight UFGT genes) were detected in F_P vs. F_G. Among these genes, three CHS, two FLS, one F3’H, one F3’5’H, one DFR, two ANS, and five UFGT genes were significantly upregulated in F_P vs. F_G. The remaining DEGs were downregulated in F_P vs. F_G. Thirty-three DEGs (seven CHS, three FLS, three F3H, three F3’H, two F3’5’H, one DFR, five ANS, and nine UFGT genes) were detected in P vs. G. Among these genes, five CHS, three FLS, two F3’H, one DFR, one ANS, and four UFGT genes were significantly upregulated in P vs. G. The remaining DEGs were downregulated in the P vs. G group. In addition, 11 DEGs, comprising three CHS (PdeCHS1, PdeCHS2, and PdeCHS7), two FLS (PdeFLS3 and PdeFLS4), one F3’H (PdeF3’H1), one DFR (PdeDFR), one ANS (PdeANS3), and three UFGT genes (PdeUFGT3, PdeUFGT4, and PdeUFGT5), were significantly upregulated in F_P vs. F_G and P vs. G.

It is well known that MYB, bHLH, and WD40 TFs play important roles in regulating anthocyanin biosynthesis. In addition, AP2/ERF, bZIP, WRKY, and Dof TFs are associated with anthocyanin biosynthesis. In this study, seven TF families were detected as differentially expressed in the RNA-seq data ( Figure 7 and Supplementary Table 10 ). A total of 75 MYB, 39 bHLH, 26 WRKY, 52 AP2, 22 Dof, 20 bZIP, and one WD40 TFs were differentially expressed in the comparison group of F_P vs. F_G. Among these TFs, 11 MYB, 5 bHLH, 6 WRKY, 6 AP2, 11 bZIP, and 1 Dof exhibited similar expression patterns as the structural genes in the anthocyanin biosynthesis pathway, which were increased in purple leaves of ZSY. Seventy-six MYB, 34 bHLH, 20 WRKY, 46 AP2, 21 Dof, 1 bZIP, and 1 WD40 TFs exhibited the opposite expression pattern to most structural genes in the anthocyanin biosynthesis pathway, which decreased in purple leaves of ZSY. Moreover, 76 MYB, 37 bHLH, 12 WD40, 30 WRKY, 38 AP2, 20 Dof, and 23 bZIP TFs were differentially expressed in the comparison group of P vs. G. Among these TFs, 19 MYB, 13 bHLH, 6 WD40, 6 WRKY, 6 AP2, and 12 bZIP TFs showed expression patterns consistent with most anthocyanin biosynthetic structural genes, which increased in the leaves of QHP. In contrast, 57 MYB, 24 bHLH, 24 WRKY, 32 AP2, 20 Dof, 11 bZIP, and 6 WD40 TFs exhibited the opposite expression pattern to most anthocyanin biosynthetic structural genes, which increased in the leaves of L2025. To explore the function of differentially expressed MYB TFs, an unrooted phylogenetic tree of MYBs, including the differentially expressed MYBs in P. deltoides and Arabidopsis MYBs, was constructed. Two P. deltoides MYBs (Podel.06G234300 and Podel.04G021100) belonging to the SG5 subgroup of R2R3-MYB TFs were identified ( Figure S1 ).

To verify the expression pattern of DEGs involved in anthocyanin biosynthesis, 12 structural genes encoding key enzymes (three CHS, one F3H, one F3’H, one DFR, one CHI, one ANS, and three UFGT) and two TFs (two SG5 R2R3-MYB TFs) were subjected to qRT-PCR analysis. The qRT-PCR results showed that the expression pattern of all the selected candidate genes was strongly correlated with the RNA-seq data ( Figure 8 and Supplementary Table 11 ).

To further explore the complex regulatory mechanism of anthocyanin biosynthesis in the leaves of colored-leaf poplar, correlation analysis of quantitative changes in anthocyanin content and related DEGs in the leaves of colored-leaf poplar was performed ( Figure 9 ). Based on the detected differentially accumulated metabolites and DEGs involved in anthocyanin biosynthesis, 10 structural genes displayed a significant positive correlation (r > 0.8) with 28 anthocyanins ( Supplementary Table 12 ). The corresponding interaction networks are shown in Figure 9A . The 10 structural genes, comprising three CHS, one F3H, one F3’H, one DFR, one ANS, and three UFGT, played a pivotal role in the anthocyanin biosynthesis pathway. The 28 anthocyanins were cyanidin O-diacetyl-hexoside-O-glyceric acid, cyanin chloride, delphinidin 3-sophoroside-5-rhamnoside, myrtillin chloride, cyanidin 3-O-glucoside, cyanidin O-malonyl-malonylhexoside, cyanidin O-syringic acid, cyanidin chloride, pseudopurpurin, idaein chloride, malvidin O-hexoside, procyanidin B2, procyanidin B3, pelargonin chloride, procyanidin A2, keracyanin chloride, cyanidin O-rutinoside, cyanidin 3-O-glucosyl-malonylglucoside, malvidin 3-galactoside, ferulylpelargonidin di-O-hexosyl-O-pentoside chloride, peonidin, chloride, procyanidin A1, procyanidin A3, pelargonidin O-acetylhexoside, pelargonidin 3-O-malonylhexoside, cyanidin O-acetylhexoside, cyanidin 3-O-malonylhexoside, and peonidin O-hexoside ( Figure 9A ). Correlations between the differentially expressed MYBs, bHLHs, and differentially accumulated metabolites showed that 10 MYB and 10 bHLH TFs exhibited a significant positive correlation with 36 and 30 anthocyanins, respectively (r > 0.8) ( Figures 9B, C ; Supplementary Tables 13 , 14 ). As crucial regulators of plant anthocyanin biosynthesis, the two R2R3-MYB TFs of the SG5 subgroup (Podel.06G234300 and Podel.04G021100), based on phylogenetic analysis and gene function annotation, were significantly positively correlated with six anthocyanins (malvidin 3-galactoside chloride, peonidin O-hexoside, pelargonidin O-acetylhexoside, pelargonidin 3-O-malonylhexoside, cyanidin 3-O-malonylhexoside, and cyanidin chloride) and 16 anthocyanins (procyanidin A1, procyanidin A2, procyanidin A3, procyanidin B2, procyanidin B3, cyanin chloride, cyanidin O-rutinoside, cyanidin O-acetylhexoside, cyanidin 3-O-glucoside, cyanidin O-malonyl-malonylhexoside, keracyanin chloride, ferulylpelargonidin di-O-hexosyl-O-pentoside, pelargonidin 3-O-malonylhexoside, malvidin 3-galactoside chloride, peonidin chloride, and pseudopurpurin), respectively. The strong correlation of Podel.04G021100 with procyanidin and cyanidin indicated that one MYB (Podel.04G021100) may be involved in regulating anthocyanin biosynthesis.