A high-throughput approach to identify high-resolution TF footprints genome-wide

To define specific candidate TF-binding sites within accessible chromatin regions, we developed MOA-seq to capture the putative footprints of native DNA-protein interactions. The assay was streamlined to be scalable for high-throughput application and includes a computational pipeline to improve the discovery of putative TF-binding sites as summarized in Fig 1. The protocol (S1 File) starts with the preservation of DNA-protein interactions by formaldehyde-crosslinking prior to tissue homogenization and nuclei extraction. To recover these small interaction regions, we took advantage of the endo- and the exo-nuclease activities of MNase, both of which are inhibited at sites of protein-bound DNA. After MNase treatment, the small DNA fragments can be recovered by reverse crosslinking, protein digestion, and size selection (S1 Fig). To increase the recovery of small (~30–80 bp) TF-bound fragments, while efficiently removing larger nucleosome-size fragments, we added higher salt concentration during the decrosslinking and performed the size selection after adapter ligation during library construction (Fig 1, Steps 1–3). Following sequencing and read mapping to a reference genome (summarized in S1 Table), we plotted the density of aligned fragment midpoints to determine MOA footprints (MFs, average 29.5 bp) and used these to improve the spatial resolution of putative TF-binding event prediction (Fig 1, Step 4). We then performed de-novo motif discovery (Fig 1, Step 5) to annotate potential cis-elements and compare them to previously defined TF motifs in plants. In summary, this MOA-seq protocol was designed to repurpose MNase from mapping nucleosomes to mapping smaller particles within ACRs, resulting in a simple, scaleable, and antibody-free approach to globaly identify putative TF-bound cis-elements at relatively high spatial resolution.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Overview of MOA-seq method.

Flowchart summarizes the five key steps that allow for genome-wide, high-throughput and -resolution transcription factor footprinting and candidate motif discovery. Crosslinked, purified nuclei are subjected to light digestion with MNase (Step 1), producing fragments that range from large nucleosome-sized (grey dsDNA) to small, protein-bound (green dsDNA) DNA-footprints. High salt is included during DNA purification and adapter ligation (Step 2) is performed prior to size selection (Step 3). These steps serve to increase the recovery of desired small footprint fragments while efficiently eliminating unwanted nucleosome-sized fragments, reducing the required sequencing depth. After MOA-seq fragments are aligned to the genome (Step 4) as illustrated around the na1 gene [63]. The resolution of the putative protein/DNA interaction sides can be further enhanced to approximately 30 bp wide footprints by defining fragment centers (MOA-footprints, orange). Sequences underlying MOA-footprints were used as input for de novo motif discovery (Step 5), helping to define the maize earshoot cistrome.

https://doi.org/10.1371/journal.pgen.1009689.g001

To assess the reproducibility of the method we analyzed the bioreplicates, summarized in Fig 2. At the genic level, MOA coverage and MFs showed good agreement between replicates (Fig 2A). This agreement was quantitatively supported genome-wide by multiple different measures of reproducibility and data quality. First, we found that ~70% of bases were shared within significantly enriched peaks that account for only ~0.6% of the genome for each replicate (Fig 2B). Second, the Fraction of Reads in Peaks (FRiP) score was 8.6% for MOA coverage, considerably higher than the 1% threshold recommended by the ENCODE consortium [26].Third, genome-wide coverage correlations were determined because they provide a statistical measure of replicate and control correspondences, independent of peak calling (Fig 2C). We found that the Pearson’s r between bio-replicates (r = 0.96) far exceeded the ENCODE recommendations for correlation values exceeding r = 0.8 [26,27]. In contrast the correlations between the MOA replicates and control were quite low (r < 0.25). Based on all of these measures, we concluded that the MOA-seq methodology as performed and described here was decidedly reproducible.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Reproducibility of MOA-seq coverage and footprints from bio-replicate analysis.

Comparison of the two bio-replicates, rep1 and rep2, each from ABCD combined (Methods). (A) MOA-seq read coverage (Mc, light orange), MOA footprint coverage (Mfc, dark orange), MNase control coverage (CTRL-c), and MNAse control footprint coverage (CTRL-fc) are shown for two representative areas around the sdg118 gene (left, Zm00001eb363830) and a DUF-like gene (right, Zm00001eb364010). For each example, the peaks for the combined datasets (Mp and Mfp) are shown. (B) Venn diagram of rep1 and rep2 shows the 1-to-1 quantification of overlapping bases within overlapping peaks. The total shared base pairs (7.42 Mbp) represent 69.7% of Rep1 and 67.8% of Rep2 bases in their respective peaks. (C) Biological replicate correlation analysis for replicate 1 (rep1), replicate 2 (rep2) and the MNase control (CTRL).

