Fresh buds of FLJ and rFLJ were collected in Spring 2013 from the Gardens of Yate Co., in Shandong Province, China. Samples were authenticated by Prof. Liying Yu (Guangxi Medical Plant Garden, China) and stored in liquid nitrogen. FLJ plants were cultivated in outdoor flowerpots near laboratory, and watered regularly to keep the soil moist. The impact of temperature and humidity on perceived air quality was studied in the ranges 20–30°C and 40–60% RH. Leaves of living FLJ plants were wiped gently with 80 μM 5-aza by cotton swabs for 24 h (about 12 h light/12 h dark in one diurnal time), while untreated leaves were used as control group. Treated and untreated leaves were from the same plants, and the experiment was repeated three times.

Genomic DNA was isolated by the cetyl trimethyl ammonium bromide (CTAB) method as described by Doyle (1991). DNA was quantified using a NanoDrop 2000 Spectrophotometer (Thermo Scientific, San Jose, CA, United States) and stored at 4°C until further analysis.

Bisulfite sequencing primers were designed based on the 5′-UTR sequences using MethPrimer (Li and Dahiya, 2002). The optimum product size (<200 bp), Tm (<55°C), and primer size (<25 bp) were taken into consideration for primer design. Primers were designed to flank the CG-rich regions of the promoter sequence. All primers for bisulfite sequencing were designed avoiding CpGs in the sequence. Conventional primers were also designed using Primer3¹ for amplification of the respective 5′-UTR regions of untreated DNA that served as the master sequence for bisulfite sequence analysis. Bisulfite DNA conversion was performed using 1 μg of genomic DNA and an EpiTech Bisulfite Kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol. PCR was performed using primers located outside the target region and designed for single-strand methylation detection. PCR products were cloned into pCR2.1 using a TA cloning kit (Invitrogen, Carlsbad, CA, United States). For each genotype, at least 20 independent clones were sequenced using the M13R primer, and data were analyzed by Cymate (Hetzl et al., 2007).

The methylation levels of selected genes were detected as described by Shinozuka et al. (2015). Briefly, 100 ng of genomic DNA was subjected to two separate treatments: methylation-sensitive digestion and mock digestion. For the methylation-sensitive digestion, the methylated DNA was digested by the methylation-sensitive enzyme MspJI, while the remaining un-methylated DNA was detected by RT-PCR. For the mock digestion, inactive MspJI enzymes were added to the reaction and the products represented the total amount of input DNA. The relative amount of each DNA fraction (methylated and un-methylated) was calculated by using the standard ΔCt method (Ordway et al., 2006), normalizing the amount of DNA in each digestion against the total amount of input DNA in the mock digestion. The amount of methylated DNA was defined as 2^-Ct(Mo-Ms)/2^-CtMo, where 2^-CtMo is the total amount of input DNA, and 2^{-Ct Ms} is the amount of un-methylated DNA.

Genomic DNA was divided into two parallel groups for MspJI digestion or mock digestion. The MspJI digestion solution included 100 ng of genomic DNA, 10 units of MspJI (New England Biolabs, Ipswich, MA, United States), 1 × NEBuffer 4 (New England Biolabs, Ipswich, MA, United States), 0.5 μM reaction activator, and deionized water to a final volume of 20 μL. The mock digestion solution included inactive MspJI enzymes. Reactions were incubated at 37°C for 4 h, and followed by heat inactivation at 65°C for 15 min.

