The investigation on metabolites, genes and open chromatins involved in colored leaves of Eucommia ulmoides ‘Ziye’

Eucommia ulmoides Oliver ‘Ziye’ has unique purple-red leaves, which contain a variety of flavonoids, so it has high ornamental and medicinal value. However, the categories of flavonoids and molecular mechanism of specific accumulation of flavonoids in ‘Ziye’ leaves is still unclear. Here, differences in metabolic level, gene expression level, chromatin accessibility and cis-regulatory elements were compared between ‘Ziye’ and ‘Huazhong 1’ with green leaf color by metabolome profiling, RNA-seq, and ATAC-seq. A total of 205 flavonoids were identified from these two varieties using ultraperformance liquid chromatography–mass spectrometry (UPLC-MS). The accumulation of most delphinidin, cyaniding, quercetin, myricetin, and isorhamnetin derivatives peaked in old leaves of ‘Ziye’. Single-molecule long-read sequencing indicated that genes in the phenylpropanoid biosynthesis and flavonoid biosynthetic pathway, as well as many transcription factors including MYB, ERF, and WRKY were highly expressed in ‘Ziye’ leaves. ATAC-seq revealed the presence of cell preferentially enriched peaks, which annotated to 6114 genes. Analysis of the genomic regions preferentially accessible in each cell type identified hundreds of overrepresented TF-binding motifs, highlighting sets of TFs such as MYB, ERF, and WRKY that are probably important for color formation of ‘Ziye’ cell. Interestingly, the TFs within each of these cell type-enriched sets also showed evidence of extensively co-regulating each other. Our work demonstrated how chromatin accessibility and TF expression level influenced the expression of flavonoid biosynthesis associated genes, resulted in flavonoids accumulation in ‘Ziye’ leaf. Our results could lay a foundation for further studies of gene expression and functional genomics in E. ulmoides.

2 medicinal books such as Compendium of Materia Medica and Shennong's Classic of Materia Medica. 47 E. ulmoides 'Ziye' is a variety with excellent purplish red leaf. It is bred on the basis of the natural 48 variation of E. ulmoides with purplish red leaf. Its purple-red leaf color is so stable that it has the 49 potential to further develop new varieties with multi character aggregation. Many studies suggested 50 that much higher flavonoids especially anthocyanins accumulation mainly contributed to purplish red 51 corloration in E. ulmoides 'Ziye' leaf , which resulted in much higher medicinal and ornamental value. 52 However, the components and contents of secondary metabolites, and the molecular mechanism of the 53 gradual purple-red color formation of E. ulmoides 'Ziye' leaf remains to be further studied. 54 It has been found that the flower and leaf colors are mainly determined by the spatially and 55 temporally  Anatomical and physiological analyses of 'Ziye' leaves 112 In order to understand the factor that caused the leaf color formation of 'Ziye', we measured the 113 content of various colored substances. We firstly determined the contents of anthocyanins and total 114 flavonoids in leaf buds, young leaves, adult leaves, and old leaves. We found that both anthocyanins 115 level and flavonoid level increased gradually with leaf development of 'Ziye' and 'Huazhong 1', and 116 both of them were much higher in 'Ziye' leaves ( Figure 1A and 1B). Based on dynamic changes of 117 anthocyanins level and flavonoid level, we selected young leaves (S3 marked in Figure 1A) and old 118 leaves (S6 marked in Figure 1A) which represented for early pigmentation stage and late pigmentation 119 stage for further metabolome, transcriptome and ATAC-seq analysis. P1 and P2 represented for early 120 pigmentation stage and late pigmentation stage of 'Ziye' leaves, meanwhile G1 and G2 of 'Huazhong 1' 121 leaves were corresponding to the P1 and P2, respectively ( Figure 2A). 122 Besides, we compared the contents of chlorophyll, carotene, and soluble sugar between the two 123 varieties. The contents of chlorophyll a and fructose were much higher in 'Huazhong 1'(G2) than those 124 in 'Ziye' (P2), but the difference was not statistically significant ( Figure 1C and 1D). In contrast, the 125 contents of sucrose and glucose were non-significantly higher in P2 than those in G2. Besides, the 126 chlorophyll b and carotene in P2 were similar to G2. 127 The leaf color of 'Ziye' arised from cytoplasm, not the plastid, and the 'Ziye' cytoplasm appeared 128 purplish red under the microscope ( Figure 2B). The chloroplast of 'Ziye' was larger than that of 129 'Huazhong 1' under the TEM ( Figure 2C). In short, anthocyanins and flavonoids enriching in the 130 cytoplasm were the major factor for purplish red color formation in leaves of 'Ziye'. 131 Pigment variation between 'Ziye' and 'Huazhong 1' leaves 132 In order to further clarify which components of anthocyanins and flavonoids were the decisive factors 133 for 'Ziye' leaf color formation, metabolome profile was carried out using UPLC-MS (Table S1) Figure 3A).

