Yeast one-hybrid protein–DNA interaction assays were conducted as previously described (Taylor-Teeples et al., 2015). The transcription factors were transformed into each yeast strain and the β-galactosidase activity was determined as previously described (Pruneda-Paz et al., 2009). Positive interactions were visually identified as incidence of yellow caused by the presence of ortho-nitrophenyl cleavage from colorless ortho-nitrophenyl-β-D-galactoside by β-galactosidase. The DNA bait strains were, similarly, tested for self-activation prior to screening, under selection but in the absence of any prey vector. A total of 34 E. coli strains harboring different A. thaliana transcription factors (Supplementary Table 1) were arrayed in 96-well plates and plasmids were prepared. proCESA4^Col (539 bp), proCESA4^Bay (539 bp), and proCESA4^Sha (540 bp) were cloned and recombined with reporter genes. Promoter sequences and primers used are described in Supplementary Table 2. Nine of overlapping fragments of CESA4^Col were independently cloned according to Pruneda-Paz et al. (2009). The oligonucleotides used to amplify promoter fragments are described in Supplementary Table 2. The screen was replicated in full to confirm the results and each clone was sequenced to re-confirm identity.

To express recombinant NST2 or SND1 protein, coding sequences were cloned and fused to glutathione S-transferase tag in the pDONR211 vector and then transferred into pDEST15 (Invitrogen). E. coli strain BL21-AI (Invitrogen) transformed with pDEST15-GST:NST2 were grown in liquid media to an OD600 of 0.4, treated with 0.2% L-arabinose to induce expression overnight and harvested by centrifugation the following day. Cells were treated with 1 mg/mL lysozyme on ice for 30 min in minimal volume of 1X PBS buffer and lysed by sonication. Cell lysates were clarified by centrifugation and incubated with 100 μL of glutathione sepharose beads (GE Healthcare, Pittsburg, PA, United States) for 30 min at 4°C with rotation. The beads were transferred to a column, washed with 10 volumes of 1X PBS. Protein was eluted in 100 mM Tris-HCl pH8.0, 100 mM NaCl and 3 mg/mL glutathione buffer and purified protein was re-suspended in 50% glycerol and stored at -80°C.

Three overlapping probes were generated for CESA4^Bay-C, CESA4^Bay-D, CESA4^Bay-E, CESA4^Sha-C, CESA4^Sha-D, CESA4^Sha-E promoter fragments using the same oligonucleotides described in Supplementary Table 2. Reactions were carried out in binding buffer (10 mM Tris, pH7.5, 50 mM KCl, 1 mM DTT, 2.5% glycerol, 5 mM MgCl2, 0.1% IGEPAL CA-630, and 0.05 μg/μl calf thymus DNA). Following the addition of 150 ng of protein from the GST purification eluate, reactions were incubated at room temperature for 30 min. Protein-DNA complexes were separated from the free DNA on 1% agarose/1X TAE gels at 4°C. The agarose gels were stained with ethidium bromide and bands visualized under UV light.

To develop NILs, RILs maintaining heterozygosity at the CESA4 locus in the F6 generation were identified (Loudet et al., 2002). One plant per RIL carrying a heterozygous CESA4 locus was identified via genotyping and selfed to obtain the F7 seeds. Within the F7 plants, pairs of lines were identified that were either homozygous for the CESA4^Bay and CESA4^Sha alleles to generate the NILs used in this study. The nearly isogenic lines analyzed in this study were developed from RIL93 and RIL350 segregating for the CESA4 region. Individual HIF plants were genotyped for the Bay-0 or Sha CESA4 allele by PCR of a 543 bp of CESA4 promoter with primers described in Supplementary Table 2. PCR products were subjected to restriction enzyme digestion with BsrDI (New England BioLabs, Ipswich, MA, United States) at 65°C for 2 h. The Sha allele incudes a BsRDI cut site (CGTTAC| NN) resulting in 401 and 142 bp fragments. The 543 bp Bay-0 allele contains no BrsDI restriction sites and remains undigested. Transcript abundance of CESA4 in the near isogenic lines (NILs) was quantified in 5 cm stems. Stem tissue was frozen in liquid nitrogen and pulverized with metal beads in a Retsch (Haan, Germany) Mixer Mill MM400. RNA extraction and cDNA synthesis was conducted as described above. Primers for CESA4 real-time PCR are described in Supplementary Table 2.