The RT-PCR solution included 10 μL 2 × I-5^TM High-Fidelity Master Mix (MLAB, Los Angeles, CA, United States), 0.4 μL 50 × ROX (Applied Biosystems, Foster City, CA, United States), 0.4 μL forward and reverse primers (10 μM), 1 μL 20 × EvaGreen (Biotium, Hayward, CA, United States), and 2 μL reaction temple (endonuclease-treated samples diluted 1:10 with double-deionized water). Primer sequences are listed in Supplementary Table S1. RT-PCR conditions were as follows: 2 min at 50°C, 3 min at 95°C, followed by 40 cycles of 10 s at 95°C and 20 s at 58°C. Melting curves were used to verify the specificity of PCR products by increasing the temperature from 60°C to 95°C in sequential steps of 0.2°C for 15 s. Amplification was performed in triplicate, and the no-template control was performed for each gene. The relative abundance of genes was determined using the comparative Ct method in ABI 7500 2.0.1 (Applied Biosystems, Foster City, CA, United States).

The current assembly of FLJ genomic sequences was investigated for the 5′-UTR region of LjPAL1 (JX068601), LjPAL2 (JX068602), LjPAL3 (JX068603), Lj4CL1 (JX068604), Lj4CL2 (JX068605), LjC4H1 (JX068606), LjC4H2 (JX068607), and LjHQT (ACZ52698) using the Blastn algorithm (Altschul et al., 1997). An e-value cut-off of 10^-30 was applied for homolog recognition. The obtained 5′-UTR regions were subjected to regulatory sequence search analysis using PlantCARE and Softberry to identify the cis-regulatory elements in each sequence (Lescot et al., 2002).

We applied RNA-seq to FLJ and rFLJ and generated over 100 million reads using the Illumina GA platform (Yuan et al., 2012). The current assembly of FLJ transcriptome was investigated for bZIP homologs in Arabidopsis thaliana using the BLASTx algorithm (Altschul et al., 1997). An e-value cut-off of 10^-5 was applied for homolog recognition. The rFLJ transcriptome was searched for orthologs from FLJ genes. An all-against-all sequence comparison was performed with BLAST (cut-off < 10^-20) to identify orthologs, which were determined from the best reciprocal hits (80% alignment length) as described by Tatusov et al. (1997). All retrieved sequences from Genscan² were used for gene prediction (Burge and Karlin, 1998). The predicted gene models were further examined and corrected by manually comparing them with related genes of other plant species. The functional and structural domains were predicted by InterProScan (Zdobnov and Apweiler, 2001) and Blast2GO (Conesa et al., 2005), respectively.

The deduced amino acid sequences were adjusted manually by BioEdit 7.0.0³ using the default parameters. The ORF of LjbZIPs was identified by BioEdit (Hall, 1999). The phylogenetic analysis of alignments was performed by ClustalW (Thompson et al., 2002) and MEGA 6.0 (Tamura et al., 2013) using neighbor-joining analysis. The reliability of tree topologies was evaluated using bootstrap support with 1,000 replicates (Pattengale et al., 2010)^. The sequences of 57 Arabidopsis bZIPs were downloaded from the Arabidopsis Information Resource⁴.

LjbZIP11 (Group A), rLjbZIP4 (Group A), LjbZIP8 (Group C), rLjbZIP18 (Group C), LjbZIP10 (Group D), and rLjbZIP1(Group D) cDNAs were amplified by PCR and cloned into pGEM T-easy cloning vector. After confirming inserted sequences, the cDNAs were double-digested with BamHI and SalI, and inserted into the pGEX-4T-1 expression vector that had been digested with BamHI and SalI. The recombinant plasmid pGEX-LjbZIP11, pGEX-rLjbZIP4, pGEX-LjbZIP8, pGEX-rLjbZIP18, pGEX-LjbZIP10, and pGEX-rLjbZIP1were transformed into Escherichia coli BL21 (DE3) chemically competent cells (Beijing TransGen Biotech, Beijing, China). Expression of the recombinant protein was induced with 0.4 mM IPTG at 16°C overnight, and cells were resuspended in PBS buffer (137 mmol/L NaCl, 2.7 mmol/L KCl, 10 mmol/L Na₂HPO₄, and 2 mmol/L KH₂PO₄) and disrupted by sonication. The recombinant protein were purified by affinity chromatography by using ProteinIso GST resin (TransGen Biotech, Beijing, China) following the manufacturer’s instruction. The purified lysate was centrifuged at 12,000 × g for 30 min at 4°C, and the supernatant was loaded onto a 8% SDS–PAGE gel after denaturation with SDS loading dye at 100°C for 5 min. The gel was stained with Coomassie Brilliant Blue G-250 and decolorized with a destaining solution (70% Ultrapure water; 10% ethanol; 20% acetic acid) (Zha et al., 2016).