138
KEGG enrichment analysis indicated that a large amount of anthocyanins were identified as 139 DEMs, of which, most were highly accumulated in P2, indicating the importance of anthocyanins 140 involved in purplish red color formation of 'Ziye' (Figure 3B and 3C). Among 27 anthocyanins 141 identified in the leaves of two E. ulmoides cultivars ( Figure 3E), all two petunidin derivatives and 142 overwhelming majority delphinidin-and cyanidin derivatives were highly accumulated in P2 of 'Ziye'. 143 In term of peonidin derivatives, 'Huazhong 1' leaves mainly accumulated peonidin-3-O-sambubioside, 144 while petunidin-3-O-rutinoside-5-O-glucoside showed similar accumulation level in two cultivators. 145 Two pelargonidin derivatives were highly accumulated in 'Huazhong 1'. Most flavonols were more 146 abundantly accumulated in 'Ziye', especially in P2, except most kaempferols which were highly 147 accumulated in G2 ( Figure 3D). 148 Overview of single-molecule long-read sequencing results 149 To identify the key genes involved in 'Ziye' leaf color formation, we performed ONT sequencing of 150 'Ziye' and 'Huazhong 1' E. ulmoides leaves at two different growth stages using the PromethION 151 platform. The number of full-length sequences through primer sequence from two ends of reads 152 obtained from each sample vary from 2,161,952 to 2,951,797 (Table S2). After polishing full-length 153 sequences, redundant-remove analysis through minimap2 software, 68,023 transcript sequences were 154 obtained. These isoforms were grouped into three categories: (1) 26,719 known isoforms that mapped 155 to the known gene set, (2) 31,592 new isoforms of annotated genes, and (3) 9,712 novel isoforms 156 belonging to 4,902 likely new candidate genes that were previously unannotated (Table S3).

157
Identification of key pathways involved in 'Ziye' leaf color formation 158 To explore the potential factors that caused the differential coloration and compound content between 159 'Ziye' and 'Huazhong 1' leaves, the KEGG enrichment was conducted. The q-value (corrected p-value) 160 of the richness factor was used to evaluate the importance of the pathways involved in purplish red leaf 161 formation ( Figure 4). 162 For G1 vs P1, most "plant-pathogen interaction", "nitrogen metabolism", and "basal transcription 163 factor" associated genes differentially expressed ( Figure 4A and Table S3). By contrast, the genes 164 involved in "phenylpropanoid biosynthesis" -the precursor for the flavonoids synthesis, "flavonoid 165 biosynthesis", and "cyanamino acid metabolism" were the most significantly enriched, indicating P2 is 166 the key period for the synthesis of flavonoids in 'Ziye' ( Figure 4B and Table S3).

167
Expression analysis of pigmentation-related flavonoid biosynthetic genes 168 Based on the annotation results, a total of 43 differentially expressed functional genes were associated 169 with pigmentation-related flavonoid biosynthesis (Table S4) that could underlie various colors. Two differentially expressed FLS were identified in E. ulmoides. 178 One of them was highly expressed in P1 and the other was highly expressed in P2. Both two 179 differentially expressed F3ʹH were differentially expressed between two color formation stages, but 180 there is no significant difference between the two leaf colors. Among FLS, EUC20656-RA showed high 181 expression abundance in both G2 and P2, while EUC17674-RA was highly accumulated in P1.