https://doi.org/10.1371/journal.pgen.1009689.g002

Next we assessed the degree to which MNase sequence-specific bias contributes potential false positive MFs by including a partially MNase-digested chromatin-free DNA input control (Fig 2A and 2C, CTRL; S2 Fig), which did not explain nor match the MOA-seq coverage profiles. Taken together, these findings suggest a good reproducibility of MOA-seq bioreplicates and consistency with previous reports that the known AT-rich bias of MNase does not substantially bias MNase-based nucleosome footprint assays [28].

MOA-seq identifies shared and unique footprints

Transcription factors preferentially bind directly upstream and downstream of the transcription start sites (TSSs) and transcription termination sites (TTSs), respectively, while being depleted in protein-coding regions. To test whether MFs showed a similar pattern, we analyzed their genome-wide distribution relative to genes, as shown in Fig 3. In our analysis, almost half (47%) of all MFs flanked genes. The second-largest fraction of MFs (31%) were located in intergenic regions, while the smallest fraction (3%) overlapped with protein-coding regions (Fig 3A). The highest enrichment of MFs was observed directly upstream and downstream the TSS and TTS, respectively (Fig 3B). Overall this pattern reflects typical TF-binding sites and open chromatin patterns around genes and intergenic enhancer regions. We previously showed that MNase-based differential nuclease sensitivity profiling could be used to map functional regions of the maize genome [29]. Because both the DNS-seq and MOA-seq make use of light-digest MNase fragments, albeit in different ways, to define accessible chromatin [30], we expected that the two would mark similar genomic regions. Indeed, by inspection and peak overlap analysis, the two methods were congruent (S3 Fig), and confirmed genome-wide with shared peaks accounting for 59% of all peaks from MOA-seq and 67% of all peaks from DNS-seq.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. MOA-seq identifies shared and unique regions in open chromatin.

Comparison of MOA-seq to open chromatin profiling techniques DNS-Seq and ATAC-Seq. (A) Genome-wide distribution of MOA-seq footprints relative to genomic annotations. (B) MOA-seq read coverage (CPM) is non-randomly distributed around genes. Genes were scaled (metagene) to 2kb and average coverage is shown within the metagene body and adjacent genomic regions +/- 1.5kb from the nearest metagene boundary (grey dash). (C) Comparative analysis, shown as bar graphs, of MOA-seq peaks overlapping (shared, blue, %) or non-overlapping (unique, yellow) with ATAC-seq (MOA/ATAC), and the reciprocal intersections (ATAC/MOA). (D) and (E) Distribution of previously published DNS-seq (Positive light-heavy digest = blue, negative = red), ATAC-Seq (grey) coverage, MOA-seq (light orange) coverage (Mc) and peaks (Mp), and MF (dark orange) coverage (Mfc) and peaks (Mfp), surrounding example genes, vp14 (D) and dwf4/brs1 (E), with direction of transcription and the TSS indicated (arrow).

https://doi.org/10.1371/journal.pgen.1009689.g003

Since MOA-seq and ATAC-seq peaks also share a similar distribution pattern, but deploy different nucleases, we analyzed the number of shared sites between them. For this, we compared our data with recently published ATAC-seq data from similar-stage earshoots [31,32], the best available dataset for comparison. While the majority (76%) of ATAC-seq peaks overlapped with one or more MOA-seq peaks, only a minority (35%) of MOA-seq peaks overlapped with one or more ATAC-seq peaks (Fig 3C). Although ATAC and MOA peaks collectively capture a similar fraction of the genome (ATAC, 0.52%; MOA-seq, 0.55%), MOA-seq did identify more, smaller, and unique peaks compared to ATAC-seq (S2 Table). However, given the differences in peak sizes and numbers, some of these differences could be due to differences in the statistical cutoffs between the segmentation algorithms. Reanalyzing the ATAC-seq data with iSeg, we found similar overlaps. The higher number of peaks called for MOA-seq does not necessarily indicate, therefore, a higher performance.

