Transcription factor binding site divergence across maize inbred lines drives transcriptional and phenotypic variation

Abstract

Regulatory elements are important constituents of plant genomes that have shaped ancient and modern crops. Their identification, function, and diversity in crop genomes however are poorly characterized, thus limiting our ability to harness their power for further agricultural advances using induced or natural variation. Here, we use DNA affinity purification-sequencing (DAP-seq) to map transcription factor (TF) binding events for 200 maize TFs belonging to 30 distinct families and heterodimer pairs in two distinct inbred lines historically used for maize hybrid plant production, providing empirical binding site annotation for 5.3% of the maize genome. TF binding site comparison in B73 and Mo17 inbreds reveals widespread differences, driven largely by structural variation, that correlate with gene expression changes. TF binding site presence-absence variation helps clarify complex QTL such as vgt1, an important determinant of maize flowering time, and DICE, a distal enhancer involved in herbivore resistance. Modification of TF binding regions via CRISPR-Cas9 mediated editing alters target gene expression and phenotype. Our functional catalog of maize TF binding events enables collective and comparative TF binding analysis, and highlights its value for agricultural improvement.

Figure 1.. Large scale DAP-seq profiling of maize TFs provides high quality genome-wide binding site data for B73v5

a. Heatmap showing Pearson correlation of TF binding sites genome wide (10 bp bins). Side annotation colors correspond to TF structural families. b. Phylogenetic tree showing sub-clade specificity among sequence motifs and target sites of maize BZIP proteins. Loci shown include known targets of maize BZIPs or their Arabidopsis homologs. c. Heatmap showing TFs whose binding sites were enriched for overlap with published maize GWAS hits for various traits.

Figure 2.. DoubleDAP-seq analysis of Q-A-R type BHLHs

a. Clustal alignment of amino acids in BHLH DNA binding domain of Q-A-R BHLHs (light green; group VIII) and group VIIb BHLHs (light blue) with QAR and HER residues shown. b. AlphaFold predicted ribbon diagrams of BHLH85 (Q-A-R), BHLH125 (H-E-R), and the heterodimer. Empirically determined DAP-seq motifs are shown below each structure. c. Number of peaks and top motifs called for DAP-seq and doubleDAP-seq datasets. d. Pearson correlation of genome-wide binding events. e. Genome browser screenshots of binding by homo- and heterodimers. f. Summary of protein-protein interactions based on DAP-seq experiments.

Figure 3.. Maize TF diversity panel reveals functional cis-regulatory modules

a. Maize TF diversity panel is represented by mostly distinct motifs. b. Stacked bar graph showing the number of DAP-CRMs (clusters of TF binding sites) and their prevalence in various gene features. c. Genome browser screenshot of TF binding peaks within the GRASSY TILLERS1 locus which contains eight distinct DAP-CRMs, each with three or more TFs. The estimated region of the known prolificacy QTL is also shown. d. Overlap of DAP-CRMs with various orthogonal functional datasets. e. Heatmap showing DAP-CRMs with specific combinations of TFs that are enriched for different GWAS traits at a significance threshold of 1e-7.

Figure 4.. Genotype-specific peaks are prevalent in B73 and Mo17 and explain genetically defined QTL

a. Stacked bar graph showing the percentage of B73-specific peaks (yellow) and shared peaks found in both B73 and Mo17 (gray). Dotted line indicates the average percentage of shared peaks. b. Percentage of B73-specific (light yellow) and shared (pink) peaks that overlap with duplicated regions (DUPs), small indels less than 50bp (INDELs), SNPs, and structural variants (SVs). Each datapoint corresponds to the percentage of peaks from an individual TF overlapping the indicated category. c. JBrowse 2 screenshot of vgt1-RAP2.7 locus showing three MADS69 binding events upstream of RAP2.7 in B73v5, one of which is located in the genetically defined vgt1 enhancer. Lower left panel shows a close-up of the region in Mo17 containing the MITE transposon and the corresponding region in B73v5 showing the MADS69 peak. d. Alignment of the MADS69 CArG-box motif with B73 and Mo17 sequences. The MITE inserts in the middle of the motif, eliminating high information content nucleotides within the motif. e. Box plots showing motif score differences for B73-specific peaks of each of the TFs in the diversity panel. For most TFs, B73 sequences have higher motif scores than the corresponding sequences in Mo17 within B73-specific peaks regions. f. JBrowse 2 genome browser screenshot showing example of 12bp indel that results in ARF TF binding in B73 but not Mo17 at the ARF15 locus. g. Heatmap showing enrichment scores of differentially expressed putative target genes of shared and B73-specific peaks in B73 vs. Mo17 in various tissues (Zhou et al., 2019). B73-specific target genes are indicated with a magenta bar and shared target genes are indicated with a dark blue bar. Heatmap color scale maps to −log 10 transformed FDR-adjusted P value of Fisher’s exact test.

Figure 5.. CRISPR induced cis-regulatory variation drives expression and phenotypic differences.

a. JBrowse 2 genome browser screenshot of ~12kb region surrounding the DICE enhancer in Mo17 and B73v3. Comparative DAP-seq data near the DICE enhancer reveals two conserved DAP-CRMs (pink highlighted area) and one Mo17-specific CRM (purple highlighted area) that appears to be a partial segmental duplication of the upstream CRM and binding sites (CRM119798). b. RNA-seq data from 11-day old seedlings from Zhou et al., 2019 showing expression levels (TPM: transcripts per million) of various BX genes located near the DICE enhancer. Gene order is same as on chromosome. Mo17 shows 51-fold greater levels of BX1 expression respectively relative to B73. c. CRISPR editing of Mo17 sequences using multiplexed guides near the DICE enhancer revealed specific TF binding sites important for BX1 expression. Relative BX1 qRT-PCR expression for six independent alleles is shown on right. Error bars represent standard deviation. **** adjusted pvalue < 0.0001, ** adjusted pvalue < 0.001. d. Schematic depicting individual enhancer components that contribute to enhanced expression of BX1 in Mo17.