Cytosine Methylation Changes the Preferred Cis-regulatory Configuration of Arabidopsis WUSCHEL-Related Homeobox 14

1. Introduction

The homeobox (HB) superfamily, defined by a conserved homeodomain (HD) of 60–66 amino acids, is widely distributed across plants and animals [2]. HB superfamily members serve as key regulators of growth and development. Plant HB superfamily TFs are involved in various processes, such as embryo patterning [3], organ development (root, shoot, floral meristems) [4], vascular development, and stress responses [5,6,7,8]. Although the HB superfamily is shared in eukaryotes, the HB TFs belonging to the WUSCHEL (WOX) family are plant-specific. WOX TFs are phylogenetically categorized into the modern, intermediate, and ancient clades based on evolutional evidences and species distributions [9]. The Arabidopsis TF — AtWOX14 (WUSCHEL-related homeobox 14) is a member of the ancient clade. AtWOX14 plays critical roles in diverse developmental processes, including bud regeneration in callus [10], lateral root and stamen development [11], vascular cell differentiation [12], lignification [13], and the accumulation of bioactive gibberellins in the stem of inflorescence [14].

To regulate the expression of target genes, transcription factors (TFs) selectively bind to specific DNA sequences. This binding process can interact with chromatin accessibility [15], histone modification [16], TF-TF interaction (cooperative or competitive) [17], as well as DNA methylation [18,19]. DNA methylation typically inhibits TF binding, for example, plant WRKY-family TFs show reduced DNA-binding affinities upon cytosine methylation [20]. However, a few TFs do enhance their binding affinity to methylated DNA [21,22,23]. For example, a comprehensive survey of 542 TFs reported that TFs critical for early embryonic development (e.g., POU5F1 and HOXB13) preferentially bind methylated DNA, highlighting a dual role of methylation in modulating TF-DNA interactions [21]. TFs frequently bind DNA as dimers, for instance, ARF forms homodimers upon binding [24], and bZIP10 and bZIP53 can heterodimerize [17]. While the previous studies systematically addressed the impact of methylation on nucleotide specificity of TFs [21], there are limited insights into how 5mC affects the configuration of dimeric binding of TFs.

Given the critical roles of AtWOX14 in plant growth and development, this study aims to thoroughly clarify its DNA-binding specificity and regulatory targets, and how these are influenced by DNA-methylation. In the absence and presence of methylation, altogether we curated 4 monomeric and 63 dimeric binding models for AtWOX14, and validated the ER0/1 dimeric binding model with EMSA. Methylation has globally changed the binding specificity of AtWOX14, not only the sequence specificity, but notably, also the configurational specificity (the spacing between the monomers) of dimerization. The specificity changes upon methylation also dictate the binding strength of AtWOX14 to its genomic targets — either being enhanced or inhibited dependent on the sequence of the cis-regulatory elements (CREs) of AtWOX14, suggesting that the methylation effect on AtWOX14’s regulatory activity can be complex.

For AtWOX14, here we use SELEX to profile the biochemical specificity, and DAP-seq to identify the genome-wide regulatory targets. This is because the randomness and complexity of SELEX input library ensure unbiased and sensitive discovery of a TF’s binding models [25,26], while DAP-seq uses genomic DNA to directly interrogate the binding sites and target genes of the TF [23]. These two methods complement each other and are preferably combined. Also, in contrast to animals whereby the majority of methylation happens on the CG dinucleotides, in plants methylation can occur in all of the CG, CHG, and CHH contexts [27]. Taking this into consideration, in this study we replaced all Cs in the SELEX/DAP-seq library with 5mCs when examining the methylation effects, using a recently developed protocol [1].

2. Results

2.1. Methylation Inhibits Monomeric But Enhances Dimeric Binding of AtWOX14