To quantify and compare levels of crystalline cellulose, the Updegraff assay as adapted and described by Kumar and Turner was used (Updegraff, 1969; Kumar and Turner, 2015). Briefly, alcohol-insoluble residue (AIR) samples were first prepared from senesced stem tissue; samples were weighed at this stage for later calculations. Then, 3 mL of an acetic/nitric acid solution (8:1:2, acetic acid: nitric acid: water) was added to AIR samples and incubated in a boiling water bath for 30 min. This step removes hemicellulose and lignin while leaving cellulose microfibrils intact. The remaining cellulosic material was then swelled in 67% sulfuric acid in a boiling water bath for 5 min to disorganize the polymers, and the subsequent monomers were finally analyzed with a sulfuric acid/anthrone colorimetric assay. Absorbance was measured at a wavelength of 620 nm with a SpectraMax M5 and corrected to a known glucose standard to calculate the percent cell wall composition of cellulose.

To investigate the histological effects of expression variation on the cell wall, 100 μm cross sections were sliced using a Leica Biosystems VT1000S Vibrating-blade microtome. Sections were incubated on the bench top for 2 min in a 2% w/v phloroglucinol/ethanol solution (Fisher Scientific, Waltham, MA), mounted in a 1:1 concentrated hydrochloric acid: water solution, and immediately imaged using a Nikon Eclipse E200 microscope and a PixeLink scope camera. For each stem, three sections were imaged, and five cells from three different portions of each stem were measured to get an average of 45 measurements per biological sample.

To qualitatively analyze cellulose crystallinity under polarized light, internode segments were first cut on a vibratome into 100 μm sections, fixed in 2% glutaraldehyde and then embedded in Epon/Araldite (Sigma-Aldrich, St. Louis, MO, United States) before slicing into 0.5 μm sections on an Reichert-Jung, Ultracut E microtome (Vienna, Austria). Polarized light microscopy takes advantage of variation in the passage of light through crystalline structures to uncover discrete differences in crystallinity configurations. Specifically, the LC-PolScope (CRI, Cambridge, MA, United States) employed for this assay uses polarized light to measure birefringent retardance and the intensity of the images generated directly correlates to the retardance value, thus qualitative differences in crystallinity can be observed. Several stem samples from each HIF were imaged and analyzed.

In order to determine a consensus VND, NST, and SND (VNS) protein binding motif, the MEME.txt motif files from O’Malley et al. (2016) were visually aligned based on nucleotide similarity and trimmed to a length of 13 bases. An average for each nucleotide, at each motif position was calculated using a bash script that referenced the aligned and trimmed motif files. The consensus matrix was then loaded into R and the Bioconductor package SeqLogo was used to generate the motif logo. A bash shell script counted each of the VNS motif variant across the A. thaliana TAIR 10 sequence assembly and for each DAP-seq set of binding sites. Sum of the counts determined the total number of VNS motif sites and percentages were calculated as the proportion of each motif variant relative to the sum and multiplied by 100. The process was repeated for the DAP-seq binding site data using the reference genome file as a guide to extract DAP-seq peak fasta sequences. ThenarrowPeak DAP-seq peak data files were then searched iteratively through each of the peak binding sites. The percentage shown for DAP-seq data represents the mean percentage across all binding site files.

Correlation coefficients and tests of significance were calculated using Pearson’s correlation tests for all CESA genes pairwise using replicate gene expression data for Bay-0 and Sha parental lines and RILs using the rcorr function from the Hmisc library in R 3.2.5. Two tailed Student’s t-tests were carried out in Excel to assess gene expression, cell thickness, and cellulose content data sets. BoxPlots were generated with the BoxPlotR Web Tool (Spitzer et al., 2014).

CESA1 (At4g32410), CESA2 (At4g39350), CESA3 (At5g05170), CESA4 (At5g44030), CESA5 (At5g09870), CESA6 (At5g64740), CESA7 (At5g17420), CESA8 (At4g18780), CESA9 (At2g21770), CESA10 (At2g25540), NST1 (At2g46770), NST2 (At3g61910), NST3/SND1 (At1g32770), PP2AA3 (At1g13320), TIP41 (At4g34270).

We explored the observation that CESA gene expression varied among different accessions of A. thaliana. West et al. (2007) measured the expression of 22,746 genes in replicated samples of the Bay-0/Sha RIL population. Above ground tissue of 211 short-day grown plants was assayed after 6 weeks of growth using the Affymetrix ATH1 GeneChip microarray. There was a continuous range of values and normal distributions were observed for the ten CESA genes known to play a role in the biosynthesis of cellulose in secondary cell walls (Figure 1A and Supplementary Figure 1). The midparent values (i.e., the values halfway between the two parents) and the median of the RILs for CESA4, CESA7, and CESA8 (Figure 1A) and the seven other CESA genes (Supplementary Figure 1) were very similar.