The G-box sequence in the promoter sequence of LjPAL2 (GenBank: JX068602) was obtained using Softberry. The oligonucleotides (5′-CTCTCCTGTCCACGT GTCATCAAC-3′ for the G-box sequence) were synthesized and labeled with biotin (Sangon Biotech, Shanghai, China) using the Lightshift Chemiluminescent EMSA kit (20148, Pierce, Thermo Fisher Scientific, United States). Complementary labeled strands were mixed together in an equimolar ratio and annealed at 25°C after denaturation at 90°C. Gel mobility shift assays were performed by incubating 40 fmol of labeled probe with 15 fmol LjbZIP proteins (LjbZIP11, LjbZIP8, LjbZIP10, rLjbZIP1, rLjbZIP4, and rLjbZIP18) and 4 pmol (or 0 pmol) competing oligonucleotides in the binding reactions at room temperature for 20 min. Mixtures were size fractionated in a non-denaturing 30% polyacrylamide gel followed by drying, transferred to nitrocellulose membranes, and detected by streptavidin-HRP/chemiluminescence for biotin-labeled probes.

The LjbZIP8 sequence was ligated into the pE3025 vector, which was digested with EcoRI and KpnI to generate pGEM-LjbZIP8 plasmids. LjbZIP-GFP was under the control of the CaMV 35S promoter. The construct was confirmed by sequencing and used for the transient particle-bombardment transformation of onion (Allium cepa) epidermis cells using a gene gun (Bio-Rad, Hercules, CA, United States). Following 24 h of post-bombardment incubation, GFP fluorescence in transformed onion cells was observed under a confocal microscope (Zeiss, Oberkochen, Germany).

To determine the transactivation activity, the ORF of LjbZIP8 was generated by PCR amplification and cloned into the EcoRI/SalI digested pGBKT7 vector to generate pBD-LjbZIP8 plasmids. The constructs were transformed into YGR2 cells by the lithium acetate-mediated method (Yuan et al., 2013). The transformed yeast strains were placed on SD/–Trp medium at 28°C for 2 days. Yeast transformants were then transferred and streaked onto solid SD/–Trp/–His/–Ade to score the growth response after 3 days. The pBD-GAL4 and pGBKT7 vectors were used as positive and negative controls, respectively.

LjbZIP8 was inserted into the binary vector pCambia1305 to produce p35Spro-LjbZIP8 plasmids that were then used to transform Agrobacterium tumefaciens EHA105. Tobacco (Nicotiana tabacum) leaf disks were transformed via an A. tumefaciens mediated leaf disk procedure (Horsch et al., 1989) and selected using 50 mg/L Hygromycin B and 200 mg/L carbenicillin. After rooting and acclimatization, the regenerated plants were cultivated in a greenhouse to set seeds by self-pollination. Transgenic plants were used for further analysis.

Total RNA was reverse-transcribed using Reverse Transcriptase MMLV (Takara, Dalian, China). PCR was performed in triplicate using an ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, United States) with SYBR Premix Ex Taq kit (TaKaRa, Dalian, China), according the manufacturer’s instructions. Gene-specific primers of NtPAL1 (M84466), NtPAL2 (D17467), NtPAL4 (EU883670.1), and LjbZIP8 (KT218632) were designed using Primer3 (Supplementary Table S1). The length of PCR products ranged from 100 to 250 bp. Ntactin was chosen as an endogenous control for studying gene expression in various samples of transgenic tobacco. The specificity of amplification was assessed by melting curve analysis, and the relative abundance of genes was determined using the comparative Ct method as suggested by ABI 7500 2.0.1 (Applied Biosystems, Foster City, CA, United States). Significance in relative expression levels of genes were determined by a paired-samples t-test. P-values below 0.05 were considered statistically significant for all tests.