To unravel the original DNA-binding specificity of AtWOX14, we first performed SELEX without methylation (Figure 1A). Recombinant AtWOX14 protein is incubated with a 101-bp dsDNA library containing a randomized region. The ligands with high affinities to AtWOX14 were then purified together with AtWOX14, and amplified for the next cycle. To examine the impact of methylation on the specificity of AtWOX14. The SELEX process is also modified to incorporate 5mCs instead of Cs during the PCR amplification [1], generating a fully methylated DNA library in each cycle (Figure 1A). After 5 cycles, the enriched DNA pools were sequenced and analyzed.

We first confirmed that SELEX has enriched binding signals of AtWOX14. First, specific 10-mer subsequences were enriched over SELEX cycles (Figure 1B). Such enrichment is exponential because a straight-line is observed after taking the logarithm (Figure S1). Enrichment-based mutual information (E-MI) analysis [28,29] has also validated the presence of TF-binding signals in the enriched library (Figure 1C).

We next investigated the impact of methylation on the preferred kmers of AtWOX14 (Figure 1D). Interestingly, when we compare all 8-mer subsequences between methylated and normal libraries, the most abundant 8-mers are more enriched in the SELEX library (Figure 1D, left). In contrast, when all 12-mer subsequences were compared, the most abundant 12-mers are enriched in the Methyl-SELEX library (Figure 1D, right). We further performed motif discovery with 8-mers and 12-mers as seeds (using Autoseed [30]), and found that 8-mers give monomeric AtWOX14 motifs, while 12-mers give dimeric motifs. Therefore, it can be concluded that methylation has inhibited the monomeric binding of AtWOX14, while enhancing its dimeric binding.

2.2. Methylation Changes the Preferred Homodimeric Configuration of AtWOX14

While methylation is known to affect TFs’ nucleotide specificity and affinity [21,23], how the incorporation of 5mC affects TFs’ dimeric configuration is yet to be explored. We noticed that dependent on methylation, the most enriched dimeric motifs are with different configurations regarding the spacings of the half-sites (Figure 1D, right). Therefore, we next de novo identified all dimeric (also monomeric) binding motifs of WOX14 (Figure 2A, B) with the reported method [30], and compared their enrichments between the SELEX and the Methyl-SELEX libraries.

The identified dimeric motifs covered all possible relative orientations of monomers (direct repeats: DR, inverted repeats: IR, or everted repeats: ER). By defining the half-site as “YAATYA”, we are able to calculate the spacing between the half-sites, and name each dimeric motif by combining the relative orientation and the spacing (name labeled to the left of each motif, Figure 2A, B).

From the normal SELEX library, we identified 34 motif models, including 2 monomeric and 32 homodimeric motifs (Figure 1C). In general, the dimeric motifs with narrower spacings displayed greater enrichments and higher information contents (IC) (Figure 1C). The ER1 motif exhibited the highest IC and enrichment, followed by ER0, DR3-5, and IR6-8 (Figure 1C). An important signature of the functional TF motifs is that they are enriched around transcription start sites (TSSs). Accordingly, we found that the identified high-IC dimeric motifs all enriched upstream of TSSs (Figure 1D). However, we noticed that the dimeric models show higher enrichments at TSSs than the monomeric models. This consists with the stronger binding signals of dimeric models in SELEX (Figure 1C), and potentially also because the higher transactivating activity of the dimeric CREs are more strongly selected for during the evolutional history. To date, the reported motifs (e.g., in public databases [31,32,33]) of WOX TFs are mainly monomeric. Here we show that the dimeric motifs are more enriched in SELEX and around the TSSs, revealing that WOX14 has a strong preference to bind as homodimers.

In the presence of methylation, we identified 2 monomeric and 29 dimeric motifs, including 8 IR, 13 DR, and 8 ER motifs (Figure 2D). Consistent with the SELEX results (Figure 1C), dimeric motifs show higher enrichments and ICs compared to the monomeric motifs, further supporting the preference of WOX14 to bind as a dimer. The high-IC motifs are also enriched around TSSs (Figure 2D). Interestingly, the most enriched dimeric motif has switched from ER0 to ER1 (Figure 2B), the motif-based enrichment analysis here is also consistent with the 12-mer based analysis (Figure 1D, right), confirming the methylation effect on the homodimeric configuration of AtWOX14.