183
There were two main branches from dihydrokaempferol into different anthocyanins. Cyanidin and 184 delphinidin formed from dihydroquercetin and dihydromyricetin and catalyzed by DFR and ANS, 185 respectively. Here, one DFR and two ANS genes were differentially expressed, and all these three genes 186 were highly expressed in P1 or P2. Later, peonidin-based anthocyanins were converted from cyanidin 187 through a series of glycosylation and methylation reactions to produce colors. In another anthocyanins 188 biosynthetic pathway, malvidin and petunidin were converted from delphinidin through a series of 189 glycosylation and methylation reactions. The methylation and glycosylation reactions in both two 190 anthocyanins synthesis branches mainly regulated by anthocyanins O-methyltransferase (AOMT). Four 191 AOMTs differentially expressed during 'Ziye' leaf pigmentation, and all of them performed high 192 accumulation level in P2. Most UFGT genes, which played key roles in 3-glucoside formation in 193 different flavonols and anthocyanins biosynthetic branches, showed higher expression levels in P1 and 194 P2.

195
The expression changes of transcription factor in the 'Ziye' pigmentation process 196 Nextly, we focused on differentially expressed transcription Factors (TFs) ( Figure 6). More than half 197 differentially expressed ERF, WRKY, NAC were highly accumulated in P2, suggesting that these TFs To examine the local chromatin landscape of 'Ziye' leaves, and compare the difference between 208 the 'Ziye' leaves and 'Huazhong 1'leaves in the open chromatin region, we performed ATAC-seq. 209 More than 86 million reads were obtained for each biological replicate through paired-end sequencing 210 (Table S5 and  and 17.22% located within 1 kb -2 kb upstream of the TSSs, 8.04% and 9.04% located within 2 kb -3 218 kb upstream of the TSSs, 10.09% (G2) and 10.02% (P2) located within the gene body, 40.96% (G2) 219 and 43.20% (P2) located in the intergenic region ( Figure 7B). This genomic distribution of 220 reproducible THSs suggested that the majority of cis-regulatory regions in the E. ulmoides genome 221 were located in the vicinity of gene core promoters, as previously observed in Arabidopsis (Lu et al., 222 2016). 223 To determine the regulatory landscape of the E. ulmoides genome between 'Ziye' and 'Huazhong 224 1'cells, we generated genome-wide cell type-preferentially enriched maps of accessible chromatin 225 regions (ACRs) from ATAC-seq. In total, 289,300 peaks and 244,953 peaks were identified in the 'Ziye' 226 and 'Huazhong 1' cells, respectively. Among these, 13770 cell preferentially enriched peaks 227 (differentially accessible chromatin regions, dACRs) were found and were annotated to 6114 genes 228 (Table S7). 30.3 % of DEGs identified by transcriptome sequencing had differential ACRs between 229 'Ziye' and 'Huazhong 1'cells. These genes were enriched in the biological process of response to light 230 stimulus, abscisic acid-activated signaling pathway, proteasome-mediated ubiquitin-dependent protein 231 catabolic process, response to stress. These results suggested that opening or closing of nuclear 232 chromatin region of genes related external environmental stimulation especially illumination were 233 essential for the color formation of 'Ziye' leaf ( Figure 8A). The determination of pigment content indicated that the formation of cytoplasmic color is mainly 253 attributed to higher accumulation level of flavonols and anthocyanin. In the process of plant color 254 formation, malvidin is responsible for blue-colored tissues, cyanidin is abundant in reddish-purple 255 (magenta), delphinidin appears blue-reddish and purple colors, petunidin has been detected in purple 256 tissues, pelargonidin for red color, and peonidin for magenta color (Robinson and Robinson, 1932; 257 Tanaka and Ohmiya, 2008). These six types of anthocyanins were all detected in 'Ziye' leaves. Among 258 these, cyaniding, petunidin and delphinidin derivatives were predominant in 'Ziye' leaves, indicating 259 that these anthocyanins may be responsible for the purplish red leaves. Besides, quercetin, 260 isorhamnetin, myricetin which were closely related to the color appearance in petals of lily cultivars 261 and Bauhinia variegate (Zhang et al., 2021). Whereas in our study, malvidin and peonidin derivatives-the 267 two common types of anthocyanins, which caused plant pigmentation, were not the main reason for the 268 purplish red formation of E.ulmoides leaves. In fact, delphinidin, petunidin, cyaniding, quercetin, 269 myricetin, and isorhamnetin derivatives were the main pigments that contributed to the purplish red 270 coloration. in the anthocyanins and flavonoids biosynthesis pathway is 3-glucoside formation by UFGT, which is 295 the key enzyme for anthocyanins and flavonoids stability. The expression of UFGT genes were 296 detected in 'Ziye' cultivars, but not in 'Huazhong 1', which consisted with previous report in 297 Arabidopsis that AtUDT78D2 (flavonoid 3-O-glucosyltransferase) catalyzes the glycosylation of both 298 flavonols and anthocyanins, were associated with anthocyanins and flavonols accumulation in purplish 299 red leaves (Yonekura-Sakakibara et al., 2008). 300 The regulation network involved in leaf color formation of 'Ziye' 301 In order to find key TF families regulating 'Ziye' color formation, we analyzed the differentially 302 enriched chromatin accessible regions to identify putative cell type-specific cis-regulatory motifs as 303 well as the TFs that bind them (Figure 3a). We found that chromatin accessibility, which was positively 304 associated with the expression level of nearby genes in both 'Ziye' and 'Huazhong 1'cells, contributed 305 significantly to the cell type-preferentially enriched expression of DEGs, as 30.3 % of DEGs had 306 differential ACRs between 'Ziye' and 'Huazhong 1'cells, including the key genes involved in flavonoid 307 biosynthesis, phenylpropaniod biosynthesis, and environment stimulation. Thus, the differentially 308 enriched THS between 'Ziye' and 'Huazhong 1' played important roles in regulating the gene 309 expression during purplish red leaf color formation process. 310 Using transcriptome sequencing data for these two cell types, we found many TFs that were 311 highly expressed in 'Ziye' cell type and whose motifs were also overrepresented in 2017). In current study, the importance of MYB and bHLH were also confirmed by our RNA-seq and 316 ATAC-seq data. Interesting, we also found that large amount of ERF and WRKY were highly 317 accumulated in 'Ziye' and many ERF and WRKY binding motifs were enriched in ACRs of 'Ziye' cells.