The differences in the distribution of MOA and ATAC-seq peaks were observed as distributed over large regions covering hundreds of thousands of bases as well as at the genic level (S4 Fig). The latter is illustrated, e.g., for the well-characterized abscisic acid biosynthetic pathway gene viviparous14 (Fig 3D, vp14) and the brassinosteroid biosynthesis gene dwarf4 (Fig 3E, dwf4). Many MOA-seq and ATAC peaks were shared in the vicinity of these two genes. However, MOA-seq identified additional peaks, not detected by ATAC, which overlapped with accessible chromatin regions also detected by DNS-seq in earshoot [7] as well as TF-binding sites such as teosinte branched1 (tb1) and knotted1 (kn1) (Fig 3D and 3E). These findings illustrate that MOA-seq and ATAC-seq map ACRs globally with considerable agreement, especially around promoter regions, whereas MOA-seq often further resolves subregions or identifies additional regions (S5 Fig), similar to that described for Arabidopsis [12].

MOA-seq maps to functional elements

Given the expectation that MOA-seq identifies regulatory regions enriched for cis-elements, we examined MFs at known and predicted TF-binding sites, as shown in Fig 4. We first analyzed the MOA-seq average coverage at ChIP-seq sites previously published for the TF encoded by fasciated ear4 (fea4) [33] which is active in both ear and tassel. Although tassel FEA4 ChIP-seq data was used, MOA-seq coverage peaked at FEA4 binding sites, compared to the surrounding area (Fig 4A). In the case of the KNOTTED1 (KN1) [34] transcription factor, both tassel and ear ChIP-seq data was available. Similar to FEA4, MOA-seq coverage peaked at KN1 binding sites present in both tassel and ear as well as sites unique to the ear. The tissue-specific binding of KN1 was mirrored by the MOA-seq coverage (Fig 4B). In addition to TF-binding sites, we examined evolutionarily-defined conserved noncoding sequences (CNS) which are enriched for cis-regulatory elements [35]. We found that MFs aligned with CNSs, even at intergenic sites far removed from known promoters (Fig 4C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. MOA-seq overlaps with TF-binding sites and correlates at promoters with mRNA levels.

Average earshoot MOA-seq coverage was plotted around the midpoints for FEA4 tassel ChIP-seq peaks (A), Knotted1 ChIP-seq peaks subdivided into Ear only, Ear & Tassel, or Tassel only (B), and conserved non-coding sequences (CNS) showing all, position-shuffled control (Shuffled), or an intergenic subset classified as being more than 1 kb from genes (C). For each plot, the average peak size was used to define a single metapeak size (metapeak, indicated in each plot) for fitting the data within the peak region; metapeak flanking regions are in real bp relative to the metapeak region. (D) MOA-seq coverage is positively correlated with gene regulation. Mean MOA-seq coverage for genes divided into quintiles (1^st highest mRNA to 5^th lowest mRNA level) according to mRNA levels in earshot tissue.

https://doi.org/10.1371/journal.pgen.1009689.g004

Previous studies demonstrated that gene expression levels and chromatin accessibility show positive correlation at proximal promoters [29,30]. To examine this relationship for MOA-seq in a genome-wide manner, we sorted the 36,441 maize genes into quintiles based on their steady-state mRNA levels from matched earshoot tissues, and inspected MOA-seq profiles around the TSSs. Across expression quintiles, we detected a positive correlation between average gene expression levels and MOA-seq read coverage, most clearly evident within ~300 bp upstream of the TSSs (Fig 4D). Together, these analyses establish compelling evidence for our starting hypothesis that MOA-seq footprints can define regulatory loci likely to be occupied by DNA-binding proteins.

MOA-seq footprints identify hundreds of putative TF motifs

Global chromatin structure assays such as ATAC-seq, DNAse-seq, and DNS-seq identify accessible chromatin regions, but their larger average footprint profiles may reduce the accuracy of footprint analysis [36]. In contrast, given the small average size of 29.5 bp for MFs (S2 Table), and their presence at known functional cis-elements, the MOA-seq-defined footprints should be ideal for de novo motif and putative TF-binding site discovery, all of which aims to elucidate genome-wide native cistromes.

We used RSAT (http://rsat.sb-roscoff.fr/) for de novo discovery of enriched motifs within the 143,009 MFs (S2 Table, iSeg-BC7 for B73v3) as input genomic sequences. We identified a total 215 significantly enriched motifs (S3 Table), summarized in the matrix similarity tree shown in Fig 5 and separated for the two types of motif discovery methods used: the 6–7 bp long "oligos" (Fig 5A, n = 140) or the "dyads" (Fig 5B, n = 75). The number of individual sites for individual MOA-identified motifs ranged from less than 200 to more than 6,000. We compared the 215 MF motifs with those previously identified both in-vivo and in-vitro by ChIP-Seq and DAP-Seq, respectively, using multiple plant motif databases. Among the full collection of motifs, 85% (119/140) of the oligo motifs and 75% (56/75) the dyad motifs showed similarities to motifs listed in at least two motif databases (S3 Table, including footprintdb [37], and softberrydb [38]). To characterize their distribution relative to genes, we split them into groups of those found within repetitive DNA, or not. The frequency of repeat-overlapping motifs ranged from 3% for dym05 to 88% for dym43, with a median value 15% for the 215 motifs. Plotting their positional tendencies relative to genes, we found that the median distance to TSS for motifs not within repeats showed a remarkable clustering within 100 bp proximal to the promoters (Fig 5C). Higher resolution TSS annotations [39] may better refine the positional locations of these motifs. For each of the 215 motifs, we produced summary catalog files for the dyad motifs (S6 File) and the oligos motifs (S7 File). These files provide reference documents, one page per motif family, listing the total number of sites, percentage found in annotated repeats, the sequence LOGOs, positional frequency distributions around TSSs, local base frequencies, and average local MOA coverage.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Global summary of maize earshot motifs discovered by MOA-seq.

(A-B) Radial similarity-cluster tree of the 140 oligo (A) and 75 dyad (B) consensus motifs identified within MOA-seq footprints. Circular sectors underlying tree branches are colored to highlight the frequency at which each motif occurred (high > 2000, medium 2000–1000, low < 1000). The inner (red dash) colored ring indicates similarity threshold for motif clustering at 90% similarity, delineating different motif clusters (alternating light and dark gray ring) flanked by motif names and consensus sequence logos at the inner and outer circles, respectively. (C) Median distance to the nearest transcription start site (TSS) plotted for each motif found in regions annotated as within repeats (IR, grey) or in non-repetitive regions (NR, red).

https://doi.org/10.1371/journal.pgen.1009689.g005

MOA-seq footprints improve the spatial resolution of TF binding event prediction

Depending on the size of the sheared DNA, ChIP-seq data usually identifies regions substantially larger than the underlying TF-binding site(s) even for TFs in multimeric complexes. Given the relatively small size of MFs, we tested whether the motifs identified at MFs were enriched at known TF-binding sites which also have the canonical motif. We found that 93% (3882/4195) and 82% (2770/3396) of the earshoot-expressed FEA4 and TB1 TF ChIP-seq peaks, respectively, overlapped with at least one MF. We then tested whether those overlapping MFs were enriched for MOA-defined motifs that closely resemble the canonical motif for the respective TF (summarized in Fig 6). The oligo motif om015 (Fig 6A) closely resembles the database motif for the bzip TF family as well as the published FEA4 consensus sequence, NCGTCA [10]. Genome-wide, om015 and a similarly abundant control motif om006 were enriched 6.9- and 1.2-fold, respectively, at FEA4 binding sites within 3 Kb gene promoters (Fig 6B, prom.). Similar enrichment trends were also observed in the entire mappable genome (Fig 6B, gen.). Examples of this overlap can be seen within FEA4 ChIP-seq binding sites upstream of putative target genes zw18-like (Fig 6C) and bx9 (Fig 6D). Similarly, a database search of the dyad motif dym33 revealed that part of its consensus sequence closely resembled the consensus motifs of the TCP TF family as well as that of one of its members, TB1, GGNCCC [40] (Fig 6E). Genome-wide, dym033 was enriched 11.7-fold at TB1 ChIP-binding sites in 3kB gene promoters, compared to only a 1.3-fold enrichment of a more abundant dym36 MOA-motif (Fig 6F). Examples of this overlap can be seen within TB1 ChIP-seq binding sites upstream of putative target genes nactf49 (Fig 6G) and downstream of d8 (Fig 6H). For the two examples (FEA4-like om015, TB1-like dym33), we found that the matching MOA motif families showed the expected enrichment at the corresponding TF ChIP peaks (S6 Fig). Notably, of the 215 motif families, om015 was ranked 2nd for enrichment at FEA4 ChIP peaks, whereas dym33 was ranked 3rd for enrichment at TB1 ChIP summits. Similarly, among the top 5 motifs enriched at KN1, the motif family om050 includes two matches to the core TGCA motif associated with KN1 binding sites [34]. More broadly comparing MOA peaks to the binding sites of 104 TFs identified by ChIP-seq we found that 66% and 69% of the base pairs in MOA coverage and footprint peaks, respectively, were shared with those under the TF ChIP-seq peaks, compared to only 4% for peaks from the MNase control (Fig 6I). In addition to this global agreement, we observed genic level overlap illustrated by the TF vs. MF profiles around the tip1 gene (Fig 6J). Given that members of TF families often share similar motifs, and that multiple members of a TF family can be co-expressed, there may exist some discrepancies between the consensus sequence of motif families defined from our global assay with MOA-footprints compared to those from individual TF ChIP assays.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. MOA-seq footprints improve the spatial resolution of TF binding event prediction.

(A-F) ChIP-seq binding peaks of (A-B) FEA4 or (E-F) TB1, ATAC-Seq coverage (Ac) and peaks (Ap), MOA-seq footprint coverage (Mfc) and peaks (Mfp), and MOA-seq motif om015 in the promoter region of example genes (A) zw18-like (GRMZM2G039505), (B) bx9 (GRMZM2G161335), (E) nactf49 (GRMZM2G347043), and (F) d8 (GRMZM2G144744). (C, G) MOA-seq motif (C) om015 (G) dym33, aligned to their best database hit (footprintdb) of bzip and TCP family, respectively. The asterisk for (G) dym33 indicates that only a portion of the larger dym33 logo is shown. (D, H) Enrichment of motifs (D) om015 and a control motif om006 within FEA4 ChIP-seq peaks, or (H) dym33 and control motif dym36 with TB1 ChIP-seq peaks, relative to the mappable B73 genome (gen.), and the 3 kb, non-ORF overlapping, region upstream of B73 gene models (prom.). (I) Genome-wide overlap analysis of base pairs shared between MOA (MOA, MFs) and ChIP-seq of 104 TFs [27] or the MNase control (CTRL). (J) Representative example of the complex landscape of TF and MFs profiles surrounding tip1 (Zm00001eb00373). TF tracks shown: MYB, TF55 rep1 GSM4095564; C2C2 TF235 rep2 GSM4095532,; bHLH, TF6 rep2 GSM4095589; TCP, TF237 rep1 GSM4095534.

https://doi.org/10.1371/journal.pgen.1009689.g006

MOA-seq footprints occur at chromatin interaction and distal enhancer regions

In addition to gene-proximal promoters, distal enhancers and long range chromatin contacts play important but not fully understood roles in gene regulation. Proximity ligation methods for detecting 3D chromatin contacts, such as HiC and ChIA-PET have been instrumental in the identification of long distance chromatin interactions. We analyzed MOA-seq coverage at previously characterized candidate enhancer and chromatin interaction sites, summarized in Fig 7. We found that the majority of maize candidate enhancers in husk and inner shoot tissue [13] and long-range chromatin interaction sites in seedlings [41] overlapped with MOA peaks (Fig 7A). To explore these at a regional level, we examined two well-characterized long-distance regulatory elements in maize around the genes tb1 and ub3, both of which are related to important agronomic traits [42,43]. We analyzed the MF profile at these long-distance chromatin interaction regions, and compared it to that of ATAC-seq. KERNEL ROW NUMBER4 (KRN4) is an intergenic quantitative trait locus that presumably contributes cis-regulatory elements to the control of the inflorescence gene UNBRANCHED3 (UB3). KRN4 is bound by the UB2 (Fig 7B, E1) transcription complexes and interacts with the UB3 promoter (Fig 7B, light blue highlight) by at least three duplex interactions (Fig 7B, A1, A2 and A3, yellow highlight) that affect UB3 expression [44,41,45]. We observed clear MF footprints at all previously annotated KRN4/UB2/UB3 interactions [44,45] in this region, even in those where ATAC-seq coverage was relatively low (Fig 7B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. MOA-seq footprints highlight known long-distance chromatin interaction.

Genome-wide overlap analysis of base pairs shared between MOA and previously defined enhancers (husk, IST, from [13]) and long-range interaction sites (seedling, from [41]). (B-C) Examples of MOA-seq coverage and peaks (Mfc, Mfp, dark orange) at known long-range interaction sites compared to ATAC-seq coverage and peaks (Ac, Ap, gray) at the (B) UB-KRN4 region and the (C) TB1 region. (B) Panel shows chromatin interaction (CI, black lines) loops inferred from HiC, UB2 ChIP-seq, and H3K27ac- and H3K4me3-mediated ChIA-PET-seq data for the KRN4 and UB3 region, three high-confidence remote sites (A1, A2 and A3; yellow columns), and the locations of KRN43.1 (light green column), and UB3 promoter (light blue column) [45]. (C) Interaction region of the -65 Kb distal enhancer region (left zoom-in panel) and promoter region (right zoom-in panel) of TB1. Motif names, positions (Mo) and consensus logos, identified within MOA-seq footprints, are indicated. Dashed lines connect consensus motifs identified in both interacting regions.The best database hit (footprintdb) with the similarity p-value are indicated.

https://doi.org/10.1371/journal.pgen.1009689.g007

The 65 Kb upstream distal enhancer region of the TB1 gene includes large-effect quantitative trait loci associated with plant morphology traits [46]. The distal control region includes a hopscotch transposon insertion site which enhances the expression of TB1 in domesticated maize compared to its wild ancestor, teosinte, affecting apical dominance and tiller bud dormancy [42]. Inspecting the MOA-seq profiles, we detected strong MFs in the proximal tb1 promoter and distal upstream region (Fig 7C, light blue highlights). Given their long-distance interaction and the known regulatory role of the hopscotch insertion site in tb1 expression, we speculated that both regions may share regulatory elements. To test this hypothesis, we inspected motifs identified within MFs for the 600 bp regions directly upstream of the hopscotch insertion (Fig 7C, distal enhancer zoomed area) and tb1 TS site (Fig 7C, proximal promoter zoomed area). Among the motifs found at these two relatively small regions, we found four different motifs, two dyad and two oligo motifs, that were indeed shared between the distal and proximal regions. Comparing these motif sequences to the footprintdb [37] database revealed significant overlaps with consensus motifs of dym45 with the C2H2-type zinc finger, dym46 with a NAC, om011 with an ABI, and om016 with the MYB TF family. Taken together, all of the findings reported here demonstrate the potential for MOA-seq to pinpoint candidate promoter TF-binding sites at promoters and intergenic loci.

Plant materials

Earshoots from B73 wild-type maize were harvested from field-grown plants during mid-morning. The tissue harvesting for materials used in this paper is the same as that used for nuclease sensitivity profiling, DNS-seq, as previously described [7]. For the materials used in this study, field-grown (Mission Road Research Facility, Dept. Biological Science, Tallahassee, FL, USA) earshoots were harvested between 9-11am on sunny days of June 8th, 9th, and 13th of 2015. Earshoot samples were rapidly harvested, measured, and immediately frozen in liquid nitrogen, pooled by date of harvest and size class and stored at -80°C. Multiple earshoots were ground frozen in liquid nitrogen, followed by subsequent aliquoting of the frozen powder for multiple preparation replicates. Size A = pool of 15–20 earshoots of length 0.5–1 cm, size B = pool of 15–20 earshoots of length 1–2 cm, size C = pool of 8 earshoots of length 2–3 cm, and size D = pool of 8 earshoots of length 3–5 cm.

MOA-seq and MNase control library sequencing

The MOA-seq bench protocol is provided (S1 File). It includes tissue fixation, nuclei isolation, MNAse digestion, library preparation, and library size selection. The size-selected indexed libraries were subjected to an equimolar pool of 10 libraries (summarized in S1 Table). The 10 libraries correspond to 2 replicates of each size class, A, B, C, and D, and technical replicates of the two B samples. The technical replicates are from the production of two different libraries made from the light digest pools for B biorep1 and B biorep2. All 10 libraries were separately size-selected, quantified, and subjected to equi-molar pooling for NSG using Illumina HiSEq 2500 paired-end 50-cycle sequencing (FSU College of Medicine, Translational Science Laboratory, https://med.fsu.edu/translationalLab/home). For the MNase control library, total maize DNA was purified using Qiagen plant (greenhouse-grown, v5 leaf blade, B73) DNA purification kit and used for partial digest titrations similar to that for MOA-seq. Specifically, we carried out a partial digest, but with reducing the MNase concentration by 1/64th to achieve the same level of partial digest used for MOA. DNA purification and library construction and size selection was done exactly as described above for MOA-seq and subjected to 100 bp PE sequencing on the MGI-2000(DNBSEQ-G400) platform. Raw sequences from the MOA-seq libraries and MNase control were deposited to the NCBI Sequence Reads Archive (https://www.ncbi.nlm.nih.gov/sra) under the BioProject ID PRJNA477338. The bio-replicate analysis showed high similarity between all (Pearson’s r > 0.9) and earshoot library bio-replicates were pooled into "rep1" or "rep2" for the two bio-replicates, or as "combined" for an earshoot pool of all 10 libraries used for peak calling.

MOA-seq data processing, genome alignment and analysis

The demultiplexed, “raw” sequence data (fastq format) was processed to trim the 3’ adaptor sequences with CutAdapt without filtering for low-quality bases at the 3’ or 5’ end (-q = 0). The minimum overlap length required between read and adapter was set to 1 (-O = 1). All other parameters were set to their defaults in CutAdapt version 1.16 [52]. Trimmed reads were aligned to the B73v3 and B73v5 reference genome assemblies with bowtie 2 using the following options: end-to-end,—no-mixed,—no-discordant, minimum fragment length = 0 (-I = 0), maximum fragment length = 1000 (-X = 1000). All other parameters were set to their defaults in Bowtie 2 version 2.3.1 [53]. Aligned reads were processed using various programs from the BEDTools suite [54,55], as described below. The bam files were combined with "samtools merge" across samples into pools for all biological and technical replicates for all the "B" sized earshoots (indicated as B_all or Bc_all) or for the pools of rep1, rep2, or all (combined). Mate-pairs of reads were specified using “samtools fixmate” and suspected PCR duplicates were marked with “samtools markdup.” Combined, deduplicated bam files were then converted to bed format with “bedtools bamtobed”. Further analyses were done on both resulting aligned fragments of all sizes, and ≤80 bp fragments. Fragment densities were calculated genomewide using “bedtools genomecov.” Coverage files in 10 bp windows were calculated utilizing “bedtools makewindows” and “bedtools intersect.” Resulting bedgraph files were then converted to bigwig format with “bedGraphToBigWig” (S1 Zip File, pooled coverage file used for peak segmentation).

To further increase the resolution of MNase-sensitive footprints, we calculated the fragment centers (also called "frenters" or "midpoints") by extracting the geometric centers of each aligned fragment with “awk,” and intersected these midpoints with 21 bp genomic windows with a step size of 5 bp. Resulting bedgraph files were then converted to bigwig format for UCSC genome browser display using “bedGraphToBigWig”.

iSeg peak calling and analysis

We applied the iSeg peak calling algorithm [56] to detect peaks in each MOA-seq read coverage file using a range of "BC" stringencies as previously applied to DNS-seq data [7]. The parameters and their values for segmenting MOA-seq data are listed below: (1) iSeg Biological Cutoff (BC) value: A BC value of 1.0 means that the height of a significant segment has to be greater than 1.0 x standard deviation of the data mean. We set the BC values to: 1.0, 2.0, 3.0, 4.0, 5.0, 7.0 and 9.0 for segmenting MOA-seq. (2) Minimum and maximum window length: These values are set allowing iSeg to scan a large number of segments within these window lengths. The settings used are minimum window length of 20 bp and maximum window length of 120 bp. (3) Standard Deviation (SD) and median absolute deviation (MAD): While not necessary in general, these values are pre-computed in our application by calculating sample SD and sample MAD for each MOA-seq file when calling its peaks. The global non-zero regions are determined by the regions in the combined COVERAGE pool that had at least one read in any sample. When computing the sample statistics, we removed the global zero regions in each MOA-seq to reduce the degree of distortion caused by sparsity. The INPUT files were bedGraph format with four tab-separated entries: chrom, chromStart, chromEnd, dataValue. The iSeg OUTPUT files are text files, tab-delimited tables with six columns (chrom, SegmentStart, SegmentEnd, meanHeight, t-statistics, p-value). If needed, adjacent book-ended peaks or those separated by 1 bp were merged to produce the final peaks BED files. We used the first three columns to generate BED format files (S1 Zip File for MOA-seq coverage peaks; S2 Zip File for MFs) and converted it into bigBed format for browser display purposes. Peak overlap analyses were done using the "Table Browser" tools from the UCSC Genome Browser [57] using the "Intersection" function with the "any overlap" setting. To optimize the comparisons across different datasets and genome assemblies (B73v3 and B73v5) we used a genomic fraction equivalency approach to select peaks that captured 0.5–0.7% of the genome for MOA coverage or control profiles and 0.1–0.2% of the genome for MOA fragment centers (MFs) or control fragment center profiles (S2 Table).

RNA purification, sequencing, and analysis

RNA was prepared from ~0.1g/samples of frozen ground powder from the same aliquots described above for MOA-seq, with replicates matching the MOA seq sample design. The mRNA was purified using RNeasy Plant MiniKit (Qiagen 74904) and submitted for library preparation and sequencing (Molecular Cloning Facility, Florida State University, Tallahassee, FL, USA). The 10 RNA-seq libraries were prepared and sequenced and available from NCBI Sequence Reads Archive BioProject PRJNA477338. For analysis of RNA-seq data, trimmed and QC filtered Fastq sequence reads were mapped to the B73v3 reference genome (AGPv3.22) using STAR v. 2.7.1a [58] in two-pass mode with additional parameters:—outSAMstrandField intronMotif,—outFilterType BySJout,—outFilterIntronMotifs RemoveNoncanonical,—quantMode TranscriptomeSAM GeneCounts. Unique reads were filtered by mapping quality (q13) and PCR duplicates removed using Samtools (v. 1.3.1). Transcript accumulation was analyzed in R (v. 3.6) using the DEseq2 software (v. 1.24.0) [59]. The combined pooled sample was used to examine the transcript abundance versus MOA coverage (Fig 4D).

Comparative analysis of MOA-seq to other genomic annotations

Several published or shared datasets were analyzed. We obtained a recently published dataset of ATAC-seq peaks from nuclei isolated from 1 cm field-grown earshoots [31,32]. For knotted1 [34] and fasciated ear4 [33] we obtained published ChIP-seq peaks and used their genomic coordinates as central features to plot the average local MOA-seq coverage. For conserved noncoding sequences, CNS, we obtained a public but unpublished dataset (courtesy of Liang, Zhikai; Schnable, James (2019): Conserved Noncoding Sequences. figshare. Dataset: https://doi.org/10.6084/m9.figshare.7804268.v1).

This CNS dataset, called CNS_ZL_v319c in this study, was produced as part of a comparison with the STAG-CNS methods previously described by Lai et al. [60].

Motif discovery, nomenclature, and analysis

Motif discovery was done using the RSAT ChIP motif discovery tool and integrated [61] [62]. The MOA-seq footprints (MFs) from iSeg peaks called at bc7 (S4 File, peaks BED file) were used to specify genomic sequences as input for RSAT. For the RSAT pipeline, peak regions were used to extract genomic DNA sequences using bedtools (vers. 2.27.1). Some of the peaks were below the minimum size limit for RSAT input (24 bp). For these, we expanded the peaks to 24 bp to retain them in the input data. The data was analyzed using default settings except: the oligo length for all analyzes was set to 6 and 7, and the markov order was set to M = 1. For the motif discovery the following TF motif databases were included: ArabidopsisPBM (11/2015), Athamap (11/2015), cisBP (Ath, 6/2015 v1.02), Cistrome by DAP-seq (6/2016), footprintDB-plants (6/2018) and JASPAR core nonredundant plant motifs (2018). NsitePL with the PlantProm database [38], which includes 3,032 previously identified plant TF-binding sites found in 576 experimentally tested promoter sequences, was further used to identify putative TFs and motifs underlying MOA-seq peaks. The full results RSAT reports for the two methods are in the supporting information zip files (S3 Zip File for dyads; S4 Zip File for oligos) and combined into a browser-compatible all motifs BED file (S4 File).