2.3. The Targets of ER0 and ER1 Are of Different Physiological Functions

We next focused on the high-affinity dimeric configurations ER0 and ER1 to investigate their specificity changes to DNA methylation. Positional comparison revealed that methylation affected the affinity of cytosines at positions 1, 5, 6, 7, and 11 (Figure 3A). Interestingly, the affinity of cytosines is not universally reduced upon methylation. Whereas the cytosines at pos1 of ER0/1 and pos11 of ER1 show decreased affinities when methylated (most preferred base changed from G to A, or C to T), the methylation effects on cytosines at the boundaries of the half-sites (pos5/6 of ER0, corresponding to pos5/7 of ER1) are more complex (Figure 3A) — their affinity increases in ER0, but decreases in ER1. This suggests that the effects of methylation on cytosine affinity can depend on the spacing of the TF dimer, which was also previously unrecognized. Similar to the dimeric motifs, the effect of methylation on cytosines of the monomeric motifs is also positional-dependent.

We further validated the dimeric binding of ER0 and ER1 to WOX14 by electrophoretic mobility shift assays (EMSA). Recombinant WOX14 protein is incubated with DNA ligands containing ER0/ER1 consensus, and also the consensus of the monomeric motif M1 as a control (Figure 3B). The results show that all ligands can bind WOX14 monomer; however, only ER0 and ER1 ligands can bind dimers of WOX14 at high protein concentrations (Figure 3B). This indicates that ER0 and ER1 are genuine homodimeric CREs of WOX14, and also suggests that the protein-level dimerization of AtWOX14 is weak, a proper CRE element on DNA is necessary to facilitate its dimerization. This could have explained why the dimeric configuration of AtWOX14 is sensitive to DNA methylation.

To explore the potential biological roles of ER0 and ER1 CREs, we conducted Gene Ontology (GO) analysis of their target genes (Figure 3C) using motif matches in the promoter region. We found that there is almost no overlap between the biological processes regulated by ER0 and ER1. Considering only the reported functions for AtWOX14 [10,11,12,14], genes downstream of ER0 were enriched in biological processes such as cell wall organization, cell wall modification, lignin metabolism and biosynthesis, gibberellin metabolism, and meristem initiation (Figure 3C), which were associated with vascular cell lignification. In contrast, genes associated with ER1 were predominantly involved in cell fate determination and auxin transport. Therefore, DNA methylation could facilitate WOX14 to utilize different dimeric binding modes (ER0 vs. ER1) and regulate distinct physiological processes.

2.4. Methylation Effects on Genome-Wide Binding Targets of AtWOX14

To explore the impact of methylation on genomic binding sites of AtWOX14, we generated native (DAP-seq), unmethylated (ampDAP-seq), and fully methylated (Methyl-ampDAP-seq) DAP-seq libraries. Similarly to Methyl-SELEX, Methyl-ampDAP was performed by amplifying genomic DNA with 5mCs instead of Cs during PCR. The Venn diagram (Figure 4A) revealed that the libraries of DAP and ampDAP are more similar to each other, while Methyl-ampDAP has identified much more (12,138) unique peaks. The similarity between DAP and ampDAP is potentially explained by the relative scarcity of CHH methylation in the genome.