318
The putative regulatory networks analysis showed three ERFs which showed highest expression level TFs. In short, chromatin opening influenced expression of both functional genes and TFs. 342 In conclusion, many ACRs specially enriched in 'Ziye' meaning the DNA sequence of these ACRs 343 would specially bound by many postfixed with 1% osmic acid for 3 h, and sequentially dehydrated using a graded acetone series (30%, 366 50%, 70%, and 90%). Spurr-embedded tissues were sectioned by an ultramicrotome (LEICA UC6, 367 Weztlar, Germany) and then further observed with a transmission electron microscope (JEM-123O, 368 Tokyo, Japan).

369
Sample preparation and metabolomic analysis 370 Metabolite analysis was performed using an ultra performance liquid chromatography-electrospray 371 ionization-tandem mass spectrometry system (UPLC-ESI-MS/MS). Different leaf samples were rapidly 372 snap-frozen in liquid nitrogen and stored at -80 ℃ until further processing. Leaf samples were 373 freeze-dried by vacuum freeze-dryer (Scientz-100F, Scientz, Ningbo, China), and crushed using a 374 mixer mill (MM 400, Retsch, Haan, Germany), subsequently. 0.1 g powdery tissue was dissolved with 375 extraction buffer, and then placed in a refrigerator at 4 ℃ overnight. The extracting solution was 376 filtrated by centrifugation at 12000 rpm for 10 min before UPLC-MS/MS analysis.

377
The extracting solution was analyzed using an UPLC-ESI-MS/MS system. We selected waters 378 ACQUITY UPLC HSS T3 C18 column (1.8 μm, 2.1 mm*100 mm) as the chromatographic column. 379 Mobile phases A was pure water with 0.1% formic acid, and mobile phases B acetonitrile with 0.1% 380 formic acid, respectively. The column oven was maintained at 40 ℃, and the injection volume was 4 381 μL. Sample measurements were performed with a gradient program that employed the starting 382 conditions of 95% A and 5% B at 10 min, a linear gradient to 5% A and 95% B was kept for 1 min. 383 Subsequently, a composition of 95% A and 5% B was adjusted within 1.1 min and kept for 2.9 min. 384 The flow velocity was set as 0.35 ml per minute. The effluent was alternatively connected to an 385 ESI-triple quadrupole-linear ion trap (QTRAP)-MS. 386 The detection of metabolites eluted from the column was executed by ESI-triple The absolute metabolite was quantified by ANOVA, and Tukey's honest significance difference test 392 was used for paired comparison via PAST v.3. x (Hammer et al., 2001).