HPLC-grade acetonitrile and methanol were obtained from Thermo Fisher Scientific (Waltham, MA, United States). Phosphoric acid was purchased from Aladdin Chemistry (Shanghai, China). Reference standards (purity > 98%), including neochlorogenic acid, CGAs, cryptochlorogenic acid, isochlorogenic acid B, isochlorogenic acid A, and isochlorogenic acid C were obtained from Tongtian Biological Technology (Shanghai, China). Deionized water was produced by Milli-Q system (Bedford, MA, United States).

Approximately 0.5 g pulverized plant samples were weighed and mixed with 50 ml of 50% MeOH aqueous solution (v/v) in a flask. The mixture was extracted in an ultrasonic bath (KQ-250DB, Kunshan Ultrasonic Instrument, Kunshan, China) for 30 min. The extraction solution was centrifuged at 4,000 rpm for 10 min, and the supernatant was diluted to 50 ml with 50% MeOH aqueous solution (v/v) and stored at 4°C. After filtration through a 0.22 μm PTFE membrane (Jinteng Technologies, Tianjin, China), the extract was ready for UPLC analysis.

Chromatographic analysis was performed using an Acquity UPLC I-Class system (Waters, Milford, MA, United States) equipped with a HSS T3 column (100 mm × 1.0 mm, particle size 1.8 mm; Waters, Milford, MA, United States). The flow rate was 0.4 mL/min, and the mobile phase consisted of 0.2% phosphoric acid in water (solution A) and acetonitrile (solution B). The elution scheme was as follows: 0–2 min, 5–10% B; 2–8 min, 10–30% B; 8–11 min, 30–45% B; 11–14 min, 45–80% B; and 14–17 min, 5% B to equilibrate the column for the next injection. The column temperature was maintained at 40°C and the injection volume was 1.0 μL. Photodiode array detector spectra were measured over a wavelength range of 210–400 nm, and chromatographic data were collected and manipulated using Empower^TM 2 (Waters, Milford, MA, United States).

All of the CGAs were quantified using chromatograms extracted at 340 nm (Supplementary Figure S8). The content of CGAs constituents were routinely quantified using standard curves of corresponding compounds. The calibration curves were constructed using seven different concentrations (0.1, 0.5, 10, 50, 100, 500, and 1000 μL/ml) of each standard and by plotting the concentration of the standard against the peak area. The areas of overlapped peaks were calculated by UPLC.

To investigate whether the content of CGAs were different between FLJ and rFLJ, we performed analysis on different tissues by UPLC. It was observed that FLJ had lower content of neochlorogenic acid, CGAs, cryptochlorogenic acid, isochlorogenic acid B, isochlorogenic acid A, and isochlorogenic acid C compared with those in rFLJ (Table 1).

Many genes (PAL1, PAL2, PAL3, 4CL1, 4CL2, C4H1, C4H2, and HQT) are involved in CGAs biosynthesis. Between FLJ and rFLJ, SNPs were identified in the genes involved in CGAs biosynthesis, but these did not result in changes in the amino acid sequences of these key enzymes (unpublished data). To acquire CpG island regions and candidate CpG loci, the genomic sequences of the 5′-UTR of genes involved in CGAs biosynthesis were obtained from the FLJ genome database (unpublished) and submitted to MethPrimer⁵. The results showed that 24 and 20 candidate CpG loci were located in the 452 bp 5′-UTR region of LjPAL2 and LjC4H1, respectively (Supplementary Table S2 and Figure S1). However, the methylation status investigated by bisulfite sequencing PCR (BSP) indicated that DNA methylation only occurred in the 5′-UTR of LjPAL2. The CpG methylation ratio of PAL2 in rFLJ (80.45%) was significantly different compared with that of rFLJ (61.24%) (Supplementary Table S3). Moreover, the expression of PAL2 orthologs in rFLJ buds was 3.25-fold higher than that in FLJ buds (Supplementary Figure S2), which was consistent with the results in our previous study (Yuan et al., 2012).