Motif discovery of the DAP-seq libraries using Homer [34] has confirmed that the peak sites of all libraries contained the monomer motif of WOX14 (Figure 4B). These motifs are consistent with the motifs derived in SELEX (Figure 2A, B). From the DAP libraries, we also confirmed that ER0 CREs were preferentially bound under methylated conditions (Figure 4C-E), because: (1) there are more peaks containing ER0 in Methyl-ampDAP than DAP/ampDAP libraries (Figure 4C); (2) the intensity of ER0 peaks is also higher than ER1 peaks in the Methyl-ampDAP but not in the DAP/ampDAP libraries (Figure 4C). In contrast to ER1, while the number of peaks containing ER0 is similar across all the DAP libraries (Figure 4C), the intensity of ER0 peaks is weaker than ER1 peaks only in the ampDAP library (Figure 4D). We next examined the individual binding sites of AtWOX14, and found that, the ER0 site near AT4G35410 showed the highest occupancy upon methylation (Figure 4C, left), whereas the ER1 site upstream of AT1G60545 showed the highest occupancy in the absence of methylation (Figure 4C, right). The ER0 site within the promoter of AT1G15640 harbors a sequence almost the same as the consensus of ER0 motif (Figure 4F), we also validated with EMSA for its capability to bind AtWOX14 dimer (Figure 4D).

Overall, the analysis of DAP-seq libraries also confirmed that upon methylation, AtWOX14 has switched its preferred dimeric mode from ER1 to ER0.

3. Discussion

This study performed SELEX and DAP-seq for AtWOX14 both in the presence and absence of cytosine methylation, and systematically studied how methylation affects the specificity and genome-wide binding sites of AtWOX14. The most interesting finding is that the preferred dimeric configuration has switched from ER0 to ER1 upon DNA methylation (Figure 4E). ER0 and ER1 were also found to regulate target genes with different functions. Therefore, it is reasonable to hypothesize that when genome-wide methylation level is low, AtWOX14 will regulate gibberellin synthesis and accumulation, thereby promoting the lignification of vascular cell [14]. However, as genome-wide methylation level increases, for example, when the level of methyl-transferase DRM increases or when the RdDM pathway is activated [35], AtWOX14 could rebalance between the ER0 and ER1 sites and emphasize more on the ER0 targets, which are involved in cell differentiation, and cell fate determination processes [10,11] (Figure 4E). Consequently, the methylation level may serve as a switch to toggle between different downstream regulatory networks of AtWOX14.

This mechanism can potentially be extrapolated to other TFs that bind as homodimers but not protein-level homodimers. This is because in general, protein-level dimerization of TFs depends on amino acids outside the DNA-binding domain (DBD) [36], which tends to be far away from the DNA and thus unlikely affected by DNA methylations [37]. However, if the protein-level interaction between TF monomers is weak, then the homodimeric binding can be heavily mediated by the DNA scaffold. In this case, it is conceivable that the chemical modifications on DNA will impose more prominent effects on the preferred dimeric configurations. Notably, dimeric binding mediated by DNA scaffold is not a rare case but represents a prevalent phenomenon, because a high-throughput study has revealed that many TFs only show a small interaction surface between the two monomers when docked to the high-affinity dimeric CRE sequence [38]. The solved TF-DNA structures also suggest that TF dimers can form purely dependent on DNA, with two monomers spaced at >10 Å (e.g., PDB: 1HJB) [39]. Dimeric binding without direct contacts can be induced by the allosteric effect [40], or if a TF alters DNA shape upon binding.

When addressing the methylation effect on a TF’s binding affinity and specificity, the previous work mainly relied on comparing between the DAP and ampDAP libraries. This approach artfully made use of the sporadic 5mCs present in the native genome, to interrogate how methylation changes the affinity of each C in a known motif model of the TF [20]. However, this method suffers from the inability to de novo build binding models in the presence of 5mC. Also, because the non-CG methylation sites on the genome are relatively sparse [41], one can easily overfit the methylation effect with limited data points. In this work, we adopted the recently developed Methyl-SELEX method [1]to prepare input DNA ligands with 100% 5mC, and extrapolated the protocol to set up Methyl-ampDAP. These methods will largely facilitate our understanding on the methylation effects of TF binding.

Consistent with previous researches [21,42], we found that the methylation effect on a cytosine’s affinity to AtWOX14 is positional dependent (Figure 3A). But interestingly, we also found that 5mC overall represses DNA binding of AtWOX14 monomers but facilitates the binding of its dimers (Figure 1D). Thus, the methylation effect is also dependent on the binding mode of a TF.