393
RNA extraction and quality assessment 394 For each sample, total RNA extraction and separation of residual DNA contamination were excuted 395 using Plant RNA Extraction Kit (TaKaRa MiniBEST, kyoto, Japan), and RNase-free DNase I (Thermo 396 Fisher Scientific, UK), respectively Then, 1.5% agarose gels, and an Agilent Bioanalyzer 2100 system 397 (Agilent Technologies, Santa Clara, CA, USA) were carried out in order to check the purity and 398 integrity of the RNA.

399
ONT library preparation, sequencing and analysis 400 1µg total RNA was prepared for cDNA libraries using cDNA-PCR Sequencing Kit (SQK-PCS109, 401 Oxford Nanopore Technologies, Oxford, UK qRT-PCR assay was performed in three independent biological replicates and three technical replicates. 430 All primers were designed by by Oligo 7 software (Table S8). 431

ATAC-seq 432
The old leaves of 'Ziye' and 'Huazhong 1' were used for ATAC-seq. For each experiment, 50000 cells 433 were centrifuged with 500 g centrifuging parameters for 5 min under 4 ℃, and then the supernatant 434 was removed. The cells were washed with 0.1 M PBS, and then suspended with cold lysis buffer. The 435 supernatant was removed using 4 ℃, 10 min and 500 g centrifuging parameters. The transposing 436 reaction system was configured with the Tn5 Transposase. The cell nucleus was suspended with the 437 transposing reaction system, and the DNA were purified after incubating at 37 ℃ for 30 minutes. The 438 PCR reaction system was configured with the purified DNA, and then the PCR amplification was 439 performed. The DNA libraries were sequenced on Illumina platform, finally. 440 Illumina Sequenced Reads Processing Raw Reads were filtered to remove adapters (the adapters 441 and reads less than 35bp) and low quality reads, and then high-quality Clean Reads provided in FASTQ 442 format were obtained for subsequent analysis. The Bowtie2 software was used to compare the 443 high-quality Reads with the reference genome to obtain the alignment efficiency of the sample Reads functional elements of each gene on the genome, the differential Peak region was annotated with the 458 functional databases mentioned above. 459 The differentially accessible chromatin regions (dACRs) that significant enriched in 'Huazhong 1' 460 cell or 'Ziye' cell were considered as cell specific type, and were further used for motif analysis. The 461 cell type-enriched Peak region from each cell type were first adjusted to the same size (200 bp). The 462 sequences present in these scaled regions were isolated using the Regulatory Sequence Analysis Tools 463 (RSAT), which also masks any repeat sequences (Medina-Rivera et al., 2015). The masked sequences 464 were run through MEME-ChIP with default parameters to identify motifs that were present in higher 465 proportions than expected by chance (E-value ≤ 0.05) (Machanick and Bailey, 2011).

466
Accession numbers 467 We deposited the RNA-seq data and ATAC-seq in the Sequence Read Archive (SRA) database 468 with an accession number SRP218063.

CONFLICT OF INTEREST 470
The authors declare that they have no conflicts of interest associated with this work. This work was supported by the National Science Foundation of China (32201643).

SUPPORTING INFORMATION 478
Additional Supporting Information may be found in the online version of this article. 479 Table S1 The identified flavonoid and their accumulation level in E. ulmoides leaves. 480 Table S2 full-length sequences statistics. 481 Table S3 Annotation of all expressed genes in E. ulmoides leaves. 482 Table S4 The expression profiles of all identified E. ulmoides genes. 483 Table S5 The evaluation statistics of ATAC sequencing data. 484 Table S6 Overview of mapping of the ATAC sequencing read 485 Table S7 The cell preferentially enriched peaks.  486  Table S8 List of primers used for qRT-PCR analysis. 487 Figure 1. Comparison of pigment contents between 'Ziye' E. ulmoides and 'Huazhong 1'E. ulmoides. 488 The ratio of (A) total anthocyanins, and (B) total flavonoids to fresh weight during leaf growth and 489 development. The ratio of (C) Chlorophyll and carotene, and (D) soluble sugar to fresh weight in old 490 leaves of 'Ziye' and 'Huazhong 1'. 491