Due to DNA methylation occurs in the upstream region of LjPAL2, we predicted the transcription factor binding sites and regulatory elements using Softberry⁶, G-box regulatory elements were found mainly in the 5′-UTR region of LjPAL2. Further analysis of the transcription factors bound to the G-box sequence identified the bZIP family (Supplementary Figure S3 and Table S4). Therefore, we speculated that bZIP transcription factors might bind to the G-box of LjPAL2. To identify bZIP proteins in FLJ, a preliminary BLASTX search was performed against the transcriptome database. Only hits with E-values below e^-30 were considered as members of this gene family. A total of 11 LjbZIP proteins and 24 rLjbZIP proteins were identified to have conserved bZIP domains, and the 11 LjbZIP sequences were submitted to GenBank with the accession numbers KT218625–KT218635 (Supplementary Table S5).

Based on sequence similarity, the identified LjbZIP sequences were clustered into 10 subgroups, according to clades with at least 80% bootstrap support (Figure 1). The validity of the phylogenetic reconstruction was confirmed by the fact that the subgroups were the same as those observed in previously constructed phylogenetic trees (Jakoby et al., 2002). LjbZIP1 was classified into subgroup S, which includes homologs that are transcriptionally activated by stress (e.g., cold, drought, and anaerobic conditions, and wounding) (Kusano et al., 1995) or are expressed in specific parts of the flower (Strathmann et al., 2001). LjbZIP2 and LjbZIP3 were classified into subgroup H, while LjbZIP8 belongs to subgroup C. LjbZIP4 and LjbZIP9 were classified into subgroup F along with At4g35040, At2g16770, and At3g51960. LjbZIP5 was classified into subgroup G, which includes GBFs from A. thaliana and their parsley homologs CPRF1, CPRF3, CPRF4a, and CPRF5 that are mainly linked to ultraviolet and blue-light signal transduction and the regulation of light-responsive promoters (Kircher et al., 1998). LjbZIP7 and LjbZIP11 were classified into subgroup A, which includes genes that play roles in ABA or stress signaling (Lescot et al., 2002). LjbZIP10 was classified into subgroup D, which includes genes that participate in defense against pathogens (Xiang et al., 1997) and control the number of floral organs (Chuang et al., 1999). In general, the gene functions of a clade are not conserved across plant species; thus, the gene function facilitates the confirmation of paralog and ortholog relationships.

To determine whether LjbZIP and rLjbZIP proteins can bind to the G-box element of the LjPAL2 Promoter, several bZIP proteins including LjbZIP11 (Group A), rLjbZIP4 (Group A), LjbZIP8 (Group C), rLjbZIP18 (Group C), LjbZIP10 (Group D), and rLjbZIP1 (Group D) were selected for EMSA analysis (Supplementary Table S6). All of the six bZIP proteins were produced in E. coli BL21 (DE3) cells, which were subsequently disrupted by sonication (Supplementary Figure S4). The G-box element in the LjPAL2 5′-UTR sequence was predicted as CACGTG using Softberry and used as a probe for EMSA analysis (Figure 2A).

Electrophoretic mobility shift assay assays showed that purified LjbZIP8 bound specifically to the G-box element and unlabelled probes inhibited the binding (Figure 2B). However, rLjbZIP1, rLjbZIP4, rLjbZIP18, LjbZIP10, and LjbZIP18 showed no binding with the G-Box element (Supplementary Figure S5). Further control experiments showed that no binding bands were detected with crude proteins of E. coli BL21 (DE3) cells with or without the empty vector. These results clearly demonstrated that LjbZIP8 could bind to the G-box element of the LjPAL2 5′-UTR.