In summary, this work comprehensively analyzed the influence of DNA methylation on DNA-binding of AtWOX14. It is also interesting to explore whether the mechanistic insights revealed by this work — DNA methylation can rewire the regulatory network by changing the preferred homodimeric configuration of a TF — can be extrapolated to other TFs.

4. Materials and Methods

4.1. Cloning and Protein Expression

The full-length CDS sequence of WOX14 was cloned into the bacterial expression vector pETG20A-His using the In-Fusion Kit (Vazyme, C112). This resulted in the generation of an N-terminal His-tagged fusion protein construct of WOX14. The construct was then transformed into Rosetta 2 (DE3) pLysS strains and induced overnight at 17°C to obtain the recombinant WOX14 protein. After overnight incubation, the bacteria were collected by centrifugation and resuspended in 150 µL of Buffer A (50 mM Tris, 300 mM NaCl, 10 mM imidazole, pH 7.5). The bacterial suspension was then subjected to freeze at -80°C overnight. The frozen bacteria were lysed at room temperature for 2 h using lysis buffer (0.5 mg/ml lysozyme, 2 mg/ml DNase I, 1 mM PMSF). Following lysis, His-tag Ni Sepharose (Sangon, C600332) was added and incubated with the lysate for 1 h. The protein-bead complexes were washed twice with Buffer A and Buffer B (50mM Tris, 300 mM NaCl, 50 mM Imidazole, pH 7.5) using a 10 kD ultra-centrifugal filter (Sigma, UFC8010) to remove contaminants. Finally, the recombinant WOX14 protein was eluted from the beads using Buffer C (50mM Tris, 300 mM NaCl, 500 mM Imidazole, pH 7.5).

4.2. SELEX

SELEX (Systematic Evolution of Ligands by Exponential Enrichment) was performed following previously described protocols [43]. The input SELEX library was designed to include a central random region of 11 bp and flanking regions that were compatible with Illumina sequencing (TruSeq adapter sequences). 6 µl fusion proteins were incubated with 6 µl of SELEX library (30 ng/µl) and 30 µl of TCAPT buffer (140 mM KCl, 5 mM NaCl, 1 mM MgCl2, 3 µM ZnSO4, 100 µM EGTA, 10 mM Tris, pH 8, 0.1% Tween) at 25 °C for 1 h. The DNA-protein complexes were enriched with His-Tag Beads (Cytiva, 17371222), followed by washing away unbound DNA ligands using the HydroSpeed plate washer (Tecan, 30190101). Subsequently, the bound DNA ligands were amplified by PCR (for Methyl-SELEX, PCR using the 5-methyl-dCTP solution) with amplification primers

(forward primer 5′-CCCTACACGACGCTCTTCC-3′,

reverse primer 5′-CAGACGTGTGCTCTTCCG-3′), and the resulting PCR products were used as the input library for the next cycle. After 5 cycles, the enriched DNA ligands were amplified with the PE primers

(forward primer5′-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTAC ACGACGCTCTTCC-3′,

reverse primer 5′-CAAGCAGAAGACGGCATACGAGAT-barcode-GTGACTGGA GTTCAGACGTGTGCTCTTCCG-3′) and Illumina sequenced.

4.3. Methyl-SELEX

Methyl-SELEX process was derived from SELEX with the addition of a DNA methylation step [1]. Specifically, the input SELEX library was amplified by PCR using 5-methyl-dCTP to generate the methylated library which was then incubated with the protein. This methylated library was incubated with 6 µl recombinant WOX14 protein and 30 µl of TCAPT buffer at 25 °C for 1 h. Methylated DNA ligands bound to WOX14 were separated from the unbound sequences through a washing step, and subsequently amplified again by PCR using 5-methyl-dCTP. The bound DNA ligands amplified by PCR were sequenced for analysis and used as input for the next cycle.