The CESA genes have been shown to be highly co-regulated within their functional classes. Three of the primary CESAs, CESA1, 3, and 6, had the most similar gene expression to each other in a meta analysis of gene expression and the same is true among the secondary CESAs: CESA4, 7, and 8 (Persson et al., 2005). Candidate genes for the transcriptional regulation of CESA gene expression have been successfully identified using the same type of analysis in various species (Brown et al., 2005; Persson et al., 2005; Yamaguchi et al., 2010; Ruprecht et al., 2011; Handakumbura et al., 2018). Among the Bay-0/Sha RILs, the expression of CESA4, 7, and 8 were significantly correlated as were the primary wall CESAs (Table 1).

Among the genes whose expression was measured using the ATH1 microarray, 69% were associated with an eQTL and each transcript was mapped to an average of 2.34 loci (West et al., 2007). To pinpoint the cause of CESA gene expression variation among RILs, we searched for eQTL for these genes (Figure 1B). An eQTL for each of the CESAs was mapped to a position of the genome outside of those genes, i.e., trans-eQTL (Figure 1B and Table 2). A trans-eQTL common to all three secondary wall CESA genes was found near the bottom of chromosome 1. Overlapping trans-eQTLs were identified for CESA4 and 8 on chromosomes 2 and 4. An eQTL unique to CESA4 with the greatest LOD score of 17.95 was found on chromosome 5, coincident with the CESA4 genomic locus. CESA4 transcript abundance in RILs varied significantly depending on the presence of the Bay-0 CESA4 promoter allele (CESA4^Bay) or the Sha CESA4 promoter allele (CESA4^Sha) (Figure 1C). No differences were observed for the expression of other secondary wall CESA genes, CESA7 or CESA8, between RILs with either CESA4^Bay or CESA4^Sha (Figure 1C).

Possible mechanisms for such a cis-eQTL include functional polymorphisms in the promoter of the gene in question, and we hypothesized that this change may disrupt interaction of SND1 or NST2, which we previously identified to interact with the CESA4 promoter using yeast one-hybrid (Taylor-Teeples et al., 2015). As such, we sequenced the promoters of Bay-0, Sha, and Col-0 and identified 16 single nucleotide polymorphisms (SNPs) between CESA4^Sha and CESA4^Col and only two SNPs between CESA4^Bay and CESA4^Col (Figure 2). To specify which of the 16 SNPs were likely responsible for the CESA4 cis-eQTL, we screened nine overlapping fragments (A–I) of CESA4^Bay and CESA4^Sha promoters by yeast one-hybrid with CESA4^Col promoter interacting proteins (Figures 2, 3A). In this assay, transcription factor proteins were fused to the Gal4 activation domain and each protein is tested for an interaction with the CESA4 promoter fragments immediately upstream of the lacZ reporter. A positive interaction can result in the cleavage from colorless ortho-nitrophenyl-β -D-galactoside by β-galactosidase resulting in a yellow color. The well-characterized wall thickening regulators SND1 and NST2 interacted with two CESA4^Bay fragments (C and E) but only a single CESA4^Sha fragment (E). One SNP in this region disrupts the second position of a perfect TERE motif in CESA4^Sha fragment C (Figure 3B). This polymorphic fragment failed to interact with SND1 and NST2 (Figure 3A). This suggests the expression differences between CESA4^Bay and CESA4^Sha could have occurred as a consequence of differential binding of NST2 or SND1 proteins to the TERE motif in fragment C.

To further explore the possibility of this regulatory mechanism, EMSA was performed to confirm the differential protein-DNA with probes corresponding to CESA4 promoter fragments C, D, and E in the presence or absence of extracts of Escherichia coli expressing GST-NST2 and GST-SND1 (Figure 3C). As anticipated, differences in mobility were observed with the TERE motif-containing fragment C of CESA4^Bay but not with the corresponding fragment of CESA4^Sha, confirming our yeast one-hybrid observations. Also consistent with our yeast one-hybrid results, fragment D from both accessions did not produce a DNA species with retarded mobility, but the TERE motif-containing fragment E did in both cases. Bacterial extracts harboring the empty GST vector did not produce any comparable shifted species. Taken together, these data suggest that the CESA4^Sha TERE motif polymorphism may result in diminished binding by SND1 and NST2 proteins.