Bioinformatic analysis revealed that LjbZIP8 contains a 432-bp ORF that encodes a 144 amino-acid protein. The predicted conserved domains showed that LjbZIP8 belongs to the bZIP GBF1 subfamily with an N-terminal proline-rich domain (Figure 3A). Multiple sequence alignments of the deduced LjbZIP8 amino acid sequence revealed that LjbZIP8 has homology with the corresponding enzyme from other plants that has a plant G-box binding 1 factor domain, particularly in the conserved domain at the N-terminal (Figure 3B).

To verify the intracellular localization of LjbZIP8, the full-length cDNA sequence of LjbZIP8 was fused in the correct ORF in front of the 5′-terminus of the GFP reporter gene under the control of the CaMV 35S promoter. Onion (Allium cepa) epidermal cells were transformed with the LjbZIP8-GFP construct or a construct containing GFP alone by particle bombardment. GFP, without an N-terminal fusion, was localized throughout the cell, while LjbZIP8-GFP was accumulated mainly in the nucleus, suggesting that LjbZIP8 is a nucleus-localized protein (Figure 4).

Most nuclei-localized proteins function as transcription factors. To investigate whether LjbZIP8 had transcriptional activity, we performed transcription activity analysis using the yeast GAL4 system. The full-length cDNA of LjbZIP8 was fused to the GAL4 DNA binding domain of the pGBKT7 vector to construct pBD-LjbZIP8, which was then used to transform the yeast strain YGR2. Yeasts transformants with pBD-GAL4 and pBD-LjbZIP8 could grow on selection media lacking tryptophan (SD/–Trp) or tryptophan, histidine, and adenine (SD/–Trp/–His/–Ade; Figure 5), while those with the empty vector pGBKT7 could not grow on the selection media, suggesting that LjbZIP8 functions as a transcriptional factor.

To further analyze the role of LjbZIP8 in vivo, tobacco was transformed with LjbZIP8, and the integration of LjbZIP8 into the tobacco genome was confirmed using PCR (Supplementary Figure S6). Real-time RT-PCR showed that the expression of LjbZIP8 was significantly increased in the transgenic plants. Three independent transgenic lines (9, 11, and 17) overexpressing LjbZIP8 were selected for further analysis. To investigate whether the overexpression of LjbZIP8 in transgenic tobacco plants affected the accumulation of CGAs, we performed analysis on transgenic leaf samples by HPLC. The content of neochlorogenic acid, CGAs, and cryptochlorogenic acid were lower in transgenic plants overexpressing LjbZIP8 compared with those in wild-type plants and transgenic plants transformed with the empty vector pCambia1305 (Supplementary Table S7).

We further investigated whether the overexpression of LjbZIP8 affected the expression of NtPALs. The expression levels of NtPAL1, NtPAL2, and NtPAL4 were decreased in transgenic tobacco plants compared with wild-type plants (Figure 6), indicating that LjbZIP8 might regulate the synthesis of CGAs by inhibiting the expression of NtPALs (Figure 7).

The DNA methylation inhibitor 5-aza was applied to the leaves of FLJ at a concentration of 80 μM and its effects on the transcriptional level of LjPAL2 and the content of CGAs were evaluated after 1 day. The results showed that LjPAL2 expression level of 5-aza-treated leaves was 5.13-fold lower than that in control leaves (Supplementary Figure S7). Furthermore, UPLC analysis showed that the content of neochlorogenic acid, CGAs, cryptochlorogenic acid, isochlorogenic acid A, and isochlorogenic acid C were lower in 5-aza-treated leaves than those in control leaves (Table 2), suggesting that lower DNA methylation rates inhibit the expression of LjPAL2 and accumulation of CGAs in FLJ (Figure 7).