4.4. Data Analyses of SELEX

The raw sequencing data generated from the Illumina NovaSeq 6000 platform was first demultiplexed according to the unique i7 index of each sample. In the pre-processing stage, the low-quality reads and PCR duplicates were removed, the adaptors were trimmed, and then the paired-end reads were merged. To de novo discover the binding motifs of WOX14 from the pre-processed SELEX reads, the Autoseed algorithm was used with parameters “-40N <Background sequence> <Signal sequence> 1 8 10 0.35 - 40 100” [30]. The previously described Enrichment-based Mutual Information (E-MI) analysis is used to evaluate the signal strength of transcription factors in the SELEX library [28,29]. The information content (IC) of motifs also was calculated. To calculate motif enrichments, motif matching analysis was performed using motifmatchr [44]. Gene Ontology (GO) analyses were conducted using the R package clusterProfiler [45].

4.5. Electrophoretic Mobility Shift Assays

The protein-DNA binding buffer was prepared with 10 mM Tris, 50 mM NaCl, 1 mM MgCl2, 4% glycerol, and 0.5 mM EDTA. Double-stranded DNA probes, either monomeric or dimeric, were synthesized and diluted to a concentration of 50nM for use in EMSA. The probe sequences used were as follows:

M1: ATGCTAGCTCCATCTGTATTGATTGTTTATGGCGGTGACGTACT;

ER0: ATGCTAGCTCCATCTGTGATTGCAATCAATGGCGGTGACGTACT;

ER1: ATGCTAGCTCCATCTGTGATTGACAATCAATGGCGGTGACGTACT;

ER0-AT1G15640: AATAAGATGTACCCAATGATTGCAATGATGAGCCCAATGA GTGC.

The purified WOX14 protein was then diluted into various gradients and incubated with 30 ng of the DNA probe in a 20 µl protein-DNA binding buffer at 25°C for 40 minutes. The resulting reaction mixture was loaded onto a 6% native PAGE gel and run in 0.5× TBE buffer at 4°C for 90 minutes at 110 V. The gel was subsequently stained with Gel-Blue (UElandy, S2019L) and imaged using a Bio-Rad scanner (Bio-Rad, 1708195EDU).

4.6. DAP and Data Pre-Processing

The ampDAP, methyl-ampDAP, and DAP protocols followed the previously described workflow [46]. To construct the DAP library, genomic DNA (gDNA) was extracted from 7-day-old Arabidopsis seeding and fragmented using the Covaris S2 sonicator. The fragmented DNA (~200bp) was then purified with DNA clean beads (Vazyme, N411) and ligated with adapters using ligation kits (Vazyme, N203, N204). For the input library of amp-DAP, the gDNA was further amplified with specific primers to obtain completely unmethylated gDNA. For the input library of methyl-amp-DAP, the gDNA was amplified with specific primers and 5-methyl-dCTP to obtain the fully methylated gDNA. Next, 6 µl of His-tagged WOX14 protein and ~200 ng of the input DAP library were incubated within30 µl of TCAPT buffer for 1 h. The protein-DNA complexes were pulled down using His-Tag Beads and washed with a HydroSpeed plate washer (Tecan, 30190101). The gDNA bound to WOX14 was amplified by PCR and sequenced on an Illumina NovaSeq 6000. For methyl-ampDAP, PCR amplification of gDNA using the 5-methyl-dCTP solution.

Raw reads were demultiplexed, and sequences with low quality or adapters were removed before aligning to the TAIR10 version of the Arabidopsis genome (https:// www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas.gz) using BWA v.0.7.17 [47]. Genome annotations were obtained from the Araport11 files (https://www.arabidopsis.org/download/list?dir=Gene s%2FAraport11_genome_release). Peak calling was performed with MACS3v.3.0.0a7 [48]. Motif discovery was performed using Homer(http://homer.ucsd.edu/homer/) [34].