Resolving the effects of the CESA4^Bay and CESA4^Sha cis-regulatory region is confounded by the entirety of the sequence variation between Bay-0 and Sha, which are maintained in different combinations among the RILs. To isolate the influence of the CESA4 locus, we identified RILs with residual heterozygosity at CESA4 and tested derived NILs. We isolated the NILs from heterozygous inbred family 93 (HIF93) and HIF350, segregating for the CESA4 promoter interval that demonstrated differential binding. Transcript abundance in developing stems containing CESA4^Bay was twofold greater than that of CESA4^Sha in both pairs of NILs (Figure 4). These results further support the differential binding between Bay-0 and Sha by NST2 and SND1, thus changing expression of CESA4.

Considering the critical role of CESA4 in plant structure, we wished to evaluate the phenotypic effect of the cis-eQTL. Bay-0 and Sha CESA4 NILs were indistinguishable at the macro level, but a subtle consequence at the cellular level remained a possibility. To explore the effect of the CESA4 cis-eQTL on cellulose, Bay-0 and Sha were first examined with stem histology to measure cell wall thickness. Stem cross sections were indistinguishable between NILs (Figure 5A). All xylem and interfascicular cells had well-defined edges and equivalent phloroglucinol staining. The CESA4^Bay NIL93 exhibited an average cell wall thickness of 3.57 μm, which was slightly greater than the 3.07 μm thickness of their CESA4^Sha NIL93 counterparts. The CESA4^Bay NIL93 samples were quite variable with a standard deviation of 0.45 μm. Cell wall thickness of CESA4^Bay and CESA4^Sha NIL350 were almost identical, 2.86 and 2.85 μm, respectively (Figure 5B).

While no significant differences in overall lignin staining or cell thickness was observed, we hypothesized that the reduction in CESA4 expression may still compromise the cellulose crystalline structure of the NILs carrying the CESA4^Sha allele. To evaluate this, we employed Updegraff assay to quantify crystalline cellulose and polarized light microscopy as a secondary indicator of crystallinity. Percent composition of crystalline cellulose was slightly higher in CESA4^Bay NIL93 than CESA4^Sha NIL93, 28.2% vs. 26.8%, but lower in CESA4^Bay NIL350 that CESA4^Sha NIL350, 27.0% vs. 29.2%. Percent cellulose crystallinity was marginally greater in the Bay-0 parental accession than Sha, 27.0% vs. 25.3% (Supplementary Figure 2). However, the NILs were not significantly different. Polarized light images of all samples mirrored brightfield images and were indistinguishable between NILs (Figure 6). Birefringent retardance was strong and consistent in xylem and interfascicular cells, revealing no difference in quantities or order of crystalline cellulose. The comparable phenotypes between NILs revealed in both the histological and chemical assays suggest the decreased CESA4 expression observed in CESA4^Sha NILs does not disrupt cellulose abundance or crystallinity.

The TERE and SNBE are compatible sequences, independently identified as binding sites of VND and NST NAC proteins (Pyo et al., 2007; Ohashi-Ito et al., 2010; Zhong et al., 2010). The TERE sequence, CTTNAAAGCNA, is consistent with the following internal sequence of the SNBE: CTTNNNNNNNA. To further resolve this binding site we explored the DNA affinity purification sequencing (DAP-seq) data previously generated for A. thaliana transcription factors and Col-0 genomic DNA (O’Malley et al., 2016). We searched the binding peaks of DAP-seq data for available VND, SND, and NST proteins: VND1, 2, 3, 4, 6, SND2, SND3, and NST1. The data included DNA sequence from libraries of DNA where methylcytosines were removed by PCR and unamplified libraries. There is a striking similarity between the top enriched motif described by the DAP-seq and the SNBE/TERE motif (Supplementary Table 3). The core DAP-seq derived motif similarity across all three groups of NACs is three positions flanking seven nucleotides of any sequence: C(G/T)TNNNNNNNA(A/C)G. We refer to this as the VND/NST/SND element (VNS, Figure 7). Interestingly, the motif appears to be palindromic depending on the variable positions. We searched the sequences of the binding sites and found the prefix CTT and AAG suffix to represent more than half of the total occurrences of the VNS, however, each possible sequence was similar to the frequency in the genome (Table 3). Therefore, there does not appear to be a preference for the variable positions. There was a low frequency of an A or a T in the first and last position of the binding site, respectively. The NAC binding site in the CESA4 promoter is consistent with both the TERE and a VNS, but with a T at the last position.