Charting epigenomic landscapes during wheat embryogenesis
We applied CUT&Tag (Cleavage Under Targets and Tagmentation) and ATAC-seq (Assay for Transposase Accessible Chromatin with high-throughput sequencing) to map the histone modification and chromatin accessibility during plant embryogenesis [26, 27]. A pilot run validated the effectiveness of the above methods. First, CUT&Tag results showed high reproducibility between experimental replicates (Pearson correlation coefficient R > 0.94) and high correlation with published chromatin immunoprecipitation (ChIP)-seq results (R > 0.8) but also high repeatability (R > 0.94 between two experimental replicates). Second, we validated that the methods were sensitive to lower nuclei input. Robust CUT&Tag results were obtained using 1000 nuclei at a sequencing depth of ~ 1/3–1/20 of ChIP-seq experiments (Additional file 1: Fig. S1a-c). ATAC-seq from about 5000 nuclei showed good signal-to-noise ratio, reasonable fragment size and genomic distribution, and fully met the demand of transcription factor (TF) footprint calling (Additional file 1: Fig. S1d-h).
Having the techniques established, we generated a ‘reference epigenome’ of wheat embryogenesis across eight characteristic stages (DPA 0, 2, 4, 6, 8, 12, 16, 22, days post anthesis). Embryo sacs were used for early embryogenesis (DPA0-4), while dissected embryos were sampled for later stages (DPA6-22). The epigenome includes seven histone modifications (H3K4me1, H3K4me3, H3K9ac, H3K9me2, H3K27ac, H3K27me3, and H3K36me3), occupancy of histone variant H2A.Z and RNA polymerase II, and chromatin accessibility, as well as transcriptomes (Fig. 1a). Epigenomic and transcriptome data contained two or three biological replicates, respectively. Globally, chromatin marks exhibited highly characteristic distribution patterns correlated with gene transcription (Additional file 1: Fig. S2a, b).
Plant embryogenesis generally divides into three stages: early embryogenesis, during which a totipotent zygote is formed following fertilization; mid embryogenesis, following which the major cell lineages and body pattern are built; and late embryogenesis, during which nutrition is accumulated in the mature embryo followed by seed dormancy [1, 2]. Unsupervised principle component analysis (PCA) and hierarchical clustering of RNA-seq and ATAC-seq revealed a continuous trajectory (Fig. 1b, e), whereas histone modifications showed a punctured stage-specific transition pattern (Additional file 1: Fig. S2c). GO enrichment of temporally correlated gene modules revealed that DNA replication was overrepresented during early embryogenesis, organ/tissue specification genes were enriched during mid-embryogenesis, and nutrient accumulation and dormancy-related genes were expressed more frequently during late embryogenesis (Fig. 1c). Known factors involved in embryo formation, such as WOX11, LEC2, ARR6, and ABI5, but not other genes, such as YUC1, showed conserved embryonic expression patterns between wheat, Brachypodium, and Arabidopsis (Fig. 1d). Such patterns indicate the evolutionary conservation and divergence of monocots and dicots embryogenesis [20, 21, 28]. Accessible chromatin was progressively established until DPA8, the transition stage for fate patterning, followed by gradually decreasing afterward (Fig. 1f). Notably, the majority of ATAC-seq peaks (an indicator of accessible chromatin) were located more than 3 Kb (20 Kb in specific cases) away from the transcriptional start site (TSS) of wheat genes (Fig. 1g and Additional file 1: Fig. S1f); the distance was significantly distant than plants with a small genome [26, 29]. Interestingly, the ATAC-seq peak-TSS distance and proportion of distal ATAC-seq peaks increased as a function of genome size in plant species. For example, the median distance was 129 bp in Arabidopsis with a genome of 125 Mb, whereas the value increased to 4620 bp in barley with a ~ 5 Gb genome (Additional file 1: Fig. S2d) [29]. Consistently, the largest B subgenome of hexaploid wheat has more distal ATAC-seq peaks than the A and D subgenomes (Fig. 1g) [25].
Taken together, we have generated a comprehensive set of epigenomic and transcriptomic datasets, representing a “reference epigenome” of wheat embryogenesis. The concordant changes of histone modification and chromatin accessibility with gene transcription facilitate elucidating the epigenomic basis of diverse regulatory events in wheat embryogenesis.
Distinct proximal and distal chromatin accessibility dynamics
Accessible chromatin exposes regulatory DNA sequences to transcription factors, which provides an effective way to elucidate transcriptional regulation [30]. In total, we identified 1,315,547 accessible chromatin regions (ACRs) from ATAC-seq, which were further categorized into gene body (g), promoter (p), and distal (d) ACRs based on the location relative to genes (Fig. 2a). As expected, both pACRs and gACRs were significantly associated with transcriptional activation (Fig. 2a, 2b, and Additional file 2: Dataset S1). Interestingly, the gain and loss of pACRs and dACRs varied during embryogenesis (Fig. 2c and Additional file 1: Fig. S3a). Many genes gained pACRs after fertilization at DPA2 but dramatically lost them during late embryogenesis at DPA12 and DPA16 (Fig. 2c). Significantly, dACRs, accounting for about 75% of total ACRs, exhibited a sharp “burst” at DPA8 and then quickly declined at DPA12 (Fig. 2c). This transient increase in chromatin accessibility is specific to the dACRs as both peak number and the signal intensity of ACRs at the proximal region were comparable among DPA6, DPA8, and DPA12 (Additional file 1: Fig. S3b, S3c).
We further analyzed the feature of ‘transient’ dACRs at DPA8 as compared to “constant” dACRs that presented at DPA8 and either of DPA6 or DPA12 (Fig. 2d, e and Additional file 1: Fig. S3d, S3e). The constant dACR region is generally enriched for active chromatin states, such as Pol II and H3K27ac marked regions, while depleted with TEs and heterochromatin-related histone modifications, such as H3K9me2 (Additional file 1: Fig. S3d, S3e). By contrast, transient dACRs harbored a higher proportion of TEs (~ 75%) compared with constant dACRs (~ 40%) and were relatively enriched for H3K9me2 (Fig. 2e, Additional file 1: Fig. S3d). Consistent with gained dACRs, more non-coding RNA (ncRNA) were detected at DPA8 as compared to DPA6 and DPA12 (Fig. 2f, Additional file 1: Fig. S3f). Of note, more than 74% of transient dACRs marked ncRNA (n = 6,412) are transcribed from TEs loci (Fig. 2g). Thus, the burst of certain dACR at DPA8 may be correlated with the transient expression of TEs during embryogenesis. Besides the temporal dynamics, subgenomes of hexaploid wheat exhibited considerable differences in chromatin accessibility, especially at the centromere and telomere, in addition to the low-collinear chromosome regions (Fig. 2h and Additional file 1: Fig. S3g).
In summary, proximal and distal chromatin accessibility undergoes distinct reprogramming during embryogenesis. Interestingly, a large proportion of distal regulation and asymmetry of chromatin accessibility among subgenomes exhibit uniqueness for the large hexaploid genome.
Species-specific chromatin remodeling during embryogenesis
To infer chromatin states, we integrated profiles of eight histone modifications and Pol II occupancy from four embryonic developmental stages (DPA 0, 2, 4, 8) using ChromHMM [31, 32]. Totally, 12 chromatin states were categorized, which were further grouped into five major functional classes: promoter (Pr), enhancer-like (EnL), transcriptional (Tr), polycomb group (PcG), and heterochromatin (Hc) (Fig. 3a). Promoter and enhancer-like classes were enriched for H3K27ac, higher chromatin accessibility, histone variant H2A.Z, and Pol II occupancy (Fig. 3a, 3b, Additional file 1: Fig. S4a). As expected, chromatin states were dramatically altered across embryo stages (Fig. 3c), around 10% of the genome changed chromatin state between adjacent developmental stages (Fig. 3d), with the promoter and enhancer-like class varying the most (Fig. 3e).
In mice and humans, H3K4me3 and H3K27me3 are extensively reprogrammed immediately after fertilization, which is pivotal for early embryogenesis [4, 9]. Thus, we investigated the dynamics of the above histone modifications during wheat embryonic development, separating proximal and distal regions as those histone modifications exhibit large distributions in the intergenic regions (Additional file 1: Fig. S2a). Genome-wide H3K4me3 exhibited a moderate change in both proximal and distal regions, as well as H2A.Z (Fig. 3f and Additional file 1: Fig. S4b, S4c). Instead, H3K27ac exhibited a sharp decline at the pro-embryo stage, then restored and maintained high until mid-embryo, followed by gradually decreasing at late and mature embryo stages (Fig. 3f and Additional file 1: Fig. S4b). For H3K27me3, a sharp erasure occurred at the pro-embryo stage, then gradually restored and kept high at late and maturation stages (Fig. 3f and Additional file 1: Fig. S4b). Interestingly, subgenomes of hexaploid wheat exhibited a different abundance of histone modifications even in collinear regions (Fig. 3g and Additional file 1: Fig. S4d, S4e). H3K27me3 exhibited higher variation among sub-genomes than active histone marks such as H3K27ac, H3K4me3, and histone variant H2A.Z (Fig. 3g and Additional file 1: Fig. S4d, S4e).
Together, we revealed species-specific dynamic patterns of histone modification reprogramming during early embryogenesis in wheat (Fig. 3f, h, and Fig S4b, S4c). Different subgenomes of polyploidy wheat have varied histone modification, especially for repressive histone modification H3K27me3 (Fig. 3g, and Fig S4d, S4e).
Reprogramming of H3K27ac, H3K27me3, and chromatin accessibility shapes gene expression dynamics in early embryogenesis
Upon fertilization, the maternal program is silenced, followed by zygotic gene activation, during which epigenetic reconfiguration occurs in both animals and higher plants [9, 23, 34]. During wheat embryogenesis, gene activation was mainly observed at DPA2 and 4, including cell cycle and cytokine signaling genes associated with zygotic activation (Fig. 4a, Additional file 1: Fig. S5a and Additional file 3: Dataset S2). Furthermore, active histone modification H3K27ac and H3K4me3 suffered an overall lopsided decrease at down-regulated genes, while nearly no change occurred at up-regulated genes (Fig. 4b and Additional file 1: Fig. S5a), which is different from the counterparts in mammals [4, 9]. H3K27me3 decreased at both up-and down-regulated genes, indicating that H3K27me3 may not be the major regulator for the divergent expression of genes. By contrast, chromatin accessibility did not change much at downregulated genes but gained at up-regulated genes (Fig. 4b). Such profiles indicate that loss of active histone modification may contribute to gene silencing, while gain of chromatin accessibility with reduction of H3K27me3 may trigger gene activation. Indeed, attenuation of either or both H3K27ac and H3K4me3 was extensively enriched for down-regulated genes (Fig. 4c, d, Additional file 1: Fig. S5b and Additional file 4: Dataset S3), with remarkable down-regulation for loss of both modifications. H3K27ac reprogramming was more effective compared with H3K4me3 (Fig. 4d). Genes subjected to both decreases in active histone marks and transcriptional levels were enriched for maternally-silenced genes, such as orthologs of floral homeotic gene AP2, signaling-related gene MAPKKK17, and stigma expressed gene IAA1 (Fig. 4e, Additional file 1: Fig. S5a, b). Whereas genes gained chromatin accessibility exhibited a significant overlap and correlation with transcriptional up-regulation (Fig. 4f, g, Additional file 1: Fig. S5b and Additional file 5: Dataset S4). These genes were enriched for DNA replication, cell cycle, heterochromatin, etc., which are essential for zygotic initiation (Fig. 4e and Additional file 1: Fig. S5a, b).
After zygotic initiation, we observe a genome-wide drop of H3K27me3 at DPA4 but sharply re-built at DPA8 (Fig. 4h, Additional file 1: Fig. S5c and Additional file 6: Dataset S5). For essential embryonic genes such as LEC1, BBM, WOX11, and ARR12, a permissive chromatin environment was observed at pro-embryo (DPA4), with decreased H3K27me3 and increased H3K27ac, and ATAC-seq signal, while a limited change of H3K4me3 from DPA2 (Fig. 4i). About 90% of the H3K27me3 erasure regions at pro-embryo (DPA4) were re-built in the transition stage (DPA8), along with a large amount of de novo gain of H3K27me3 regions (Fig. 4j). The erasure of H3K27me3 was associated with the cytokinin pathway and essential embryonic genes, while the few augment was involved in the auxin-related genes (Fig. 4j and Additional file 1: Fig. S5d, e). In general, H3K27me3 was negatively correlated with transcription (Fig. 4k and Additional file 1: Fig. S2b). However, decreasing H3K27me3 alone could not explain gene activation from DPA2 to DPA4 (Fig. 4k, l). Of note, H3K27me3 reduction was highly overlapped with adding H3K27ac and chromatin accessibility (Fig. 4l). Overlapped genes showed significant activation, with a more dramatic effect by gaining chromatin accessibility (Fig. 4l and Additional file 7: Dataset S6). Furthermore, there were 561 genes whose transcriptional activation exhibited highly synchronous patterns with a gain of chromatin accessibility after DPA4 (Fig. 4m and Additional file 8: Dataset S7). This synchronization likely ensures the proper timing of gene activation. For example, essential gene LEC1 was first activated followed by pluripotency gene ARR12 and later differentiation genes XMXL5 and IPS2 (Fig. 4m and Additional file 8: Dataset S7). Together, this suggests that resetting H3K27me3 lifts the chromatin layer barrier on embryo patterning genes, which is progressively activated along with the gain of chromatin accessibility.
By contrast, the global H3K27me3 was readily found at DAP6 and fully restored at DAP8 (Fig. 4h, i and Fig S5c). It suggests that the re-gain of H3K27me3 might be important for proper embryo development. In Arabidopsis, Polycomb repressive complex 2 (PRC2) is responsible for deposition of H3K27me3; consisting of four subunits including the single-copy gene encoded FERTILIZATION-INDEPENDENT ENDOSPRM (FIE) [35]. In wheat, there are seven genes encoding FIE orthologues (Additional file 1: Fig. S6a), with a group of homoeologs (A/B/D triads) highly expressed during embryogenesis (Additional file 1: Fig. S6a). In situ hybridization further confirmed the high expression of TraesCS7D02G305100 (one of the triads) in embryo at DPA6, while gradually reduced at DPA8 and DPA12 (Additional file 1: Fig. S6b). We further generated mutations of this group of FIE coding triads by genome editing [36, 37]. Various lines with different combinations of mutations of the triads (fie_aabbdd) were genotyped by sequencing (Additional file 1: Fig. S6c). In heterozygous line D47 (fie_aabbDd), all the triads harbored severe mutation with frameshift, while for homozygous lines (fie_aabbdd) B69 and C87, each had a weak mutation for either A or B triads, respectively (Additional file 1: Fig. S6c). Consistently, offspring of D47 line showed severe embryo development defects at the early DPA6 stage (Additional file 1: Fig. S6d), the homozygous seeds of D47 were likely arrested (Additional file 1: Fig. S6e), whereas the homozygous line C89 was fertile but showed delayed embryonic development (Additional file 1: Fig. S6d) and reduced germination rate (Additional file 1: Fig. S6f). In wildtype KN199, the abaxial-adaxial axis of the embryo becomes obvious and the scutellum arises from the abaxial side at DPA8; while the scutellum was still not visible at DPA12 in line C87 (Fig S6d). Thus, PRC2-mediated H3K27me3 deposition is important for wheat embryo development.
Collectively, reprogramming of histone modification, particularly H3K27ac and H3K27me3, and establishing chromatin accessibility at early embryo development correlates with the proper transition from maternal to zygotic process and re-gain of PRC2 deposition of H3K27me3 is vital for the embryo development.
Transcription factors mediated regulatory network orchestrates embryonic patterning
During mid-embryogenesis, the chromatin environment was relatively accessible, marked by high ATAC-seq and H3K27ac signals (Fig. 1f, 3f, h, and Additional file 1: Fig. S3a, S4b). Embryo patterning is gradually formed along with a precisely programmed transcription regulation pyramid [38]. The complex relationship between TFs-target genes and TFs per se could be well exemplified by the continuous transcription trajectories (Fig. 5a, b and Additional file 9: Dataset S8). We calculated the pseudotime of stages covering early- to mature-embryo based on the PCA distance (Fig. 5b). As expected, an excellent association between gene expression and pACR was observed (Fig. 5c), indicating that individual pACRs and TF binding govern the expression of these genes. A comparison of genes expressed in the early time point (clusters 1 and 2) and the latter ones (cluster 3 and 4) revealed a functional diversity, in which the early elevated genes were involved in cell division, whereas the latter ones in cell specification (Fig. 5d). Moreover, motifs of TFs binding sites, whose orthologs can facilitate inducement of cell totipotency in Arabidopsis, were enriched in early built ACRs, such as the binding motifs of LEC1, MYB118, WUS, etc. (Fig. 5e). By contrast, the TFs functioning in seed dormancy mainly occupied the later elevated ACRs, such as ABI5 (Fig. 5e). This result demonstrates the evolutionary conservation of several key factors in embryo patterning in monocot and eudicot.
We further dissect the architecture of the gene regulation network. In addition to binding motif specificity, we filtered target genes to TFs by co-expression pattern (Method; Fig. 5f). As a result, 1158 paired TF-cis motif interactions were identified, including 154 target genes, 695 ACRs, and 191 TFs (Fig. 5f). To better understand how TFs are coordinated during embryo pattern formation, we extracted the TF-TF networks (Fig. 5g). We identified four pivotal classes of TFs at the center of the network, TCPs, ARFs, MYBs, and WOXs. Some TFs, such as WOXs and ARFs, were reported to regulate embryo patterning in other species. Several known genes regulating organ and tissue specification in Arabidopsis were governed by these TF modules, such as IAA14, LAX2, PLT1, and SMXL5. We selected a regulatory module containing LEC1- Myb118- ZHD5- LEC2- BBM for validation. We found LEC1, Myb118, and ZHD15 genes were highly expressed in pro-embryo stages, while LEC2 was activated later, with the peak expression occurring at the transition stage, and BBM gradually induced at mid-embryo, indicating a sequential expressed pattern (Fig. 5h). The luciferase activity from the reporter assay further validated the regulatory circuit along LEC1, MYB118, ZHD5, LEC2, and BBM module (Fig. 5i). Interestingly, TaLEC2 could activate TaBBM expression via the B3 domain binding motifs within pACR of TaBBM (Fig. 5i), which is not reported in model plant Arabidopsis.
Thus, a cohesive and sequential wave of TF regulation governs embryo pattern formation during mid-embryogenesis in wheat.
Chromatin condensation associated with totipotency reduction and organogenesis restriction during embryonic maturation
Along with embryo maturation, cells differentiate and gradually lose totipotency (Additional file 1: Fig. S7a), but unlike animals, extensive organogenesis is prohibited in plants [1,2,3,4]. We asked how the totipotency capability is attenuated and why the terminal organ is not formed during embryo maturation. First, we identified organ identity genes by their specific expression in terminal organs such as roots and leaves. For comparison, totipotent genes were characterized by high expression in mid-embryo or embryonic callus (Fig. 6a and Additional file 10: Dataset S9). GO enrichment of different groups of genes verified the accuracy of gene identity (Additional file 1: Fig. S7b). Next, we evaluated the chromatin landscape of each gene cluster in different organs. Depletion of H3K27me3 on organ identity genes was observed at the respective organ, whereas obvious ATAC-seq signals were detected at both totipotency and organ identity genes in corresponding organs (Fig. 6b, c, Additional file 1: Fig. S7c, d). Surprisingly, low/no H3K27me3 was enriched in totipotent genes even in differentiated organs (Fig. 6c). In embryonic callus, from which new plant will regenerate, the promoter region of both totipotent genes and organ identity genes were opened with high chromatin accessibility and low H3K27me3 level (Fig. 6b, c). Thus, chromatin accessibility and H3K27me3 coordinately regulate the organ identity gene’s activation potential, while chromatin accessibility but not H3K27me3 mediates repression is indispensable for totipotency gene silencing. Indeed, several regeneration-related factors, such as REV [39], SAUR41 [40], CKX7 [41], and DDM1 [42], were gradually diminished during embryo maturation along with closing chromatin status but activated during the callus induction process (Fig. 6d). Last, we investigated the gene regulatory logic underlying plant organ formation. WRKYs and AtHB TFs footprint were enriched in regulatory regions of root and leaf identity genes, respectively, whereas bHLH and GATA TFs footprint were enriched in both gene sets (Additional file 1: Fig. S7e). Remarkably, nearly half of the identified TFs did not show organ-specific expression but a wide-spectrum pattern, indicating that TF expression per se cannot determine the specificity of organ identity (Additional file 1: Fig. S7f and Additional file 11: Dataset S10). For example, although WRKY75 was expressed during embryo maturation and it owns the capability to bind and active AMT1;1, an ammonium transporter coding gene in the root (Fig. 6e, Additional file 1: Fig. S7f), AMT1;1 was silenced due to the H3K27me3 repression and compaction chromatin environment at the promoter (Fig. 6f).
Taken together, H3K27me3 and chromatin compaction mediated epigenetic regulation associated with the proper embryo maturation process, with reduced totipotency capability but not extensive organogenesis in wheat.
Epigenomic regulation contributes to subgenome divergence of polyploid genome and stage-specific transcriptional divergence
As wheat is an allohexaploid plant, we further examined the epigenomic contributions to the evolution and polyploidization. We clustered genes based on evolutionary age (Additional file 1: Fig. S8a) and found that most of the genus Triticum specific genes, especially genes from A. tauschii (DD), did not participate in embryonic development (Additional file 1: Fig. S8a, S8b). For the homeolog triads present in all three subgenomes, chromatin accessibility was associated with the bias expression pattern (Additional file 1: Fig. S9) in addition to histone modifications as reported [43]. As for the bias-expressed homoeologs, suppressed triads were generally more than dominant triads, and D-subgenome suppressed triads were relatively less than A or B-subgenome suppressed. More balanced triads were expressed in the middle period than in early and late embryogenesis (Additional file 1: Fig. S8c). Such observation was generally consistent with the previous report but different in detail [20].
Polyploidization may drive neo- or sub-functionalization of genes to confer phenotypic plasticity, which can be represented by the varied expression pattern. Instead of identifying genes of differential expression at a specific embryonic stage by comparing hexaploid wheat with tetraploid and diploid [20], we searched for genes exhibiting varied expression profiles during embryogenesis. Indexed by Pearson correlation degree between allohexaploid (AABBDD) and ancestors (AA/BB/DD/AABB), genes were categorized into dysfunction, middle, and conserved clusters (Fig. 7a, Additional file 1: Fig. S10a and Additional file 12: Dataset S11). As expected, conserved genes were mostly homoeolog triads, with sub-genome balanced expression, and showed higher sequence conservation and stronger negative selection pressure (Fig. 7b, c and Additional file 1: Fig. S10b-f). However, a considerable proportion (8%, 1096 genes) of dysfunctional genes changed expression in hexaploid wheat as compared to ancestors, but with slight sequence variation and subject to strong negative selection (Fig. 7c). Promoters of these genes, especially around 1.5 ~ 3 Kb upstream of the TSS, were preferably inserted by TEs (Fig. 7d). Intriguingly, TE insertion regions tended to be accessible chromatin and enriched for H3K27ac (Fig. 7e). For example, TaUGT91C1 showed a B-subgenome high expression in hexaploid wheat, which is different from its ancestors (Fig. 7f). Accordingly, the promoter harbored TEs with higher chromatin accessibility are found in the B sub-genome compared to their A or D counterparts.
A previous study on embryogenesis revealed a conserved hourglass model in hexaploid wheat and the tetraploid and diploid ancestry, but slightly different for the phylotypic stage [20]. Here, we find a developmental stage-specific pattern for transcriptional divergence. A sharp mid-developmental transition was presented, separating two conserved early and late developmental stages in the comparison between homoeologs in hexaploid per se and between counterparts of hexaploid and its ancestors (Fig. 7a, g, Additional file 1: Fig. S10a, Additional file 12: Dataset S11 and Additional file 13: Dataset S12). We thus hypothesized that different regulatory factors might drive such divergence. Indeed, H3K27ac mainly contributed to the middle stage rather than the early- and late-stage, whereas H3K27me3 and chromatin accessibility behaved oppositely (Fig. 7h). In addition, chromVAR analysis indicated that key TFs function mainly in the middle embryogenesis (Fig. 7h and Additional file 14: Dataset S13). Thus, the stage-specific pattern of transcriptional divergence may be governed by chromatin layer regulation at the early and late stages, while cis-motifs and TFs function at the middle stage (Fig. 7i).