Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Chromatin accessibility landscapes revealed the subgenome-divergent regulation networks during wheat grain development

Chromatin accessibility landscapes revealed the subgenome-divergent regulation networks during... aBIOTECH https://doi.org/10.1007/s42994-023-00095-8 aBIOTECH RESEARCH ARTICLE Chromatin accessibility landscapes revealed the subgenome- divergent regulation networks during wheat grain development 1 1 1 1 1 Hongcui Pei , Yushan Li , Yanhong Liu , Pan Liu , Jialin Zhang , 1 1& Xueni Ren , Zefu Lu Institute of Crop Sciences, Chinese Academy of Agricultural Sciences, Beijing 100081, China Received: 30 September 2022 / Accepted: 22 January 2023 Abstract Development of wheat (Triticum aestivum L.) grain mainly depends on the processes of starch synthesis and storage protein accumulation, which are critical for grain yield and quality. However, the regulatory network underlying the transcriptional and physiological changes of grain development is still not clear. Here, we combined ATAC-seq and RNA-seq to discover the chromatin accessibility and gene expression dynamics during these processes. We found that the chromatin accessibility changes are tightly asso- ciated with differential transcriptomic expressions, and the proportion of distal ACRs was increased gradually during grain development. Specific transcription factor (TF) binding sites were enriched at different stages and were diversified among the 3 subgenomes. We further predicted the potential interactions between key TFs and genes related with starch and storage protein biosynthesis and found different copies of some key TFs played diversified roles. Overall, our findings have provided numerous resources and illustrated the regulatory network during wheat grain development, which would shed light on the improvement of wheat yields and qualities. Keywords Wheat, Chromatin accessibility, Subgenome-divergence, Regulatory network, Grain development INTRODUCTION development is important to understand the molecular basis of wheat yield and quality and discover novel key Wheat (Triticum aestivum L.) is one of the most regulatory factors. important staple crops for global food security, provid- Starch and protein are the main nutrient components ing more than a fifth of the calories and protein con- of wheat grain. A variety of starch biosynthesis genes sumed by humans (http://faostat.fao.org). Wheat grains have been characterized, including ADP-glucose are mainly composed of pericarp, endosperm and pyrophosphorylases (AGPases), ADP-glucose (ADPG) embryo, of which the endosperm is the main source of transporter, granule-bound starch synthases (GBSSs), white flour. The endosperm is a storage organ that starch synthases (SSs), Starch branching enzymes accumulates a large quantity of starch and proteins (SBEs), debranching enzymes (DBEs), starch/a-glucan during grain development (Gu et al. 2015; Liu et al. phosphorylases (PHOs), disproportionating enzymes 2018; She et al. 2010). Therefore, elucidation of the (DPEs), and protein targeting to starch (PTST), which transcriptional regulatory networks during wheat grain coordinately work in starch biosynthesis during wheat grain development (Huang et al. 2021; Kang et al. 2013; Wang et al. 2019; Yamamori et al. 2000). Starch is Hongcui Pei and Yushan Li have contributed equally to this work. synthesized by plants using the enzyme AGPases that reacts glucose 1-phosphate with ATP to form ADPG, the & Correspondence: luzefu@caas.cn (Z. Lu) The Author(s) 2023 aBIOTECH ADPG were subsequent synthetized to amylose and RESULTS amylopectin with different enzymes. Most of the storage proteins in grains are gluten, including high molecular Chromatin accessibility changes are associated weight glutenin subunits (HMW-GSs), low molecular with differential transcriptomic expressions weight glutenin subunits (LMW-GSs) and gliadins, during grain development which influence the dough property (Liu et al. 2005; Ren et al. 2022). The expressions of those key genes in Wheat grains are mainly composed of endosperm and starch and/or gluten biosynthesis pathways are tightly embryo, we performed ATAC-seq and RNA-seq with endosperm or seeds in 5, 9, 15, and 20 days after pol- related with grain yield and quality. A group of tran- scriptional factors, such as TaRSR1 and TabZIP28, were lination (DAP) to explore the chromatin accessibility found to be involved in regulating the expression of and transcriptomic expression dynamics during wheat genes related to starch synthesis in wheat (Liu et al. grain development, and used leaves as control. A total of 2016; Song et al. 2020), and TaGAMyb, TaFUSCA3 and 103,516, 108,846, 101,620, 105,578 and 104,407 high PBF-D regulate different HMW-GS gene (Glu-1Dy, Glu- confident accessible chromatin regions (ACRs) were 1Bx7 and Glu-1) expression through binding specific identified in leaf and grain in DAP5, DAP9, DAP15 and DNA motifs (Guo et al. 2015; Sun et al. 2017; Zhu et al. DAP20, respectively (Fig. 1A). All the ACRs from the five 2018). Besides, many TFs, such as TaNAC019 control tissues are mostly enriched around the TSS regions glutenin and starch accumulation by targeting both (Supplementary Fig. 1A) and showed high correlation HMW-GS genes and the starch synthase gene TaSSIIa to coefficient ([ 0.80) between two biological replicates synergetic regulation of wheat yield and quality (Gao (Supplementary Fig. 1B–D). Similar with the ACR dis- et al. 2021; Liu et al. 2020). tribution in other large genome species, a significant Wheat grain development can be typically defined faction of ACRs distributed 10 Kb away from genes into three stages: the pregrain-filling phase for 10 days apart from most of ACRs located near the gene (Fig. 1B). after pollination, grain-filling phase for 10–20 days after We then grouped the ACRs as genic ACRs (genic ACRs, pollination and desiccation phase then follows beyond overlapping with genebody), proximal ACRs and distal 30 days after pollination (Wan et al. 2008). The transi- ACRs according the distance between ACR centers and tion of different grain development phase is accompa- their nearest genes. For definition of proximal ACRs and nied by dramatic transcriptional and physiological distal ACRs, we employed different distance cutoffs (2 changes. However, little is known about their regulatory Kb, 4 Kb, 5 Kb, 6 Kb, 8 Kb, 10 Kb), and those within each mechanisms and contribution of each subgenome in cutoff away from genes as proximal ACRs, and those larger than cutoff away from genes as distal ACRs allohexaploid wheat. Accessible chromatin regions are putative cis-regulatory elements (CREs) that are tar- (Supplementary Fig. 2). Interestingly, the proportion of geted by TFs and play vital roles in transcriptional distal ACRs in grain samples were all much higher than regulation. It is reported that the CRE variants are in leaf except in DAP15 with 6 Kb as the cutoff (Sup- associated with many important agronomic traits plementary Fig. 2). When using 2 Kb to distinguish (Adamski et al. 2021; Deplancke et al. 2016; Kouzarides them, more than 50% of ACRs were predominately 2007; Liu et al. 2021; Lu et al. 2017, 2019; Rodgers- enriched in the distal regions of genes in DAP20 Melnick et al. 2016; Tian et al. 2021; Wang et al. 2021b). (Fig. 1C). These results suggested that during grain In this study, we unraveled the regulatory network development, more distal regions were predominantly underlying wheat grain development by combing ATAC- opened and these may be critical for the expression of seq (Assay for Transposase Accessible Chromatin genes regulating starch and storage protein sequencing) and RNA-seq with samples from a series biosynthesis. grain developmental stages. Our studies have not only We further examined the sequence variations around revealed the chromatin accessibility landscapes of ACRs using the reported wheat whole-genome genetic wheat grain development, but also provided valuable variation map (VMap 1.0) (Zhou et al. 2020). To avoid insights into the dynamic regulatory network mediated the influence of biased signals around genes, we only by important transcription factors and the divergences picked up ACR whose centers are[ 2 Kb away from of the three subgenomes during grain development in genes. The results showed that lower sequence varia- the allohexaploid wheat. tion around ACR centers (Fig. 1D), which suggested those ACRs are functionally important. Interestingly, we found that variations were higher in DAP20 than the earlier stages, indicating a lower selection pressure of those genes. The Author(s) 2023 aBIOTECH AC B Genic Proximal Distal 0.8 0.6 0.4 0.2 Leaf DAP5 01 12 34 5 6 78 90 10⁶ 10⁸ DAP9 DAP15 Distance between genes and peak centers (Kbp) DAP20 E F Genic Proximal Distal 0.010 ACR present 1 *** *** ACR absent 0.8 0.009 0.6 *** *** *** 0.008 0.4 0.2 0.007 0.006 0 -2.0Kbp TSS 2.0Kbp Leaf DAP5 DAP9 DAP15 DAP20 GI diffACRs DEGs H Scaled TIS Scaled TPM Cluster organic acid metabolic process oxoacid metabolic process −1.5 0 1.5 −1.5 0 1.5 carboxylic acid metabolic process dACRs in Prmoter DEGs carbohydrate metabolic process Leafs cellular macromolecule biosynthetic process cellular amide metabolic process p < 2.2e-16 peptide metabolic process Leafs amide biosynthetic process DAP5s peptide biosynthetic process gene_expression p < 3.2e-5 glucan biosynthetic process organonitrogen compound biosynthetic process DAP5s carbohydrate biosynthetic process DAP9s oligosaccharide metabolic process disaccharide metabolic process p < 4.0e-4 DAP9s sucrose metabolic process protein modification response to light stimulus p-value DAP5&9s cellular response to radiation 525 3e−04 DAP5&9s cellular response tolight stimulus p < 9.8e-8 cellular response to abiotic stimulus seed development 2e−04 12753 embryo development DAP15s DAP15s protein phosphorylation 1e−04 p < 0.01 phosphorus metabolic process regulation of transcription regulation of RNA metabolic process 13809 regulation of primary metabolic process DAP20s DAP20s regulation of nucleic acid templated transcription GeneRatio regulation of nitrogen compound metabolic process p < 3.5e-10 1211 0.1 phosphorylation DAP15&20s response to water 0.2 27650 response to oxygen containing compound 0.3 DAP15&20s 4285 response to inorganic substance response to acid chemical p < 2.2e-16 13326 response to abiotic stimulus AllGs AllGs p < 0.3 Fig. 1 Chromatin accessibilities are relative to the differential gene expression during wheat grain development. A UpSet Plot showing the ACR distribution across five tissues (DAP, days after pollination). B Density plot showing the distribution between ACR centers and their nearest genes. Black dot line indicates 10 Kb away from genes. C Distribution of ACRs on different genomic features. Genic, ACRs with their centers in gene bodies; Proximal, ACRs with their centers \ = 2 Kb upstream of TSS or \ = 2 Kb downstream of TES; Distal, ACR with their centers [ 2 Kb away from genes. D Distribution of sequence variations around ACR centers. SNPs from Zhou et al. were used. E The expression of genes with ACRs in their 2 Kb nearby regions (ACR present) and without ACRs in their 2 Kb nearby regions (ACR absent) across the five tissues. F Distribution of eight differential ACR (diffACR) clusters on different genomic features. The definition of each group is same with C. G The overlapping information of eight diffACR clusters related genes and eight DEG clusters. H Chromatin accessibilities of regions from 2 Kb upstream to 100 bp downstream of TSS and associated genes’ expression levels. Chromatin accessibilities were represented by Tn5 transposome integration sites (TISs) and expression levels by transcripts per million (TPM). I Gene ontology (GO) enrichment analysis for the overlapped gene groups in G The Author(s) 2023 Leaf DAP5 DAP9 DAP15 DAP20 DAP5 DAP9 DAP15 DAP20 Endosperm Whole seed Leaf DAP5 DAP9 DAP15 DAP20 input Leaf DAP5 DAP9 DAP15 DAP20 Leafs DAP5s Leaf DAP5 DAP9s DAP9 DAP5&9s DAP15 DAP20 DAP15s Leaf DAP20s DAP5 DAP9 DAP15&20s DAP15 AllGs DAP20 Leafs DAP5s DAP9s DAP5&9s DAP15s DAP20s DAP15&20s Relative SNP density No. of Intersections log2(TPM+1) Density 0.0 0.1 0.2 0.3 0.4 0.5 Percentage overlapping genome feature Percentage overlapping genome feature aBIOTECH The chromatin accessibility changes are usually Subgenome-divergent regulation during wheat associated with differential transcriptomic expressions grain development (Chen et al. 2019; Pajoro et al. 2014). To study the relationship between chromatin accessibilities and gene Transcription factors (TFs) play important roles in the regulation of gene expression, while ACRs are predom- expressions, we classified the total wheat genes based on the presences of ACRs within 2 Kb away from genes. inantly TF-binding sites (Liu and Bergmann 2021; Marand et al. 2021). To identify the key TFs responsible The results showed that the expression of genes asso- ciated with ACRs were significantly higher than those for the gene expression dynamics during wheat grain development, we analyzed the enrichments of the 40 without ACRs in all five tissues (Fig. 1E). To identify the putative CREs responsible for the transcriptomic JASPAR cluster motifs (Castro-Mondragon et al. 2022)in expression changes, we identified the differential ACRs the eight diffACRs clusters. The results showed that ARF, (diffACRs) by pair-wise comparisons and then grouped BPC and ERF were enriched in leaf-specific ACRs; ABI, the diffACRs into eight specific clusters, including leaf- ATHB, NAC, MYB, SPL, STZ and bZIP were enriched in specific diffACRs (Leafs), DAP5 specific (DAP5s), DAP9 DAP5s, DAP9s and DAP5&9s; AGL, ARF, E2FA, and REF specific (DAP9s), DAP5 and DAP9 specific (DAP5&9s), were enriched in DAP15s, DAP20s and DAP15&20; Dof was enriched in all grain specific ACRs (Fig. 2A), indi- DAP15 specific (DAP15s), DAP20 specific (DAP20s), DAP15 and DAP20 specifics (DAP15&20s) and all grain cating the specific roles of different TF families in dif- ferent grain developmental stages. specific (AllGs). In total, we identified 31,000, 8346, 9754, 13,829,15,974, 17,253, 39,725 and 1287 specific Bread wheat (Triticum aestivum L.) consists of A, B diffACRs in leafs, DAP5s, DAP9s, DAP5&9s, DAP15s, and D subgenomes derived from three diploid species, DAP20s, DAP15&20s and AllGs groups, respectively and exhibits improved grain yield and nutritional con- (Supplementary Fig. 3A and Supplemental Data Set1). tent relative to diploid and tetraploid wheat (Yang et al. The diffACRs were more enriched in the distal regions 2014). Transcriptomic expression or dynamic modifi- than the total ACRs, indicating accessibilities of distal cation asymmetry of wheat chromosomes contributed to wide adaptability of wheat (Li et al. 2019; Ramirez- regions were more volatile during grain development (Fig. 1F). We further divided DEGs into eight clusters Gonzalez et al. 2018; Wang et al. 2021a). To study the dynamic patterns of gene expression and promoter using the similar strategy and identified the overlaps between DEGs and diffACRs associated genes, which accessibilities of homeologs during wheat grain devel- opment, we conducted the gene expression bias and revealed higher overlapping rates compared to the random shuffled gene sets in all groups except AllGs promoter accessibility bias analysis by ternary plotting (Fig. 2B). 12,218, 11,802, 11,488, 11,319 and 11,597 of (Fig. 1G, Supplementary Fig. 3B, and Supplemental Data Set2). The chromatin accessibilities of promoter regions triad genes showed balanced expression, and more than (– 2 Kb upstream to 100 bp downstream of TSS sites) 75%, 79%, 78%, 76%, 81% and 74% triad genes and the expression of each overlapped gene cluster showed balanced promoter accessibilities in leaf, DAP5, showed similar trends (Fig. 1H), which further con- DAP9, DAP15 and DAP20, respectively (Supplementary firmed that the chromatin accessibility changes are Fig. 4A and B). It should be noted that more unbalanced associated with differential transcriptomic expressions genes were identified in the growing seeds (DAP5, DAP9, DAP15 and DAP20) than in leaves (Supplemental during wheat grain development. To clarify the biological processes mediated by each Data Set3). Among the subgenome differentially expressed genes, the counts of B-suppressed genes and gene group, we performed gene ontology (GO) enrich- ment analysis. The results showed that leaf-specific promoters were more than the other groups, indicating genes were highly enriched in phosphorylation and B subgenome genes were in weaker positions during response to light or abiotic stimulus, while the DAP5, these processes. We further found that only 56%–62% DAP9 and DAP5&9 specific genes were highly enriched genes were balanced in expression among subgenomes in carbohydrate, sucrose metabolic process, in addition, in all tissues and few genes were consistent among amide biosynthesis and nitrogen compound metabolic different stages (Fig. 2C). Similar patterns were also process regulation were enriched in the overlapped found with promoter accessibilities, which indicate genes specific in DAP15 and/or DAP20 (Fig. 1I). These highly diverse regulation among subgenomes during grain development (Fig. 2D). In addition, we conducted enriched GO terms were consistent with the grain development processes. a correlation analysis between gene expression and promoter accessibility, which showed the unbalanced genes were positively correlated with differential The Author(s) 2023 bZIP/RAV bZIP/RAV * * cluster19_NAC cluster19_NAC REF/RAV REF/RAV * * cluster20_NAC cluster20_NAC GT-4 GT-4 * * ARR/GATA ARR/GATA KUA/SRM KUA/SRM cluster31_ARF cluster31_ARF PHL/KAN/HHO PHL/KAN/HHO TGA TGA * * PEND PEND HMG HMG * * cluster12_MYB cluster12_MYB Dof Dof * * cluster39 cluster39 FLC/SOC1/SVP FLC/SOC1/SVP * * cluster16_MYB cluster16_MYB Dof Dof * * KUA/SRM KUA/SRM FLC/SOC1/SVP FLC/SOC1/SVP * * PEND PEND GT-1/AGL GT-1/AGL * * * * PHL/KAN/HHO PHL/KAN/HHO RVE RVE * * * * GT-4 GT-4 WRKY WRKY * * IDD/OBP IDD/OBP TRP TRP * * TCX/SOL TCX/SOL TCX/SOL TCX/SOL * * GT-2 GT-2 GT-2 GT-2 * * cluster12_MYB cluster12_MYB IDD/OBP IDD/OBP * * WRKY WRKY cluster16_MYB cluster16_MYB * * * * cluster39 cluster39 ARR/GATA ARR/GATA * * TRP TRP cluster20_NAC cluster20_NAC * * RVE RVE * * TGA TGA ATHB/HAT ATHB/HAT * * HMG HMG STZ/SIZF STZ/SIZF * * * * cluster19_NAC cluster19_NAC SPL SPL * * * * cluster31_ARF cluster31_ARF ABI3/FUS3 ABI3/FUS3 * * * * CMTA3/FAR1 CMTA3/FAR1 PLT PLT * * bHLH/bZIP bHLH/bZIP FY FY * * E2FA E2FA ATHB/HAT ATHB/HAT * * CDC CDC BPC BPC * * cluster8_ERF cluster8_ERF Foxn1 Foxn1 * * * * ARF/ERF ARF/ERF CMTA3/FAR1 CMTA3/FAR1 * * * * cluster2_ERF cluster2_ERF bHLH/bZIP bHLH/bZIP * * * * cluster32_ARF cluster32_ARF cluster32_ARF cluster32_ARF * * LBD LBD STZ/SIZF STZ/SIZF * * BPC BPC * * SPL SPL * * ARF ARF * * ABI3/FUS3 ABI3/FUS3 Foxn1 Foxn1 * * * * cluster1_TCP cluster1_TCP LFY LFY * * cluster6_TCP cluster6_TCP cluster1_TCP cluster1_TCP * * LBD LBD cluster6_TCP cluster6_TCP * * * * CDC CDC PLT PLT * * * * E2FA E2FA REF/RAV REF/RAV * * * * GT-1/AGL GT-1/AGL cluster26_ARF cluster26_ARF * * bZIP/RAV bZIP/RAV cluster8_ERF cluster8_ERF * * ARF/ERF ARF/ERF * * cluster2_ERF cluster2_ERF * * 20 20 40 40 60 60 80 80 10 100 0 20 20 40 40 60 60 8 80 0 1 100 00 20 20 40 40 60 60 80 80 100 100 20 20 40 40 60 60 8 80 0 1 100 00 20 20 40 40 60 60 80 80 100 100 aBIOTECH Fig. 2 Subgenome-divergent A A regulation during wheat grain Leafs Leafs development A Motif density DAP5s DAP5s of 40 clustered motifs from DAP9s DAP9s JASPAR database (Castro- DAP5_9s DAP5_9s Mondragon et al. 2022)in DAP15s DAP15s eight diffACR clusters. DAP20s DAP20s B Ternary plot showing DAP15&20s DAP15&20s expression bias of the syntenic -2.5 -2.5 0 0 2.5 2.5 genes from three subgenomes. Z-score Z-score The balanced category: (Obeserved/expected) (Obeserved/expected) balanced expression of the triad genes. A dominant, B B B CD CD RNA-seq (gene expression) RNA-seq (gene expression) ATAC-seq (TSS region) ATAC-seq (TSS region) dominant and D dominant: Categor Category y higher expression from A, B or Balanced Balanced D than the other two A suppressed A suppressed orthologs. A suppressed, B B suppressed B suppressed suppressed and D suppressed: D suppressed D suppressed lower expression from A, B or A dominant A dominant D than the other two B dominant B dominant D dominant D dominant orthologs. C Sankey diagram showing expression dynamic of the syntenic genes during Leaf Leaf D DA AP P5 5 D DA AP P9 9 D DA AP P15 15 D DA AP P2 20 0 Leaf Leaf D DA AP P5 5 D DA AP P9 9 D DA AP P1 15 5 D DA AP P20 20 grain development. D Sankey E E Leaf Leaf DAP5 DAP5 DAP9 DAP9 DAP1 DAP15 5 DAP2 DAP20 0 diagram showing the B B B B B B B B B B chromatin accessibility 100 100 100 100 100 100 100 100 100 100 bZIP28 bZIP28 80 80 NA NAC019 C019 80 80 80 80 80 80 80 80 dynamics of the promoter NA NAC019 C019 GAMyb GAMyb NA NAC019 C019 60 60 60 60 60 60 60 60 60 60 NA NAC100 C100 10 100 NA NAC0 C0 C0 C019 C019 0 NA NAC019 C019 P PBF PBF regions (from 2 Kb upstream bZIP28 bZIP28 P28 P28 PBF PBF PBF PBF GAM GA G G G GAMyb GAMyb A AM AM AM 40 40 40 40 bZIP28 bZIP28 P28 P P28 P P P 8 8 8 40 40 40 40 40 40 bZIP28 bZIP28 8 8 8 8 bZIP28 bZIP28 IP IP28 P2 P 8 bZIP28 bZIP28 ZIP28 ZI ZIP28 28 8 8 8 PBF PBF NA NAC100 C100 N N N N N N N N N N N N NA NA A A A A AC C C C C C C019 C019 PBF PB PBF PB P PBF PBF F F GAM G G G GAMyb GAMyb AM Myb yb b NA N N NA A A A AC100 C100 N N N N NA NA A A A AC1 C C C C1 C C100 C100 N N N N NA NA A AC C100 C100 to 100 bp downstream of TSS) PBF PBF 20 20 20 20 20 20 20 20 GAMyb GAMyb NA NAC100 C100 RSR1 RSR1 20 20 RSR1 RSR1 GAMyb GAMyb GAMyb GAMyb RSR1 RSR1 RSR1 RSR1 RSR1 RSR1 RSR1 RSR1 of the syntenic genes during A A D D A A D D A A D D A A D D A A D D grain development. E The FF subgenome diversified expressions of key A A transcription factors. F Motif ** ** ** ** ** ** ** ** ** ** ** ** * * Leafs Leafs B B ** ** ** ** * * * * density of eight diffACR D D ** ** ** ** * * * * A A clusters in each subgenome (* ** ** * * DAP5s DAP5s B B indicate motif density [ 1) ** ** * * D D ** ** * * A A * * * * * * * * * * ** ** ** ** ** ** * * DAP9s DAP9s B B D D A A ** ** * * ** ** * * * * ** ** ** ** * * DAP5&9s DAP5&9s B B D D A A * * * * * * DAP15s DAP15s B B * * * * * * * * * * * * * * * * * * ** ** ** ** ** ** ** ** D D * * * * * * * * * * * * * * * * * * ** ** ** ** ** ** ** ** A A * * * * * * DAP20s DAP20s B B * * * * * * ** ** ** ** * * * * ** ** ** ** * * * * * * * * D D * * * * * * ** ** ** ** * * * * ** ** ** ** * * * * * * * * A A * * * * * * * * * * ** ** * * * * ** ** * * * * * * * * DAP15&20s DAP15&20s B B * * * * ** ** * * * * D D * * * * ** ** * * * * A A * * * * * * * * * * * * * * Alls Alls B B * * * * * * * * * * * * D D * * * * * * * * * * * * -3 -3 0 0 3 3 Z-score Z-score (Obeserved/expected) (Obeserved/expected) promoter accessible regions in all samples (Supple- binding divergence among the three subgenomes, we mentary Fig. 4C). calculated the densities of TF-binding motifs in different The diversity of CREs was considered as key reasons grain specific ACRs among three subgenomes. Although of transcriptomic expression differentiations (Buen- the number of diffACRs in three subgenomes did not rostro et al. 2013; Lu et al. 2017). The previous reported show obvious differences (Supplementary Fig. 4D), a key TFs controlling grain development showed dynamic large proportion of TF-binding sites showed bias subgenome expression differences in leaf, DAP5, DAP9, enrichments among the three subgenomes (Fig. 2F). For DAP15 and DAP20 (Fig. 2E). To identify the specific TF- example, E2FA, CDC, LBD, BPC, ARF and Foxn1 binding The Author(s) 2023 20 20 40 40 60 60 80 80 100 100 20 20 40 40 60 60 80 80 100 100 2 20 0 40 40 60 60 80 80 100 100 20 20 40 40 60 60 80 80 100 100 20 20 4 40 0 60 60 80 80 100 100 aBIOTECH sites showed higher enrichment on A subgenome in development and may contribute to the improved grain leaf-specific ACRs. MYB, NAC and TGA binding sites yield and nutritional content after polyploidization. showed higher enrichment on A subgenome in DAP9s and DAP5&9s diffACRs. bZIP, BPC and KUA TF-binding Starch and gluten regulatory networks mediated sites showed higher enrichment on A subgenome in by key TFs DAP15, DAP20 and DAP15&20s groups. ARF, MYB and TGA binding sites were subgenome divergently enriched During the filling stages, expression of genes involved in in DAP15s, DAP20s and AllGs groups. The densities of the starch and gluten biosynthesis pathways were TF-binding sites were also varied during grain devel- gradually increased. However, the exact cis-regulatory opment. For example, NAC, KUA and MYB binding motifs elements responsible for the tissue-specific expression were enriched in A subgenome at the DAP9s, while bias were still elusive. We picked up the genes related with enriched in B and D subgenomes at DAP15s (Fig. 2E). gluten protein accumulation and starch biosynthesis Gene Ontology (GO) enrichment analysis showed that regulation including HMW-GS, LMW-GS, alpha-gliadin, the balanced and unbalanced genes among subgenomes gamma-gliadin and AGPS1, AGPL1, GBSSI, SBEI, SBEII, during seed development were enriched in different SSII, and SuSy in starch biosynthesis (Supplemental Data pathways, the balanced genes were enriched in DNA Set4). We found that the promoter accessibilities of replication, gene expression and protein folding process, genes involved in starch biosynthesis were higher in while A suppressed genes were mainly enriched in seeds than leaves and increased after DAP5, while those vitamin biosynthetic process, thiamine biosynthetic of glutenin and gliadin genes were increased after DAP9 process and lipid catabolic process, B-suppressed genes (Fig. 3A), which is consistent with the transcriptomic were mainly enriched in carbohydrate homeostasis expression patterns (Fig. 3B). Interestingly, the pro- process, carbohydrate phosphorylation process and moter accessibilities of genes in DAP20 were starting RNA modification process, D suppressed genes were decrease while those genes were still highly expressed mainly enriched in organic substance metabolic process (Fig. 3A and B). Meanwhile, for some genes highly and nitrogen compound biosynthetic process, A domi- expressed after DAP9, higher promoter accessibilities nant genes were mainly enriched in phosphatidylinosi- could be found at DAP5 (Fig. 3A and B), indicating the tol biosynthetic process and carbohydrate catabolic changes of promoter accessibilities are prior to the process, B dominant genes were mainly enriched in transcriptomic expression changes. Taken together, carbohydrate biosynthetic process and metal ion these results indicated that the chromatin accessibilities transport process, D dominant genes were mainly were tightly associated with the expression changes of those grain-filling-related genes. enriched in protein phosphorylation and aminoglycan metabolic process. These enriched GO terms were con- To find the key TFs controlling grain development, we sistent with the grain development processes (Supple- further calculated the densities of TF-binding motifs in mentary Fig. 5). To prove the specific binding of TFs in the promoters of genes related with gluten protein different tissues, we combined the published genome- accumulation and starch biosynthesis, which suggested wide profiling of TFBSs (transcription factor-binding that NAC, bZIP and MYB family TFs play important roles sites) in common wheat (DAP-seq data) (Zhang et al. (Supplementary Fig. 7). Consistently, the key TFs rela- 2022) and our ATAC-seq data of different tissues to ted with grain development, including TaNAC019, Tab- better identify the potential in vivo TFBSs in different ZIP28, and TaGAMyb showed increased promoter tissues. The results showed that more that 60% TFs’ accessibilities and transcriptomic expression patterns in (NAC6A-1, NAC6B-1, NAC-6D-1 and MYB-7A-3) binding grains (Fig. 3C and Supplementary Fig. 8). We then sites were different among various tissues (Supple- constructed a regulatory network based on the chro- mentary Fig. 6A). In addition, we analyzed the sub- matin accessibility data and expression pattern of those genome-divergent regulation by TFs (NAC6A-1, NAC6B- key TFs and genes involved in starch biosynthesis and 1, NAC-6D-1 and MYB-7A-3), and found that nearly 80% protein accumulations. The results showed that the traids showed subgenome-divergent regulations (Sup- three copies of TaSPA were involved in the networks in plementary Fig. 6B). These results indicated tissue- DAP9 and DAP15; TaNAC019 and TaNAC100 in DAP15 specific and subgenome-divergent regulation by multi- and DAP20; while the others mostly in all the stages ple TFs during wheat grain development. The high after DAP9 (Fig. 3D). Subgenome specific regulation variations of TF-binding sites enrichment bias among could be found with FUSCA3, of which the copy from D subgenomes at different stages suggested the three subgenome were involved in networks in DAP9 and subgenomes worked coordinately during wheat grain DAP15, while the other two in all the three stages. Similarly, the copy of TaGAMYB from D subgenome were The Author(s) 2023 aBIOTECH -200 2 36 TaSBEIIa TaSuSy2 TaSBEIb A B C Scaled TIS log10(TPM+1) TraesCS2A02G293400 TraesCS2B02G194200 TraesCS7A02G549300 TaAGPL1-A [0 - 77] [0 - 48] [0 - 99] TaAGPL1-B Leaf TaAGPL1-D [0 - 77] [0 - 99] [0 - 48] DAP5 TaAGPS1a-A TaAGPS1a-B [0 - 48] [0 - 77] [0 - 99] DAP9 TaAGPS1a-D TaAGPS2-D [0 - 77] [0 - 99] [0 - 48] DAP15 TaGBSSI-A TaGBSSI-D [0 - 48] [0 - 77] [0 - 99] DAP20 TaISA1-B TaISA1-D [0 - 249875] [0 - 668244] [0 - 25185] Leaf TaISA2-B [0 - 25185] [0 - 249875] [0 - 668244] TaPUL-A DAP5 TaPUL-B TaPUL-D [0 - 249875] [0 - 668244] [0 - 25185] DAP9 TaSBEIa-A TaSBEIa-D [0 - 249875] [0 - 668244] [0 - 25185] DAP15 TaSBEIb-A TaSBEIb-B [0 - 249875] [0 - 25185] [0 - 668244] DAP20 TaSBEIc-D TaSBEIIa-A TaSBEIIa-B TaGBSSI TabZIP28 TaNAC019 TaSBEIIa-D TraesCS7A02G070100 TraesCS2D02G146100 TraesCS3B02G092800 TaSBEIIb-A TaSBEIIb-B [0 - 28] [0 - 125] TaSBEIIb-D [0 - 149] Leaf TaSSIc-A [0 - 125] [0 - 28] [0 - 149] TaSSIc-D DAP5 TaSSIIa-A [0 - 125] [0 - 149] [0 - 28] TaSSIIa-B DAP9 TaSSIIa-D [0 - 125] [0 - 149] [0 - 28] TaSSIIIa-D DAP15 TaSSIIIb-A TaSSIIIb-B [0 - 125] [0 - 149] [0 - 28] DAP20 TaSuSy1-A [0 - 25580] [0 - 383539] [0 - 53686] TaSuSy1-B Leaf TaSuSy1-D [0 - 25580] TaSuSy2-A [0 - 383539] [0 - 53686] DAP5 TaSuSy2-B [0 - 25580] TaSuSy2-D [0 - 383539] [0 - 53686] DAP9 TaSuSy4-A TraesCS7D02G064300TraesCS7D02G117800 TaSuSy4-D [0 - 25580] [0 - 383539] [0 - 53686] DAP15 TaSuSy5-A TaSuSy5-B [0 - 25580] DAP20 [0 - 383539] [0 - 53686] TaSuSy5-D Omega-gliadin TaNAC100 D TaNAC019 Interaction in B DAP9 B DAP15 DAP20 D Gamma-gliadin FUSCA3 SPR Alpha-gliadin B HMW-glutenin TaGAMYB TabZIP28 LMW-glutenin E15 E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E1 E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E1 1 1 15 15 1 1 1 15 15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 15 15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 E9 E E E5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 9 9 9 9 9 9 9 9 9 9 Glutenin Gliadin Starch synthesis A PBF-D TaSPA Fig. 3 Construction of regulatory networks for starch and gluten biosynthesis by integration the chromatin accessibility and expression data. A and B Promoter chromatin accessibility dynamics (A) and expression dynamics (B) of genes involved in glutenin and starch accumulation in leaf and grain samples. Regions from 2 Kb upstream to 100 bp downstream of transcriptional start sites were considered as promoters. C Genome browser tracks showing the expression and promoter chromatin accessibility patterns of starch synthesis genes (TaSBEIIa, TaSuSy2, TaSBEIb and TaGBSSI) and key TF genes (TaZIP28 and TaNAC019). D Regulatory networks involving the starch and gluten biosynthesis genes and key TFs. TF families are shown as circles with different colors. Starch biosynthesis genes are shown as grey square; Gliadin genes as orange triangle; Glutenin genes as blue diamond. TF-to-gene edges are shown with three colors represent three grain phases. DAP9, green; DAP15 blue; DAP20, pink in DAP9 and A subgenome copy in DAP9 and DAP15 and protein synthesis. To further verify the predicted while B subgenome copy in all the three stages. The transcriptional regulation network, we conducted the results also revealed most of these key TFs were dual-luciferase transcriptional activity assay (Fig. 4A–B involved in both starch synthesis and protein accumu- and Supplemental Data Set5). The assay showed that lations, indicating a complex network between starch the reporter activities driven by the promoters of The Author(s) 2023 Leaf DAP5 DAP9 DAP15 DAP20 Leaf DAP5 DAP9 DAP15 DAP20 RNA-seq ATAC-seq RNA-seq ATAC-seq aBIOTECH TaHWM-1Dpro TaHWM TaAGPL1-1Bpro TaAGPL1 TaISA2-1Bpro TaISA2 TF motif TF motif TF motif 2 2 1 1 cluster_19_NAC cluster_19_NAC cluster_19_NAC 1 0 0 123456789 1011121314 1 2 3 4 5 6 7 8 9 1011121314 123456789 1011121314 2 2 1 1 cluster_20_NAC cluster_20_NAC cluster_20_NAC 0 0 112 24 112 24 112 24 2 2 Position 1 1 cluster_12_MYB cluster_12_MYB 0 0 111 22 111 22 2 2 1 1 cluster_16_MYB cluster_16_MYB 018 16 018 16 Position Position BC TaAGPL1-1Bpro Effector 0.04 *** 0.004 TaHWM-1Dpro 35S GFP NOS *** 0.03 0.003 0.02 35S MYB_B-GFP NOS 0.002 0.01 TaISA2-1Bpro 35S NAC100_A-GFP NOS TaHWM-1Dpro 0.001 ** *** 0.0008 0.003 35S NAC100_D-GFP NOS 0.0006 0.002 Reporter 0.0004 0.001 0.0002 TaHWM-1Dpro LUC 35S REN NOS 0.0000 0.000 TaAGPL1-1Bpro LUC 35S REN NOS TaISA2-1Bpro LUC 35S REN NOS TaISA2-1Bpro TaISA2-1Bpro TaAGPL1-1Bpro TaAGPL1-1Bpro 0.0015 0.0020 0.015 0.015 ** **** 0.0015 **** *** 0.0010 0.010 0.010 *** ** *** 0.0010 * **** *** 0.0005 0.005 0.005 **** 0.0005 0.000 0.0000 0.0000 0.000 Fig. 4 Dual-luciferase transcriptional activity assay to assess the capability of key TFs to transactivate predicted target gene expression. A The distribution of key TF-binding motifs in the promoters of starch and protein synthesis genes. B Schematic diagrams of the effector and reporter constructs for the dual-luciferase transcriptional activity assay. C Reporter assay showed that NAC100-A or NAC100- D regulates both starch and protein synthesis pathways. D NAC100 and MYB-B coordinately regulate TaISA2 and TaAGPL1 genes. LUC/ REN indicates the signal ratio of LUC (firefly luciferase) to REN (Renilla reniformis luciferase) activity. Data are represented as means ± SD from three replicates. (Student’s t-test; *P \ 0.05, **P \ 0.01) TaAGPL1-1B, TaHWM-1D and TaISA2-1B were markedly (TaAGPL1-1B and TaISA2-1B) (Fig. 4D). In summary, we activated in the presence of NAC100-A or NAC100-D, have explored the chromatin accessibility landscapes suggested that NAC100 regulates both starch and pro- and identified the key TFs and regulatory network that tein synthesis pathways (Fig. 4C). Moreover, NAC100 were involved in the starch and protein biosynthesis and MYB could coordinately regulate same targets during wheat grain development. The Author(s) 2023 GFP+GFP GFP+NAC100_A GFP+MYB_B GFP NAC100_A+MYB_B NAC100_A GFP NAC100_A GFP+GFP GFP+NAC100_D GFP+MYB_B NAC100_D+MYB_B GFP NAC100_D GFP GFP+GFP NAC100_D GFP+NAC100_A GFP+MYB_B NAC100_A+MYB_B GFP+GFP GFP+NAC100_D GFP+MYB_B NAC100_D+MYB_B LUC/REN Bits Bits LUC/REN Bits Bits Bits Bits LUC/REN LUC/REN LUC/REN LUC/REN Bits Bits Bits Bits aBIOTECH DISCUSSION METHODS Wheat grain development is critical for the quality and Plant material and growth conditions yield of wheat, a series of genes regulated storage protein and starch biosynthesis have been discovered in Wheat (Triticum aestivum; BBAADD, 2n = 6x = 42) the previous study (Xiao et al. 2022). Transcriptional cultivar AK58 was transferred to the greenhouse with regulation of gluten and starch genes involves complex conditions day/night cycle set at 16 h (22 C)/8 h network between cis- and trans-acting factors, but little (19 C) after vernalization for 1 mouth in cold chamber is known about this regulation. The accessible chro- (4 C). Wheat endosperm was excised by removing the matin regions often associated with TFs related CREs seed embryo at 5 and 9 days after fertilization. Wheat (Lu et al. 2017; Marand et al. 2021; Pei et al. 2022). seed was collected at 15 and 20 days after fertilization. Here, we constructed a chromatin accessibility atlas for leaf and four different stages of seed during grain RNA-seq development. This regulatory atlas enabled specific ACR and TF identification at different phase of wheat grain The collected wheat endosperm and seed were flash- development. We found NAC, bZIP and MYB binding frozen with liquid N2 immediately. Total RNA was TM motifs are highly enriched in the grain specific ACRs as extracted with TRIzol Reagent (Invitrogen; previous reported (Gao et al. 2021; Guo et al. 2015; Liu 15,596–026) following the manufacturer’s instructions. et al. 2020; Song et al. 2020). More importantly, some RNA-Seq libraries were constructed and sequenced by novel TFs, including REF, Dof and HMG were also dis- Berry Genomics (Beijing, China), the libraries were covered highly enriched in the these specific ACRs. sequenced on the Illumina NovaSeq 6000 platform and Further studies on TFs in those families will help to produced 150-bp paired-end reads. understand the molecular basis of grain development. Chromatin accessibilities’ changes are associated ATAC-seq with genes’ differential expression. In this study, we combined ATAC-seq with RNA-seq and constructed the ATAC-seq was performed as described previously (Lu dynamic regulatory network during wheat grain devel- et al. 2017). For each sample, approximately 1 g flash- opment, and further validated the regulatory network frozen wheat samples were chopped with a razor blade by dual-luciferase transcriptional activity assays. We in 1 mL ice-prechilled lysis buffer (15 mM Tris–HCl pH found that the phase-specific ACRs are associated with 7.5, 20 mM NaCl, 80 mM KCl, 0.5 mM spermine, 5 mM phase-specific DEGs, and the function of phase-specific 2-Mercaptoethanol, 0.2% TritonX-100). The chopped genes is related to specific the biological process of slurry containing crude nuclei extract was filtered twice specific tissue. The seed-specific genes were highly through a 40 lm filter. The crude nuclei were stained enriched in carbohydrate and sucrose metabolic pro- with DAPI (sigma, catalog number: D9542) and loaded cess, and seed development processes. We also found to a flow cytometer (BD FACSCanto) for selected. Nuclei the changes of promoter accessibilities of starch and pellets were obtained after centrifuged and washed with protein biosynthesis genes were prior to their expres- Tris-Mg buffer (10 mM Tris–HCl pH 8.0, 5 mM MgCl ). sion changes, indicating those RNAs were likely accu- The obtained nuclei were incubated with 3.5 lLTn5 mulated rather than rapidly synthesized. Besides, we transposomes in 40 lL TTBL buffer (TruePrep DNA found the homeolog dynamic transcriptomic expression Library Prep Kit V2 for Illumina, Vazyme Biotech co., ltd, divergence during wheat grain development (Xiang TD501) for each sample at 37 C for 30 min without et al. 2019) are associated with the divergence of pro- rotation. We purified the integration products with a TM moter accessibilities and the TF-binding sites in the NEB Monarch DNA Cleanup Kit (T1030S) and then three subgenomes at different phases of wheat grain amplified for 10–13 cycles using the NEBNext Ultra II development. Therefore, those accessible regions in Q5 master mix (M0544L). PCR cycles were determined different phases could help to accurately discover the as described previously (Lu et al. 2017). Amplified key cis- and trans-acting factors during wheat grain libraries were purified with Hieff NGS DNA Selection development. In summary, our data provided a tran- Beads (Yeasen, 12601ES03) to remove free index scriptional regulation resource to understand the com- primers. plex networks and identified the key regulatory elements and TFs underlying wheat grain development. The Author(s) 2023 aBIOTECH Supplementary InformationThe online version contains Processing of RNA-seq data supplementary material available at https://doi.org/10.1007/ s42994-023-00095-8. Raw reads were preprocessed by fastp for filtering low- quality reads and adapter trimming. The remaining Acknowledgements We would like to acknowledge Ting Li from the flow cytometry core of Institute of Genetics and Develop- reads were aligned to the wheat reference genome using mental Biology, Chinese Academy of Sciences for supporting the hisat2 with default parameters. All aligned reads were nuclear sorting. This project was financially supported by the sorted by SAMtools v1.3.1 (Swiezewski et al. 2009). To Outstanding Young Scientist Foundation of NSFC (Overseas), the compare the expression from different phases, TPM was Central Public-interest Scientific Institution Basic Research Found (S2022ZD02), the Fundamental Research Funds from Institute of calculated by TPMCalculator. Differential expressed Crop Sciences, Chinese Academy of Agricultural Sciences genes (DEGs) were defined as adjusted P-value \ 0.01 (S2020YC07 and S2021YC03) and CAAS Agricultural Science and and log2(FoldChange) [ 1, which determined by Technology Innovation Program, China (CAAS-ZDRW202002). DESeq2. Author contributions ZL conceived and designed the project. HP and YSL performed ATAC-seq experiments. YHL, PL, and XR were Processing of ATAC-seq data responsible for wheat materials. HP, YSL and JLZ analyzed the data. ZL, YSL, YHL, and PL contributed to project discussion. HP Raw reads were trimmed with fastp with default and ZL wrote the manuscript draft, and YSL, YHL and PL revised it. parameters. Trimmed reads were aligned to the wheat reference genome using Bowtie v2.2.4 with the follow- Data availability All the raw sequencing data generated during the current study are available in the Gene Expression Omnibus ing parameters: ‘‘bowtie2-X 1000 –very-sensitive’’. (GEO) database (https://www.ncbi.nlm.nih.gov/geo) under Aligned reads were sorted using SAMtools v1.3.1 accession number GSE214739. (Swiezewski et al. 2009) and clonal duplicates were removed using Picard version v2.16.0 http://broad Declarations institute.github.io/picard/). Peak calling was performed Conflict of interest The authors declare that they have no con- as described previously. A black list was first generated flict of interest. using the control samples as input with MACS2. The raw peaks were first called with MACS2 with the following Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, parameter ‘‘–keep-dup all –nomodel –extsizes 150 –shift adaptation, distribution and reproduction in any medium or for- -75’’ (Zhang et al. 2008); then split into 150 bp bins mat, as long as you give appropriate credit to the original with 50 bp overlapping; bins passed filtering by a Tn5 author(s) and the source, provide a link to the Creative Commons integration site density cutoff were selected and merged licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s with bedtools with ‘‘– d 150’’; peaks overlapped with the Creative Commons licence, unless indicated otherwise in a credit blacklist peaks, or homologous to plant organelles DNA line to the material. If material is not included in the article’s (NCBI) were discard and the rest were considered as Creative Commons licence and your intended use is not permitted high-quality ACRs. Bigwig files were used to visualize by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To the peaks in the integrative genomics viewer. view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. Dual-luciferase reporter gene assay Dual-luciferase reporter gene assay was performed in N. benthamiana leaves by Agrobacterium tumefaciens-me- References diated transient expression system. The promoters of TaAGPL1-1B, TaHWM-1D and TaISA2-1B were amplified Adamski NM, Simmonds J, Brinton JF, Backhaus AE, Chen Y, Smedley M, Hayta S, Florio T, Crane P, Scott P et al (2021) from AK58 genomic DNA and inserted into the pGreenII Ectopic expression of Triticum polonicum VRT-A2 underlies 0800-LUC vector as reporter plasmids (Hellens et al. elongated glumes and grains in hexaploid wheat in a dosage- 2005). The CDS of NAC100-A, NAC100-D and MYB-B dependent manner. Plant Cell 33:2296–2319 were cloned in frame into the pEG1300-GFP vector as Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ (2013) Transposition of native chromatin for fast and sensitive effector plasmids. N. benthamiana plants at the six-leaf epigenomic profiling of open chromatin, DNA-binding pro- stage were co-infiltrated with A. tumefaciens strain teins and nucleosome position. Nat Methods 10:1213–1218 GV3101 harboring different combinations of these Castro-Mondragon JA, Riudavets-Puig R, Rauluseviciute I, Lemma plasmids. The dual-luciferase reporter assay system RB, Turchi L, Blanc-Mathieu R, Lucas J, Boddie P, Khan A, Manosalva Perez N et al (2022) JASPAR 2022: the 9th release (Yeasen, China) was used to measure LUC and REN of the open-access database of transcription factor binding activities and three independent replications were profiles. Nucl Acids Res 50:D165–D173 conducted. The Author(s) 2023 aBIOTECH Chen S, Lake BB, Zhang K (2019) High-throughput sequencing of Lu Z, Marand AP, Ricci WA, Ethridge CL, Zhang X, Schmitz RJ the transcriptome and chromatin accessibility in the same (2019) The prevalence, evolution and chromatin signatures cell. Nat Biotechnol 37:1452–1457 of plant regulatory elements. Nat Plants 5:1250–1259 Deplancke B, Alpern D, Gardeux V (2016) The genetics of Marand AP, Chen ZL, Gallavotti A, Schmitz RJ (2021) A cis- transcription factor DNA binding variation. Cell 166:538–554 regulatory atlas in maize at single-cell resolution. Cell Gao Y, An K, Guo W, Chen Y, Zhang R, Zhang X, Chang S, Rossi V, Jin 184:3041 F, Cao X et al (2021) The endosperm-specific transcription Pajoro A, Madrigal P, Muino JM, Matus JT, Jin J, Mecchia MA, factor TaNAC019 regulates glutenin and starch accumulation Debernardi JM, Palatnik JF, Balazadeh S, Arif M et al (2014) and its elite allele improves wheat grain quality. Plant Cell Dynamics of chromatin accessibility and gene regulation by 33:603–622 MADS-domain transcription factors in flower development. Gu A, Hao P, Lv D, Zhen S, Bian Y, Ma C, Xu Y, Zhang W, Yan Y Genome Biol 15:R41 (2015) Integrated proteome analysis of the wheat embryo Pei H, Teng W, Gao L, Gao H, Ren X, Liu Y, Jia J, Tong Y, Wang Y, Lu Z and endosperm reveals central metabolic changes involved in (2022) Low-affinity SPL binding sites contribute to sub- the water deficit response during grain development. J Agric genome expression divergence in allohexaploid wheat. Sci Food Chem 63:8478–8487 China Life Sci 1–16 Guo W, Yang H, Liu Y, Gao Y, Ni Z, Peng H, Xin M, Hu Z, Sun Q, Yao Y Ramirez-Gonzalez RH, Borrill P, Lang D, Harrington SA, Brinton J, (2015) The wheat transcription factor TaGAMyb recruits Venturini L, Davey M, Jacobs J, van Ex F, Pasha A et al (2018) histone acetyltransferase and activates the expression of a The transcriptional landscape of polyploid wheat. Science high-molecular-weight glutenin subunit gene. Plant J 361:6403 84:347–359 Ren J, Jiang Z, Li W, Kang X, Bai S, Yang L, Li S, Zhang D (2022) Hellens RP, Allan AC, Friel EN, Bolitho K, Grafton K, Templeton MD, Characterization of glutenin genes in bread wheat by third- Karunairetnam S, Gleave AP, Laing WA (2005) Transient generation RNA sequencing and the development of a Glu- expression vectors for functional genomics, quantification of 1Dx5 marker specific for the extra cysteine residue. J Agric promoter activity and RNA silencing in plants. Plant Methods Food Chem 70:7211–7219 1:13 Rodgers-Melnick E, Vera DL, Bass HW, Buckler ES (2016) Open Huang L, Tan H, Zhang C, Li Q, Liu Q (2021) Starch biosynthesis in chromatin reveals the functional maize genome. Proc Natl cereal endosperms: An updated review over the last decade. Acad Sci USA 113:E3177-3184 Plant Commun 2:100237 She KC, Kusano H, Koizumi K, Yamakawa H, Hakata M, Imamura T, Kang G, Liu G, Peng X, Wei L, Wang C, Zhu Y, Ma Y, Jiang Y, Guo T Fukuda M, Naito N, Tsurumaki Y, Yaeshima M et al (2010) A (2013) Increasing the starch content and grain weight of novel factor FLOURY ENDOSPERM2 is involved in regulation common wheat by overexpression of the cytosolic AGPase of rice grain size and starch quality. Plant Cell 22:3280–3294 large subunit gene. Plant Physiol Biochem 73:93–98 Song Y, Luo G, Shen L, Yu K, Yang W, Li X, Sun J, Zhan K, Cui D, Liu D Kouzarides T (2007) Chromatin modifications and their function. et al (2020) TubZIP28, a novel bZIP family transcription Cell 128:693–705 factor from Triticum urartu, and TabZIP28, its homologue Li Z, Wang M, Lin K, Xie Y, Guo J, Ye L, Zhuang Y, Teng W, Ran X, from Triticum aestivum, enhance starch synthesis in wheat. Tong Y et al (2019) The bread wheat epigenomic map reveals New Phytol 226:1384–1398 distinct chromatin architectural and evolutionary features of Sun F, Liu X, Wei Q, Liu J, Yang T, Jia L, Wang Y, Yang G, He G (2017) functional genetic elements. Genome Biol 20:139 Functional characterization of TaFUSCA3, a B3-superfamily Liu A, Bergmann DC (2021) How to build a crop plant: defining transcription factor gene in the wheat. Front Plant Sci 8:1133 the cis-regulatory landscape of maize. Cell 184:2804–2806 Swiezewski S, Liu F, Magusin A, Dean C (2009) Cold-induced Liu G, Wu Y, Xu M, Gao T, Wang P, Wang L, Guo T, Kang G (2016) silencing by long antisense transcripts of an Arabidopsis Virus-induced gene silencing identifies an important role of Polycomb target. Nature 462:799–802 the TaRSR1 transcription factor in starch synthesis in bread Tian H, Li Y, Wang C, Xu X, Zhang Y, Zeb Q, Zicola J, Fu Y, Turck F, Li wheat. Int J Mol Sci 17(10):1557 L et al (2021) Photoperiod-responsive changes in chromatin Liu J, Chen Z, Wang Z, Zhang Z, Xie X, Wang Z, Chai L, Song L, Cheng accessibility in phloem companion and epidermis cells of X, Feng M et al (2021) Ectopic expression of VRT-A2 underlies Arabidopsis leaves. Plant Cell 33:475–491 the origin of Triticum polonicum and Triticum petropavlovskyi Wan Y, Poole RL, Huttly AK, Toscano-Underwood C, Feeney K, with long outer glumes and grains. Mol Plant Welham S, Gooding MJ, Mills C, Edwards KJ, Shewry PR et al 14(9):1472–1488 (2008) Transcriptome analysis of grain development in Liu L, He ZH, Yan J, Zhang Y, Xia XC, Pena RJ (2005) Allelic hexaploid wheat. BMC Genom 9:121 variation at the Glu-1 and Glu-3 loci, presence of the 1B.1R Wang M, Li Z, Zhang Y, Zhang Y, Xie Y, Ye L, Zhuang Y, Lin K, Zhao F, translocation, and their effects on mixographic properties in Guo J et al (2021a) An atlas of wheat epigenetic regulatory Chinese bread wheats. Euphytica 142:197–204 elements reveals subgenome divergence in the regulation of Liu Y, Han C, Deng X, Liu D, Liu N, Yan Y (2018) Integrated development and stress responses. Plant Cell 33:865–881 physiology and proteome analysis of embryo and endosperm Wang X, Aguirre L, Rodriguez-Leal D, Hendelman A, Benoit M, highlights complex metabolic networks involved in seed Lippman ZB (2021b) Dissecting cis-regulatory control of germination in wheat (Triticum aestivum L.). J Plant Physiol quantitative trait variation in a plant stem cell circuit. Nat Plants 7:419–427 229:63–76 Liu Y, Hou J, Wang X, Li T, Majeed U, Hao C, Zhang X (2020) The Wang Y, Hou J, Liu H, Li T, Wang K, Hao C, Liu H, Zhang X (2019) NAC transcription factor NAC019-A1 is a negative regulator TaBT1, affecting starch synthesis and thousand kernel weight, of starch synthesis in wheat developing endosperm. J Exp Bot underwent strong selection during wheat improvement. J Exp 71:5794–5807 Bot 70:1497–1511 Lu Z, Hofmeister BT, Vollmers C, DuBois RM, Schmitz RJ (2017) Xiang DQ, Quilichini TD, Liu ZY, Gao P, Pan YL, Li Q, Nilsen KT, Combining ATAC-seq with nuclei sorting for discovery of cis- Venglat P, Esteban E, Pasha A et al (2019) The Transcriptional regulatory regions in plant genomes. Nucl Acids Res 45:e41 landscape of polyploid wheats and their diploid ancestors The Author(s) 2023 aBIOTECH during embryogenesis and grain development. Plant Cell subgenome-convergent and -divergent transcription in com- 31:2888–2911 mon wheat. Nat Commun 13:6940 Xiao J, Liu B, Yao Y, Guo Z, Jia H, Kong L, Zhang A, Ma W, Ni Z, Xu S Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al (2022) Wheat genomic study for genetic improvement of Nusbaum C, Myers RM, Brown M, Li W et al (2008) Model- traits in China. Sci China Life Sci 65:1718–1775 based analysis of ChIP-Seq (MACS). Genome Biol 9:R137 Yamamori M, Fujita S, Hayakawa K, Matsuki J, Yasui T (2000) Zhou Y, Zhao X, Li Y, Xu J, Bi A, Kang L, Xu D, Chen H, Wang Y, Wang Genetic elimination of a starch granule protein, SGP-1, of YG et al (2020) Triticum population sequencing provides wheat generates an altered starch with apparent high insights into wheat adaptation. Nat Genet 52:1412–1422 amylose. Theor Appl Genet 101:21–29 Zhu J, Fang L, Yu J, Zhao Y, Chen F, Xia G (2018) 5-Azacytidine Yang C, Zhao L, Zhang H, Yang Z, Wang H, Wen S, Zhang C, Rustgi S, treatment and TaPBF-D over-expression increases glutenin von Wettstein D, Liu B (2014) Evolution of physiological accumulation within the wheat grain by hypomethylating the responses to salt stress in hexaploid wheat. Proc Natl Acad Glu-1 promoters. Theor Appl Genet 131:735–746 Sci USA 111:11882–11887 Zhang Y, Li Z, Liu J, Zhang Y, Ye L, Peng Y, Wang H, Diao H, Ma Y, Wang M et al (2022) Transposable elements orchestrate The Author(s) 2023 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png aBIOTECH Springer Journals

Chromatin accessibility landscapes revealed the subgenome-divergent regulation networks during wheat grain development

Loading next page...
 
/lp/springer-journals/chromatin-accessibility-landscapes-revealed-the-subgenome-divergent-dMylATJpMW
Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2023
ISSN
2096-6326
eISSN
2662-1738
DOI
10.1007/s42994-023-00095-8
Publisher site
See Article on Publisher Site

Abstract

aBIOTECH https://doi.org/10.1007/s42994-023-00095-8 aBIOTECH RESEARCH ARTICLE Chromatin accessibility landscapes revealed the subgenome- divergent regulation networks during wheat grain development 1 1 1 1 1 Hongcui Pei , Yushan Li , Yanhong Liu , Pan Liu , Jialin Zhang , 1 1& Xueni Ren , Zefu Lu Institute of Crop Sciences, Chinese Academy of Agricultural Sciences, Beijing 100081, China Received: 30 September 2022 / Accepted: 22 January 2023 Abstract Development of wheat (Triticum aestivum L.) grain mainly depends on the processes of starch synthesis and storage protein accumulation, which are critical for grain yield and quality. However, the regulatory network underlying the transcriptional and physiological changes of grain development is still not clear. Here, we combined ATAC-seq and RNA-seq to discover the chromatin accessibility and gene expression dynamics during these processes. We found that the chromatin accessibility changes are tightly asso- ciated with differential transcriptomic expressions, and the proportion of distal ACRs was increased gradually during grain development. Specific transcription factor (TF) binding sites were enriched at different stages and were diversified among the 3 subgenomes. We further predicted the potential interactions between key TFs and genes related with starch and storage protein biosynthesis and found different copies of some key TFs played diversified roles. Overall, our findings have provided numerous resources and illustrated the regulatory network during wheat grain development, which would shed light on the improvement of wheat yields and qualities. Keywords Wheat, Chromatin accessibility, Subgenome-divergence, Regulatory network, Grain development INTRODUCTION development is important to understand the molecular basis of wheat yield and quality and discover novel key Wheat (Triticum aestivum L.) is one of the most regulatory factors. important staple crops for global food security, provid- Starch and protein are the main nutrient components ing more than a fifth of the calories and protein con- of wheat grain. A variety of starch biosynthesis genes sumed by humans (http://faostat.fao.org). Wheat grains have been characterized, including ADP-glucose are mainly composed of pericarp, endosperm and pyrophosphorylases (AGPases), ADP-glucose (ADPG) embryo, of which the endosperm is the main source of transporter, granule-bound starch synthases (GBSSs), white flour. The endosperm is a storage organ that starch synthases (SSs), Starch branching enzymes accumulates a large quantity of starch and proteins (SBEs), debranching enzymes (DBEs), starch/a-glucan during grain development (Gu et al. 2015; Liu et al. phosphorylases (PHOs), disproportionating enzymes 2018; She et al. 2010). Therefore, elucidation of the (DPEs), and protein targeting to starch (PTST), which transcriptional regulatory networks during wheat grain coordinately work in starch biosynthesis during wheat grain development (Huang et al. 2021; Kang et al. 2013; Wang et al. 2019; Yamamori et al. 2000). Starch is Hongcui Pei and Yushan Li have contributed equally to this work. synthesized by plants using the enzyme AGPases that reacts glucose 1-phosphate with ATP to form ADPG, the & Correspondence: luzefu@caas.cn (Z. Lu) The Author(s) 2023 aBIOTECH ADPG were subsequent synthetized to amylose and RESULTS amylopectin with different enzymes. Most of the storage proteins in grains are gluten, including high molecular Chromatin accessibility changes are associated weight glutenin subunits (HMW-GSs), low molecular with differential transcriptomic expressions weight glutenin subunits (LMW-GSs) and gliadins, during grain development which influence the dough property (Liu et al. 2005; Ren et al. 2022). The expressions of those key genes in Wheat grains are mainly composed of endosperm and starch and/or gluten biosynthesis pathways are tightly embryo, we performed ATAC-seq and RNA-seq with endosperm or seeds in 5, 9, 15, and 20 days after pol- related with grain yield and quality. A group of tran- scriptional factors, such as TaRSR1 and TabZIP28, were lination (DAP) to explore the chromatin accessibility found to be involved in regulating the expression of and transcriptomic expression dynamics during wheat genes related to starch synthesis in wheat (Liu et al. grain development, and used leaves as control. A total of 2016; Song et al. 2020), and TaGAMyb, TaFUSCA3 and 103,516, 108,846, 101,620, 105,578 and 104,407 high PBF-D regulate different HMW-GS gene (Glu-1Dy, Glu- confident accessible chromatin regions (ACRs) were 1Bx7 and Glu-1) expression through binding specific identified in leaf and grain in DAP5, DAP9, DAP15 and DNA motifs (Guo et al. 2015; Sun et al. 2017; Zhu et al. DAP20, respectively (Fig. 1A). All the ACRs from the five 2018). Besides, many TFs, such as TaNAC019 control tissues are mostly enriched around the TSS regions glutenin and starch accumulation by targeting both (Supplementary Fig. 1A) and showed high correlation HMW-GS genes and the starch synthase gene TaSSIIa to coefficient ([ 0.80) between two biological replicates synergetic regulation of wheat yield and quality (Gao (Supplementary Fig. 1B–D). Similar with the ACR dis- et al. 2021; Liu et al. 2020). tribution in other large genome species, a significant Wheat grain development can be typically defined faction of ACRs distributed 10 Kb away from genes into three stages: the pregrain-filling phase for 10 days apart from most of ACRs located near the gene (Fig. 1B). after pollination, grain-filling phase for 10–20 days after We then grouped the ACRs as genic ACRs (genic ACRs, pollination and desiccation phase then follows beyond overlapping with genebody), proximal ACRs and distal 30 days after pollination (Wan et al. 2008). The transi- ACRs according the distance between ACR centers and tion of different grain development phase is accompa- their nearest genes. For definition of proximal ACRs and nied by dramatic transcriptional and physiological distal ACRs, we employed different distance cutoffs (2 changes. However, little is known about their regulatory Kb, 4 Kb, 5 Kb, 6 Kb, 8 Kb, 10 Kb), and those within each mechanisms and contribution of each subgenome in cutoff away from genes as proximal ACRs, and those larger than cutoff away from genes as distal ACRs allohexaploid wheat. Accessible chromatin regions are putative cis-regulatory elements (CREs) that are tar- (Supplementary Fig. 2). Interestingly, the proportion of geted by TFs and play vital roles in transcriptional distal ACRs in grain samples were all much higher than regulation. It is reported that the CRE variants are in leaf except in DAP15 with 6 Kb as the cutoff (Sup- associated with many important agronomic traits plementary Fig. 2). When using 2 Kb to distinguish (Adamski et al. 2021; Deplancke et al. 2016; Kouzarides them, more than 50% of ACRs were predominately 2007; Liu et al. 2021; Lu et al. 2017, 2019; Rodgers- enriched in the distal regions of genes in DAP20 Melnick et al. 2016; Tian et al. 2021; Wang et al. 2021b). (Fig. 1C). These results suggested that during grain In this study, we unraveled the regulatory network development, more distal regions were predominantly underlying wheat grain development by combing ATAC- opened and these may be critical for the expression of seq (Assay for Transposase Accessible Chromatin genes regulating starch and storage protein sequencing) and RNA-seq with samples from a series biosynthesis. grain developmental stages. Our studies have not only We further examined the sequence variations around revealed the chromatin accessibility landscapes of ACRs using the reported wheat whole-genome genetic wheat grain development, but also provided valuable variation map (VMap 1.0) (Zhou et al. 2020). To avoid insights into the dynamic regulatory network mediated the influence of biased signals around genes, we only by important transcription factors and the divergences picked up ACR whose centers are[ 2 Kb away from of the three subgenomes during grain development in genes. The results showed that lower sequence varia- the allohexaploid wheat. tion around ACR centers (Fig. 1D), which suggested those ACRs are functionally important. Interestingly, we found that variations were higher in DAP20 than the earlier stages, indicating a lower selection pressure of those genes. The Author(s) 2023 aBIOTECH AC B Genic Proximal Distal 0.8 0.6 0.4 0.2 Leaf DAP5 01 12 34 5 6 78 90 10⁶ 10⁸ DAP9 DAP15 Distance between genes and peak centers (Kbp) DAP20 E F Genic Proximal Distal 0.010 ACR present 1 *** *** ACR absent 0.8 0.009 0.6 *** *** *** 0.008 0.4 0.2 0.007 0.006 0 -2.0Kbp TSS 2.0Kbp Leaf DAP5 DAP9 DAP15 DAP20 GI diffACRs DEGs H Scaled TIS Scaled TPM Cluster organic acid metabolic process oxoacid metabolic process −1.5 0 1.5 −1.5 0 1.5 carboxylic acid metabolic process dACRs in Prmoter DEGs carbohydrate metabolic process Leafs cellular macromolecule biosynthetic process cellular amide metabolic process p < 2.2e-16 peptide metabolic process Leafs amide biosynthetic process DAP5s peptide biosynthetic process gene_expression p < 3.2e-5 glucan biosynthetic process organonitrogen compound biosynthetic process DAP5s carbohydrate biosynthetic process DAP9s oligosaccharide metabolic process disaccharide metabolic process p < 4.0e-4 DAP9s sucrose metabolic process protein modification response to light stimulus p-value DAP5&9s cellular response to radiation 525 3e−04 DAP5&9s cellular response tolight stimulus p < 9.8e-8 cellular response to abiotic stimulus seed development 2e−04 12753 embryo development DAP15s DAP15s protein phosphorylation 1e−04 p < 0.01 phosphorus metabolic process regulation of transcription regulation of RNA metabolic process 13809 regulation of primary metabolic process DAP20s DAP20s regulation of nucleic acid templated transcription GeneRatio regulation of nitrogen compound metabolic process p < 3.5e-10 1211 0.1 phosphorylation DAP15&20s response to water 0.2 27650 response to oxygen containing compound 0.3 DAP15&20s 4285 response to inorganic substance response to acid chemical p < 2.2e-16 13326 response to abiotic stimulus AllGs AllGs p < 0.3 Fig. 1 Chromatin accessibilities are relative to the differential gene expression during wheat grain development. A UpSet Plot showing the ACR distribution across five tissues (DAP, days after pollination). B Density plot showing the distribution between ACR centers and their nearest genes. Black dot line indicates 10 Kb away from genes. C Distribution of ACRs on different genomic features. Genic, ACRs with their centers in gene bodies; Proximal, ACRs with their centers \ = 2 Kb upstream of TSS or \ = 2 Kb downstream of TES; Distal, ACR with their centers [ 2 Kb away from genes. D Distribution of sequence variations around ACR centers. SNPs from Zhou et al. were used. E The expression of genes with ACRs in their 2 Kb nearby regions (ACR present) and without ACRs in their 2 Kb nearby regions (ACR absent) across the five tissues. F Distribution of eight differential ACR (diffACR) clusters on different genomic features. The definition of each group is same with C. G The overlapping information of eight diffACR clusters related genes and eight DEG clusters. H Chromatin accessibilities of regions from 2 Kb upstream to 100 bp downstream of TSS and associated genes’ expression levels. Chromatin accessibilities were represented by Tn5 transposome integration sites (TISs) and expression levels by transcripts per million (TPM). I Gene ontology (GO) enrichment analysis for the overlapped gene groups in G The Author(s) 2023 Leaf DAP5 DAP9 DAP15 DAP20 DAP5 DAP9 DAP15 DAP20 Endosperm Whole seed Leaf DAP5 DAP9 DAP15 DAP20 input Leaf DAP5 DAP9 DAP15 DAP20 Leafs DAP5s Leaf DAP5 DAP9s DAP9 DAP5&9s DAP15 DAP20 DAP15s Leaf DAP20s DAP5 DAP9 DAP15&20s DAP15 AllGs DAP20 Leafs DAP5s DAP9s DAP5&9s DAP15s DAP20s DAP15&20s Relative SNP density No. of Intersections log2(TPM+1) Density 0.0 0.1 0.2 0.3 0.4 0.5 Percentage overlapping genome feature Percentage overlapping genome feature aBIOTECH The chromatin accessibility changes are usually Subgenome-divergent regulation during wheat associated with differential transcriptomic expressions grain development (Chen et al. 2019; Pajoro et al. 2014). To study the relationship between chromatin accessibilities and gene Transcription factors (TFs) play important roles in the regulation of gene expression, while ACRs are predom- expressions, we classified the total wheat genes based on the presences of ACRs within 2 Kb away from genes. inantly TF-binding sites (Liu and Bergmann 2021; Marand et al. 2021). To identify the key TFs responsible The results showed that the expression of genes asso- ciated with ACRs were significantly higher than those for the gene expression dynamics during wheat grain development, we analyzed the enrichments of the 40 without ACRs in all five tissues (Fig. 1E). To identify the putative CREs responsible for the transcriptomic JASPAR cluster motifs (Castro-Mondragon et al. 2022)in expression changes, we identified the differential ACRs the eight diffACRs clusters. The results showed that ARF, (diffACRs) by pair-wise comparisons and then grouped BPC and ERF were enriched in leaf-specific ACRs; ABI, the diffACRs into eight specific clusters, including leaf- ATHB, NAC, MYB, SPL, STZ and bZIP were enriched in specific diffACRs (Leafs), DAP5 specific (DAP5s), DAP9 DAP5s, DAP9s and DAP5&9s; AGL, ARF, E2FA, and REF specific (DAP9s), DAP5 and DAP9 specific (DAP5&9s), were enriched in DAP15s, DAP20s and DAP15&20; Dof was enriched in all grain specific ACRs (Fig. 2A), indi- DAP15 specific (DAP15s), DAP20 specific (DAP20s), DAP15 and DAP20 specifics (DAP15&20s) and all grain cating the specific roles of different TF families in dif- ferent grain developmental stages. specific (AllGs). In total, we identified 31,000, 8346, 9754, 13,829,15,974, 17,253, 39,725 and 1287 specific Bread wheat (Triticum aestivum L.) consists of A, B diffACRs in leafs, DAP5s, DAP9s, DAP5&9s, DAP15s, and D subgenomes derived from three diploid species, DAP20s, DAP15&20s and AllGs groups, respectively and exhibits improved grain yield and nutritional con- (Supplementary Fig. 3A and Supplemental Data Set1). tent relative to diploid and tetraploid wheat (Yang et al. The diffACRs were more enriched in the distal regions 2014). Transcriptomic expression or dynamic modifi- than the total ACRs, indicating accessibilities of distal cation asymmetry of wheat chromosomes contributed to wide adaptability of wheat (Li et al. 2019; Ramirez- regions were more volatile during grain development (Fig. 1F). We further divided DEGs into eight clusters Gonzalez et al. 2018; Wang et al. 2021a). To study the dynamic patterns of gene expression and promoter using the similar strategy and identified the overlaps between DEGs and diffACRs associated genes, which accessibilities of homeologs during wheat grain devel- opment, we conducted the gene expression bias and revealed higher overlapping rates compared to the random shuffled gene sets in all groups except AllGs promoter accessibility bias analysis by ternary plotting (Fig. 2B). 12,218, 11,802, 11,488, 11,319 and 11,597 of (Fig. 1G, Supplementary Fig. 3B, and Supplemental Data Set2). The chromatin accessibilities of promoter regions triad genes showed balanced expression, and more than (– 2 Kb upstream to 100 bp downstream of TSS sites) 75%, 79%, 78%, 76%, 81% and 74% triad genes and the expression of each overlapped gene cluster showed balanced promoter accessibilities in leaf, DAP5, showed similar trends (Fig. 1H), which further con- DAP9, DAP15 and DAP20, respectively (Supplementary firmed that the chromatin accessibility changes are Fig. 4A and B). It should be noted that more unbalanced associated with differential transcriptomic expressions genes were identified in the growing seeds (DAP5, DAP9, DAP15 and DAP20) than in leaves (Supplemental during wheat grain development. To clarify the biological processes mediated by each Data Set3). Among the subgenome differentially expressed genes, the counts of B-suppressed genes and gene group, we performed gene ontology (GO) enrich- ment analysis. The results showed that leaf-specific promoters were more than the other groups, indicating genes were highly enriched in phosphorylation and B subgenome genes were in weaker positions during response to light or abiotic stimulus, while the DAP5, these processes. We further found that only 56%–62% DAP9 and DAP5&9 specific genes were highly enriched genes were balanced in expression among subgenomes in carbohydrate, sucrose metabolic process, in addition, in all tissues and few genes were consistent among amide biosynthesis and nitrogen compound metabolic different stages (Fig. 2C). Similar patterns were also process regulation were enriched in the overlapped found with promoter accessibilities, which indicate genes specific in DAP15 and/or DAP20 (Fig. 1I). These highly diverse regulation among subgenomes during grain development (Fig. 2D). In addition, we conducted enriched GO terms were consistent with the grain development processes. a correlation analysis between gene expression and promoter accessibility, which showed the unbalanced genes were positively correlated with differential The Author(s) 2023 bZIP/RAV bZIP/RAV * * cluster19_NAC cluster19_NAC REF/RAV REF/RAV * * cluster20_NAC cluster20_NAC GT-4 GT-4 * * ARR/GATA ARR/GATA KUA/SRM KUA/SRM cluster31_ARF cluster31_ARF PHL/KAN/HHO PHL/KAN/HHO TGA TGA * * PEND PEND HMG HMG * * cluster12_MYB cluster12_MYB Dof Dof * * cluster39 cluster39 FLC/SOC1/SVP FLC/SOC1/SVP * * cluster16_MYB cluster16_MYB Dof Dof * * KUA/SRM KUA/SRM FLC/SOC1/SVP FLC/SOC1/SVP * * PEND PEND GT-1/AGL GT-1/AGL * * * * PHL/KAN/HHO PHL/KAN/HHO RVE RVE * * * * GT-4 GT-4 WRKY WRKY * * IDD/OBP IDD/OBP TRP TRP * * TCX/SOL TCX/SOL TCX/SOL TCX/SOL * * GT-2 GT-2 GT-2 GT-2 * * cluster12_MYB cluster12_MYB IDD/OBP IDD/OBP * * WRKY WRKY cluster16_MYB cluster16_MYB * * * * cluster39 cluster39 ARR/GATA ARR/GATA * * TRP TRP cluster20_NAC cluster20_NAC * * RVE RVE * * TGA TGA ATHB/HAT ATHB/HAT * * HMG HMG STZ/SIZF STZ/SIZF * * * * cluster19_NAC cluster19_NAC SPL SPL * * * * cluster31_ARF cluster31_ARF ABI3/FUS3 ABI3/FUS3 * * * * CMTA3/FAR1 CMTA3/FAR1 PLT PLT * * bHLH/bZIP bHLH/bZIP FY FY * * E2FA E2FA ATHB/HAT ATHB/HAT * * CDC CDC BPC BPC * * cluster8_ERF cluster8_ERF Foxn1 Foxn1 * * * * ARF/ERF ARF/ERF CMTA3/FAR1 CMTA3/FAR1 * * * * cluster2_ERF cluster2_ERF bHLH/bZIP bHLH/bZIP * * * * cluster32_ARF cluster32_ARF cluster32_ARF cluster32_ARF * * LBD LBD STZ/SIZF STZ/SIZF * * BPC BPC * * SPL SPL * * ARF ARF * * ABI3/FUS3 ABI3/FUS3 Foxn1 Foxn1 * * * * cluster1_TCP cluster1_TCP LFY LFY * * cluster6_TCP cluster6_TCP cluster1_TCP cluster1_TCP * * LBD LBD cluster6_TCP cluster6_TCP * * * * CDC CDC PLT PLT * * * * E2FA E2FA REF/RAV REF/RAV * * * * GT-1/AGL GT-1/AGL cluster26_ARF cluster26_ARF * * bZIP/RAV bZIP/RAV cluster8_ERF cluster8_ERF * * ARF/ERF ARF/ERF * * cluster2_ERF cluster2_ERF * * 20 20 40 40 60 60 80 80 10 100 0 20 20 40 40 60 60 8 80 0 1 100 00 20 20 40 40 60 60 80 80 100 100 20 20 40 40 60 60 8 80 0 1 100 00 20 20 40 40 60 60 80 80 100 100 aBIOTECH Fig. 2 Subgenome-divergent A A regulation during wheat grain Leafs Leafs development A Motif density DAP5s DAP5s of 40 clustered motifs from DAP9s DAP9s JASPAR database (Castro- DAP5_9s DAP5_9s Mondragon et al. 2022)in DAP15s DAP15s eight diffACR clusters. DAP20s DAP20s B Ternary plot showing DAP15&20s DAP15&20s expression bias of the syntenic -2.5 -2.5 0 0 2.5 2.5 genes from three subgenomes. Z-score Z-score The balanced category: (Obeserved/expected) (Obeserved/expected) balanced expression of the triad genes. A dominant, B B B CD CD RNA-seq (gene expression) RNA-seq (gene expression) ATAC-seq (TSS region) ATAC-seq (TSS region) dominant and D dominant: Categor Category y higher expression from A, B or Balanced Balanced D than the other two A suppressed A suppressed orthologs. A suppressed, B B suppressed B suppressed suppressed and D suppressed: D suppressed D suppressed lower expression from A, B or A dominant A dominant D than the other two B dominant B dominant D dominant D dominant orthologs. C Sankey diagram showing expression dynamic of the syntenic genes during Leaf Leaf D DA AP P5 5 D DA AP P9 9 D DA AP P15 15 D DA AP P2 20 0 Leaf Leaf D DA AP P5 5 D DA AP P9 9 D DA AP P1 15 5 D DA AP P20 20 grain development. D Sankey E E Leaf Leaf DAP5 DAP5 DAP9 DAP9 DAP1 DAP15 5 DAP2 DAP20 0 diagram showing the B B B B B B B B B B chromatin accessibility 100 100 100 100 100 100 100 100 100 100 bZIP28 bZIP28 80 80 NA NAC019 C019 80 80 80 80 80 80 80 80 dynamics of the promoter NA NAC019 C019 GAMyb GAMyb NA NAC019 C019 60 60 60 60 60 60 60 60 60 60 NA NAC100 C100 10 100 NA NAC0 C0 C0 C019 C019 0 NA NAC019 C019 P PBF PBF regions (from 2 Kb upstream bZIP28 bZIP28 P28 P28 PBF PBF PBF PBF GAM GA G G G GAMyb GAMyb A AM AM AM 40 40 40 40 bZIP28 bZIP28 P28 P P28 P P P 8 8 8 40 40 40 40 40 40 bZIP28 bZIP28 8 8 8 8 bZIP28 bZIP28 IP IP28 P2 P 8 bZIP28 bZIP28 ZIP28 ZI ZIP28 28 8 8 8 PBF PBF NA NAC100 C100 N N N N N N N N N N N N NA NA A A A A AC C C C C C C019 C019 PBF PB PBF PB P PBF PBF F F GAM G G G GAMyb GAMyb AM Myb yb b NA N N NA A A A AC100 C100 N N N N NA NA A A A AC1 C C C C1 C C100 C100 N N N N NA NA A AC C100 C100 to 100 bp downstream of TSS) PBF PBF 20 20 20 20 20 20 20 20 GAMyb GAMyb NA NAC100 C100 RSR1 RSR1 20 20 RSR1 RSR1 GAMyb GAMyb GAMyb GAMyb RSR1 RSR1 RSR1 RSR1 RSR1 RSR1 RSR1 RSR1 of the syntenic genes during A A D D A A D D A A D D A A D D A A D D grain development. E The FF subgenome diversified expressions of key A A transcription factors. F Motif ** ** ** ** ** ** ** ** ** ** ** ** * * Leafs Leafs B B ** ** ** ** * * * * density of eight diffACR D D ** ** ** ** * * * * A A clusters in each subgenome (* ** ** * * DAP5s DAP5s B B indicate motif density [ 1) ** ** * * D D ** ** * * A A * * * * * * * * * * ** ** ** ** ** ** * * DAP9s DAP9s B B D D A A ** ** * * ** ** * * * * ** ** ** ** * * DAP5&9s DAP5&9s B B D D A A * * * * * * DAP15s DAP15s B B * * * * * * * * * * * * * * * * * * ** ** ** ** ** ** ** ** D D * * * * * * * * * * * * * * * * * * ** ** ** ** ** ** ** ** A A * * * * * * DAP20s DAP20s B B * * * * * * ** ** ** ** * * * * ** ** ** ** * * * * * * * * D D * * * * * * ** ** ** ** * * * * ** ** ** ** * * * * * * * * A A * * * * * * * * * * ** ** * * * * ** ** * * * * * * * * DAP15&20s DAP15&20s B B * * * * ** ** * * * * D D * * * * ** ** * * * * A A * * * * * * * * * * * * * * Alls Alls B B * * * * * * * * * * * * D D * * * * * * * * * * * * -3 -3 0 0 3 3 Z-score Z-score (Obeserved/expected) (Obeserved/expected) promoter accessible regions in all samples (Supple- binding divergence among the three subgenomes, we mentary Fig. 4C). calculated the densities of TF-binding motifs in different The diversity of CREs was considered as key reasons grain specific ACRs among three subgenomes. Although of transcriptomic expression differentiations (Buen- the number of diffACRs in three subgenomes did not rostro et al. 2013; Lu et al. 2017). The previous reported show obvious differences (Supplementary Fig. 4D), a key TFs controlling grain development showed dynamic large proportion of TF-binding sites showed bias subgenome expression differences in leaf, DAP5, DAP9, enrichments among the three subgenomes (Fig. 2F). For DAP15 and DAP20 (Fig. 2E). To identify the specific TF- example, E2FA, CDC, LBD, BPC, ARF and Foxn1 binding The Author(s) 2023 20 20 40 40 60 60 80 80 100 100 20 20 40 40 60 60 80 80 100 100 2 20 0 40 40 60 60 80 80 100 100 20 20 40 40 60 60 80 80 100 100 20 20 4 40 0 60 60 80 80 100 100 aBIOTECH sites showed higher enrichment on A subgenome in development and may contribute to the improved grain leaf-specific ACRs. MYB, NAC and TGA binding sites yield and nutritional content after polyploidization. showed higher enrichment on A subgenome in DAP9s and DAP5&9s diffACRs. bZIP, BPC and KUA TF-binding Starch and gluten regulatory networks mediated sites showed higher enrichment on A subgenome in by key TFs DAP15, DAP20 and DAP15&20s groups. ARF, MYB and TGA binding sites were subgenome divergently enriched During the filling stages, expression of genes involved in in DAP15s, DAP20s and AllGs groups. The densities of the starch and gluten biosynthesis pathways were TF-binding sites were also varied during grain devel- gradually increased. However, the exact cis-regulatory opment. For example, NAC, KUA and MYB binding motifs elements responsible for the tissue-specific expression were enriched in A subgenome at the DAP9s, while bias were still elusive. We picked up the genes related with enriched in B and D subgenomes at DAP15s (Fig. 2E). gluten protein accumulation and starch biosynthesis Gene Ontology (GO) enrichment analysis showed that regulation including HMW-GS, LMW-GS, alpha-gliadin, the balanced and unbalanced genes among subgenomes gamma-gliadin and AGPS1, AGPL1, GBSSI, SBEI, SBEII, during seed development were enriched in different SSII, and SuSy in starch biosynthesis (Supplemental Data pathways, the balanced genes were enriched in DNA Set4). We found that the promoter accessibilities of replication, gene expression and protein folding process, genes involved in starch biosynthesis were higher in while A suppressed genes were mainly enriched in seeds than leaves and increased after DAP5, while those vitamin biosynthetic process, thiamine biosynthetic of glutenin and gliadin genes were increased after DAP9 process and lipid catabolic process, B-suppressed genes (Fig. 3A), which is consistent with the transcriptomic were mainly enriched in carbohydrate homeostasis expression patterns (Fig. 3B). Interestingly, the pro- process, carbohydrate phosphorylation process and moter accessibilities of genes in DAP20 were starting RNA modification process, D suppressed genes were decrease while those genes were still highly expressed mainly enriched in organic substance metabolic process (Fig. 3A and B). Meanwhile, for some genes highly and nitrogen compound biosynthetic process, A domi- expressed after DAP9, higher promoter accessibilities nant genes were mainly enriched in phosphatidylinosi- could be found at DAP5 (Fig. 3A and B), indicating the tol biosynthetic process and carbohydrate catabolic changes of promoter accessibilities are prior to the process, B dominant genes were mainly enriched in transcriptomic expression changes. Taken together, carbohydrate biosynthetic process and metal ion these results indicated that the chromatin accessibilities transport process, D dominant genes were mainly were tightly associated with the expression changes of those grain-filling-related genes. enriched in protein phosphorylation and aminoglycan metabolic process. These enriched GO terms were con- To find the key TFs controlling grain development, we sistent with the grain development processes (Supple- further calculated the densities of TF-binding motifs in mentary Fig. 5). To prove the specific binding of TFs in the promoters of genes related with gluten protein different tissues, we combined the published genome- accumulation and starch biosynthesis, which suggested wide profiling of TFBSs (transcription factor-binding that NAC, bZIP and MYB family TFs play important roles sites) in common wheat (DAP-seq data) (Zhang et al. (Supplementary Fig. 7). Consistently, the key TFs rela- 2022) and our ATAC-seq data of different tissues to ted with grain development, including TaNAC019, Tab- better identify the potential in vivo TFBSs in different ZIP28, and TaGAMyb showed increased promoter tissues. The results showed that more that 60% TFs’ accessibilities and transcriptomic expression patterns in (NAC6A-1, NAC6B-1, NAC-6D-1 and MYB-7A-3) binding grains (Fig. 3C and Supplementary Fig. 8). We then sites were different among various tissues (Supple- constructed a regulatory network based on the chro- mentary Fig. 6A). In addition, we analyzed the sub- matin accessibility data and expression pattern of those genome-divergent regulation by TFs (NAC6A-1, NAC6B- key TFs and genes involved in starch biosynthesis and 1, NAC-6D-1 and MYB-7A-3), and found that nearly 80% protein accumulations. The results showed that the traids showed subgenome-divergent regulations (Sup- three copies of TaSPA were involved in the networks in plementary Fig. 6B). These results indicated tissue- DAP9 and DAP15; TaNAC019 and TaNAC100 in DAP15 specific and subgenome-divergent regulation by multi- and DAP20; while the others mostly in all the stages ple TFs during wheat grain development. The high after DAP9 (Fig. 3D). Subgenome specific regulation variations of TF-binding sites enrichment bias among could be found with FUSCA3, of which the copy from D subgenomes at different stages suggested the three subgenome were involved in networks in DAP9 and subgenomes worked coordinately during wheat grain DAP15, while the other two in all the three stages. Similarly, the copy of TaGAMYB from D subgenome were The Author(s) 2023 aBIOTECH -200 2 36 TaSBEIIa TaSuSy2 TaSBEIb A B C Scaled TIS log10(TPM+1) TraesCS2A02G293400 TraesCS2B02G194200 TraesCS7A02G549300 TaAGPL1-A [0 - 77] [0 - 48] [0 - 99] TaAGPL1-B Leaf TaAGPL1-D [0 - 77] [0 - 99] [0 - 48] DAP5 TaAGPS1a-A TaAGPS1a-B [0 - 48] [0 - 77] [0 - 99] DAP9 TaAGPS1a-D TaAGPS2-D [0 - 77] [0 - 99] [0 - 48] DAP15 TaGBSSI-A TaGBSSI-D [0 - 48] [0 - 77] [0 - 99] DAP20 TaISA1-B TaISA1-D [0 - 249875] [0 - 668244] [0 - 25185] Leaf TaISA2-B [0 - 25185] [0 - 249875] [0 - 668244] TaPUL-A DAP5 TaPUL-B TaPUL-D [0 - 249875] [0 - 668244] [0 - 25185] DAP9 TaSBEIa-A TaSBEIa-D [0 - 249875] [0 - 668244] [0 - 25185] DAP15 TaSBEIb-A TaSBEIb-B [0 - 249875] [0 - 25185] [0 - 668244] DAP20 TaSBEIc-D TaSBEIIa-A TaSBEIIa-B TaGBSSI TabZIP28 TaNAC019 TaSBEIIa-D TraesCS7A02G070100 TraesCS2D02G146100 TraesCS3B02G092800 TaSBEIIb-A TaSBEIIb-B [0 - 28] [0 - 125] TaSBEIIb-D [0 - 149] Leaf TaSSIc-A [0 - 125] [0 - 28] [0 - 149] TaSSIc-D DAP5 TaSSIIa-A [0 - 125] [0 - 149] [0 - 28] TaSSIIa-B DAP9 TaSSIIa-D [0 - 125] [0 - 149] [0 - 28] TaSSIIIa-D DAP15 TaSSIIIb-A TaSSIIIb-B [0 - 125] [0 - 149] [0 - 28] DAP20 TaSuSy1-A [0 - 25580] [0 - 383539] [0 - 53686] TaSuSy1-B Leaf TaSuSy1-D [0 - 25580] TaSuSy2-A [0 - 383539] [0 - 53686] DAP5 TaSuSy2-B [0 - 25580] TaSuSy2-D [0 - 383539] [0 - 53686] DAP9 TaSuSy4-A TraesCS7D02G064300TraesCS7D02G117800 TaSuSy4-D [0 - 25580] [0 - 383539] [0 - 53686] DAP15 TaSuSy5-A TaSuSy5-B [0 - 25580] DAP20 [0 - 383539] [0 - 53686] TaSuSy5-D Omega-gliadin TaNAC100 D TaNAC019 Interaction in B DAP9 B DAP15 DAP20 D Gamma-gliadin FUSCA3 SPR Alpha-gliadin B HMW-glutenin TaGAMYB TabZIP28 LMW-glutenin E15 E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E1 E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E1 1 1 15 15 1 1 1 15 15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 15 15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 E9 E E E5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 9 9 9 9 9 9 9 9 9 9 Glutenin Gliadin Starch synthesis A PBF-D TaSPA Fig. 3 Construction of regulatory networks for starch and gluten biosynthesis by integration the chromatin accessibility and expression data. A and B Promoter chromatin accessibility dynamics (A) and expression dynamics (B) of genes involved in glutenin and starch accumulation in leaf and grain samples. Regions from 2 Kb upstream to 100 bp downstream of transcriptional start sites were considered as promoters. C Genome browser tracks showing the expression and promoter chromatin accessibility patterns of starch synthesis genes (TaSBEIIa, TaSuSy2, TaSBEIb and TaGBSSI) and key TF genes (TaZIP28 and TaNAC019). D Regulatory networks involving the starch and gluten biosynthesis genes and key TFs. TF families are shown as circles with different colors. Starch biosynthesis genes are shown as grey square; Gliadin genes as orange triangle; Glutenin genes as blue diamond. TF-to-gene edges are shown with three colors represent three grain phases. DAP9, green; DAP15 blue; DAP20, pink in DAP9 and A subgenome copy in DAP9 and DAP15 and protein synthesis. To further verify the predicted while B subgenome copy in all the three stages. The transcriptional regulation network, we conducted the results also revealed most of these key TFs were dual-luciferase transcriptional activity assay (Fig. 4A–B involved in both starch synthesis and protein accumu- and Supplemental Data Set5). The assay showed that lations, indicating a complex network between starch the reporter activities driven by the promoters of The Author(s) 2023 Leaf DAP5 DAP9 DAP15 DAP20 Leaf DAP5 DAP9 DAP15 DAP20 RNA-seq ATAC-seq RNA-seq ATAC-seq aBIOTECH TaHWM-1Dpro TaHWM TaAGPL1-1Bpro TaAGPL1 TaISA2-1Bpro TaISA2 TF motif TF motif TF motif 2 2 1 1 cluster_19_NAC cluster_19_NAC cluster_19_NAC 1 0 0 123456789 1011121314 1 2 3 4 5 6 7 8 9 1011121314 123456789 1011121314 2 2 1 1 cluster_20_NAC cluster_20_NAC cluster_20_NAC 0 0 112 24 112 24 112 24 2 2 Position 1 1 cluster_12_MYB cluster_12_MYB 0 0 111 22 111 22 2 2 1 1 cluster_16_MYB cluster_16_MYB 018 16 018 16 Position Position BC TaAGPL1-1Bpro Effector 0.04 *** 0.004 TaHWM-1Dpro 35S GFP NOS *** 0.03 0.003 0.02 35S MYB_B-GFP NOS 0.002 0.01 TaISA2-1Bpro 35S NAC100_A-GFP NOS TaHWM-1Dpro 0.001 ** *** 0.0008 0.003 35S NAC100_D-GFP NOS 0.0006 0.002 Reporter 0.0004 0.001 0.0002 TaHWM-1Dpro LUC 35S REN NOS 0.0000 0.000 TaAGPL1-1Bpro LUC 35S REN NOS TaISA2-1Bpro LUC 35S REN NOS TaISA2-1Bpro TaISA2-1Bpro TaAGPL1-1Bpro TaAGPL1-1Bpro 0.0015 0.0020 0.015 0.015 ** **** 0.0015 **** *** 0.0010 0.010 0.010 *** ** *** 0.0010 * **** *** 0.0005 0.005 0.005 **** 0.0005 0.000 0.0000 0.0000 0.000 Fig. 4 Dual-luciferase transcriptional activity assay to assess the capability of key TFs to transactivate predicted target gene expression. A The distribution of key TF-binding motifs in the promoters of starch and protein synthesis genes. B Schematic diagrams of the effector and reporter constructs for the dual-luciferase transcriptional activity assay. C Reporter assay showed that NAC100-A or NAC100- D regulates both starch and protein synthesis pathways. D NAC100 and MYB-B coordinately regulate TaISA2 and TaAGPL1 genes. LUC/ REN indicates the signal ratio of LUC (firefly luciferase) to REN (Renilla reniformis luciferase) activity. Data are represented as means ± SD from three replicates. (Student’s t-test; *P \ 0.05, **P \ 0.01) TaAGPL1-1B, TaHWM-1D and TaISA2-1B were markedly (TaAGPL1-1B and TaISA2-1B) (Fig. 4D). In summary, we activated in the presence of NAC100-A or NAC100-D, have explored the chromatin accessibility landscapes suggested that NAC100 regulates both starch and pro- and identified the key TFs and regulatory network that tein synthesis pathways (Fig. 4C). Moreover, NAC100 were involved in the starch and protein biosynthesis and MYB could coordinately regulate same targets during wheat grain development. The Author(s) 2023 GFP+GFP GFP+NAC100_A GFP+MYB_B GFP NAC100_A+MYB_B NAC100_A GFP NAC100_A GFP+GFP GFP+NAC100_D GFP+MYB_B NAC100_D+MYB_B GFP NAC100_D GFP GFP+GFP NAC100_D GFP+NAC100_A GFP+MYB_B NAC100_A+MYB_B GFP+GFP GFP+NAC100_D GFP+MYB_B NAC100_D+MYB_B LUC/REN Bits Bits LUC/REN Bits Bits Bits Bits LUC/REN LUC/REN LUC/REN LUC/REN Bits Bits Bits Bits aBIOTECH DISCUSSION METHODS Wheat grain development is critical for the quality and Plant material and growth conditions yield of wheat, a series of genes regulated storage protein and starch biosynthesis have been discovered in Wheat (Triticum aestivum; BBAADD, 2n = 6x = 42) the previous study (Xiao et al. 2022). Transcriptional cultivar AK58 was transferred to the greenhouse with regulation of gluten and starch genes involves complex conditions day/night cycle set at 16 h (22 C)/8 h network between cis- and trans-acting factors, but little (19 C) after vernalization for 1 mouth in cold chamber is known about this regulation. The accessible chro- (4 C). Wheat endosperm was excised by removing the matin regions often associated with TFs related CREs seed embryo at 5 and 9 days after fertilization. Wheat (Lu et al. 2017; Marand et al. 2021; Pei et al. 2022). seed was collected at 15 and 20 days after fertilization. Here, we constructed a chromatin accessibility atlas for leaf and four different stages of seed during grain RNA-seq development. This regulatory atlas enabled specific ACR and TF identification at different phase of wheat grain The collected wheat endosperm and seed were flash- development. We found NAC, bZIP and MYB binding frozen with liquid N2 immediately. Total RNA was TM motifs are highly enriched in the grain specific ACRs as extracted with TRIzol Reagent (Invitrogen; previous reported (Gao et al. 2021; Guo et al. 2015; Liu 15,596–026) following the manufacturer’s instructions. et al. 2020; Song et al. 2020). More importantly, some RNA-Seq libraries were constructed and sequenced by novel TFs, including REF, Dof and HMG were also dis- Berry Genomics (Beijing, China), the libraries were covered highly enriched in the these specific ACRs. sequenced on the Illumina NovaSeq 6000 platform and Further studies on TFs in those families will help to produced 150-bp paired-end reads. understand the molecular basis of grain development. Chromatin accessibilities’ changes are associated ATAC-seq with genes’ differential expression. In this study, we combined ATAC-seq with RNA-seq and constructed the ATAC-seq was performed as described previously (Lu dynamic regulatory network during wheat grain devel- et al. 2017). For each sample, approximately 1 g flash- opment, and further validated the regulatory network frozen wheat samples were chopped with a razor blade by dual-luciferase transcriptional activity assays. We in 1 mL ice-prechilled lysis buffer (15 mM Tris–HCl pH found that the phase-specific ACRs are associated with 7.5, 20 mM NaCl, 80 mM KCl, 0.5 mM spermine, 5 mM phase-specific DEGs, and the function of phase-specific 2-Mercaptoethanol, 0.2% TritonX-100). The chopped genes is related to specific the biological process of slurry containing crude nuclei extract was filtered twice specific tissue. The seed-specific genes were highly through a 40 lm filter. The crude nuclei were stained enriched in carbohydrate and sucrose metabolic pro- with DAPI (sigma, catalog number: D9542) and loaded cess, and seed development processes. We also found to a flow cytometer (BD FACSCanto) for selected. Nuclei the changes of promoter accessibilities of starch and pellets were obtained after centrifuged and washed with protein biosynthesis genes were prior to their expres- Tris-Mg buffer (10 mM Tris–HCl pH 8.0, 5 mM MgCl ). sion changes, indicating those RNAs were likely accu- The obtained nuclei were incubated with 3.5 lLTn5 mulated rather than rapidly synthesized. Besides, we transposomes in 40 lL TTBL buffer (TruePrep DNA found the homeolog dynamic transcriptomic expression Library Prep Kit V2 for Illumina, Vazyme Biotech co., ltd, divergence during wheat grain development (Xiang TD501) for each sample at 37 C for 30 min without et al. 2019) are associated with the divergence of pro- rotation. We purified the integration products with a TM moter accessibilities and the TF-binding sites in the NEB Monarch DNA Cleanup Kit (T1030S) and then three subgenomes at different phases of wheat grain amplified for 10–13 cycles using the NEBNext Ultra II development. Therefore, those accessible regions in Q5 master mix (M0544L). PCR cycles were determined different phases could help to accurately discover the as described previously (Lu et al. 2017). Amplified key cis- and trans-acting factors during wheat grain libraries were purified with Hieff NGS DNA Selection development. In summary, our data provided a tran- Beads (Yeasen, 12601ES03) to remove free index scriptional regulation resource to understand the com- primers. plex networks and identified the key regulatory elements and TFs underlying wheat grain development. The Author(s) 2023 aBIOTECH Supplementary InformationThe online version contains Processing of RNA-seq data supplementary material available at https://doi.org/10.1007/ s42994-023-00095-8. Raw reads were preprocessed by fastp for filtering low- quality reads and adapter trimming. The remaining Acknowledgements We would like to acknowledge Ting Li from the flow cytometry core of Institute of Genetics and Develop- reads were aligned to the wheat reference genome using mental Biology, Chinese Academy of Sciences for supporting the hisat2 with default parameters. All aligned reads were nuclear sorting. This project was financially supported by the sorted by SAMtools v1.3.1 (Swiezewski et al. 2009). To Outstanding Young Scientist Foundation of NSFC (Overseas), the compare the expression from different phases, TPM was Central Public-interest Scientific Institution Basic Research Found (S2022ZD02), the Fundamental Research Funds from Institute of calculated by TPMCalculator. Differential expressed Crop Sciences, Chinese Academy of Agricultural Sciences genes (DEGs) were defined as adjusted P-value \ 0.01 (S2020YC07 and S2021YC03) and CAAS Agricultural Science and and log2(FoldChange) [ 1, which determined by Technology Innovation Program, China (CAAS-ZDRW202002). DESeq2. Author contributions ZL conceived and designed the project. HP and YSL performed ATAC-seq experiments. YHL, PL, and XR were Processing of ATAC-seq data responsible for wheat materials. HP, YSL and JLZ analyzed the data. ZL, YSL, YHL, and PL contributed to project discussion. HP Raw reads were trimmed with fastp with default and ZL wrote the manuscript draft, and YSL, YHL and PL revised it. parameters. Trimmed reads were aligned to the wheat reference genome using Bowtie v2.2.4 with the follow- Data availability All the raw sequencing data generated during the current study are available in the Gene Expression Omnibus ing parameters: ‘‘bowtie2-X 1000 –very-sensitive’’. (GEO) database (https://www.ncbi.nlm.nih.gov/geo) under Aligned reads were sorted using SAMtools v1.3.1 accession number GSE214739. (Swiezewski et al. 2009) and clonal duplicates were removed using Picard version v2.16.0 http://broad Declarations institute.github.io/picard/). Peak calling was performed Conflict of interest The authors declare that they have no con- as described previously. A black list was first generated flict of interest. using the control samples as input with MACS2. The raw peaks were first called with MACS2 with the following Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, parameter ‘‘–keep-dup all –nomodel –extsizes 150 –shift adaptation, distribution and reproduction in any medium or for- -75’’ (Zhang et al. 2008); then split into 150 bp bins mat, as long as you give appropriate credit to the original with 50 bp overlapping; bins passed filtering by a Tn5 author(s) and the source, provide a link to the Creative Commons integration site density cutoff were selected and merged licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s with bedtools with ‘‘– d 150’’; peaks overlapped with the Creative Commons licence, unless indicated otherwise in a credit blacklist peaks, or homologous to plant organelles DNA line to the material. If material is not included in the article’s (NCBI) were discard and the rest were considered as Creative Commons licence and your intended use is not permitted high-quality ACRs. Bigwig files were used to visualize by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To the peaks in the integrative genomics viewer. view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. Dual-luciferase reporter gene assay Dual-luciferase reporter gene assay was performed in N. benthamiana leaves by Agrobacterium tumefaciens-me- References diated transient expression system. The promoters of TaAGPL1-1B, TaHWM-1D and TaISA2-1B were amplified Adamski NM, Simmonds J, Brinton JF, Backhaus AE, Chen Y, Smedley M, Hayta S, Florio T, Crane P, Scott P et al (2021) from AK58 genomic DNA and inserted into the pGreenII Ectopic expression of Triticum polonicum VRT-A2 underlies 0800-LUC vector as reporter plasmids (Hellens et al. elongated glumes and grains in hexaploid wheat in a dosage- 2005). The CDS of NAC100-A, NAC100-D and MYB-B dependent manner. Plant Cell 33:2296–2319 were cloned in frame into the pEG1300-GFP vector as Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ (2013) Transposition of native chromatin for fast and sensitive effector plasmids. N. benthamiana plants at the six-leaf epigenomic profiling of open chromatin, DNA-binding pro- stage were co-infiltrated with A. tumefaciens strain teins and nucleosome position. Nat Methods 10:1213–1218 GV3101 harboring different combinations of these Castro-Mondragon JA, Riudavets-Puig R, Rauluseviciute I, Lemma plasmids. The dual-luciferase reporter assay system RB, Turchi L, Blanc-Mathieu R, Lucas J, Boddie P, Khan A, Manosalva Perez N et al (2022) JASPAR 2022: the 9th release (Yeasen, China) was used to measure LUC and REN of the open-access database of transcription factor binding activities and three independent replications were profiles. Nucl Acids Res 50:D165–D173 conducted. The Author(s) 2023 aBIOTECH Chen S, Lake BB, Zhang K (2019) High-throughput sequencing of Lu Z, Marand AP, Ricci WA, Ethridge CL, Zhang X, Schmitz RJ the transcriptome and chromatin accessibility in the same (2019) The prevalence, evolution and chromatin signatures cell. Nat Biotechnol 37:1452–1457 of plant regulatory elements. Nat Plants 5:1250–1259 Deplancke B, Alpern D, Gardeux V (2016) The genetics of Marand AP, Chen ZL, Gallavotti A, Schmitz RJ (2021) A cis- transcription factor DNA binding variation. Cell 166:538–554 regulatory atlas in maize at single-cell resolution. Cell Gao Y, An K, Guo W, Chen Y, Zhang R, Zhang X, Chang S, Rossi V, Jin 184:3041 F, Cao X et al (2021) The endosperm-specific transcription Pajoro A, Madrigal P, Muino JM, Matus JT, Jin J, Mecchia MA, factor TaNAC019 regulates glutenin and starch accumulation Debernardi JM, Palatnik JF, Balazadeh S, Arif M et al (2014) and its elite allele improves wheat grain quality. Plant Cell Dynamics of chromatin accessibility and gene regulation by 33:603–622 MADS-domain transcription factors in flower development. Gu A, Hao P, Lv D, Zhen S, Bian Y, Ma C, Xu Y, Zhang W, Yan Y Genome Biol 15:R41 (2015) Integrated proteome analysis of the wheat embryo Pei H, Teng W, Gao L, Gao H, Ren X, Liu Y, Jia J, Tong Y, Wang Y, Lu Z and endosperm reveals central metabolic changes involved in (2022) Low-affinity SPL binding sites contribute to sub- the water deficit response during grain development. J Agric genome expression divergence in allohexaploid wheat. Sci Food Chem 63:8478–8487 China Life Sci 1–16 Guo W, Yang H, Liu Y, Gao Y, Ni Z, Peng H, Xin M, Hu Z, Sun Q, Yao Y Ramirez-Gonzalez RH, Borrill P, Lang D, Harrington SA, Brinton J, (2015) The wheat transcription factor TaGAMyb recruits Venturini L, Davey M, Jacobs J, van Ex F, Pasha A et al (2018) histone acetyltransferase and activates the expression of a The transcriptional landscape of polyploid wheat. Science high-molecular-weight glutenin subunit gene. Plant J 361:6403 84:347–359 Ren J, Jiang Z, Li W, Kang X, Bai S, Yang L, Li S, Zhang D (2022) Hellens RP, Allan AC, Friel EN, Bolitho K, Grafton K, Templeton MD, Characterization of glutenin genes in bread wheat by third- Karunairetnam S, Gleave AP, Laing WA (2005) Transient generation RNA sequencing and the development of a Glu- expression vectors for functional genomics, quantification of 1Dx5 marker specific for the extra cysteine residue. J Agric promoter activity and RNA silencing in plants. Plant Methods Food Chem 70:7211–7219 1:13 Rodgers-Melnick E, Vera DL, Bass HW, Buckler ES (2016) Open Huang L, Tan H, Zhang C, Li Q, Liu Q (2021) Starch biosynthesis in chromatin reveals the functional maize genome. Proc Natl cereal endosperms: An updated review over the last decade. Acad Sci USA 113:E3177-3184 Plant Commun 2:100237 She KC, Kusano H, Koizumi K, Yamakawa H, Hakata M, Imamura T, Kang G, Liu G, Peng X, Wei L, Wang C, Zhu Y, Ma Y, Jiang Y, Guo T Fukuda M, Naito N, Tsurumaki Y, Yaeshima M et al (2010) A (2013) Increasing the starch content and grain weight of novel factor FLOURY ENDOSPERM2 is involved in regulation common wheat by overexpression of the cytosolic AGPase of rice grain size and starch quality. Plant Cell 22:3280–3294 large subunit gene. Plant Physiol Biochem 73:93–98 Song Y, Luo G, Shen L, Yu K, Yang W, Li X, Sun J, Zhan K, Cui D, Liu D Kouzarides T (2007) Chromatin modifications and their function. et al (2020) TubZIP28, a novel bZIP family transcription Cell 128:693–705 factor from Triticum urartu, and TabZIP28, its homologue Li Z, Wang M, Lin K, Xie Y, Guo J, Ye L, Zhuang Y, Teng W, Ran X, from Triticum aestivum, enhance starch synthesis in wheat. Tong Y et al (2019) The bread wheat epigenomic map reveals New Phytol 226:1384–1398 distinct chromatin architectural and evolutionary features of Sun F, Liu X, Wei Q, Liu J, Yang T, Jia L, Wang Y, Yang G, He G (2017) functional genetic elements. Genome Biol 20:139 Functional characterization of TaFUSCA3, a B3-superfamily Liu A, Bergmann DC (2021) How to build a crop plant: defining transcription factor gene in the wheat. Front Plant Sci 8:1133 the cis-regulatory landscape of maize. Cell 184:2804–2806 Swiezewski S, Liu F, Magusin A, Dean C (2009) Cold-induced Liu G, Wu Y, Xu M, Gao T, Wang P, Wang L, Guo T, Kang G (2016) silencing by long antisense transcripts of an Arabidopsis Virus-induced gene silencing identifies an important role of Polycomb target. Nature 462:799–802 the TaRSR1 transcription factor in starch synthesis in bread Tian H, Li Y, Wang C, Xu X, Zhang Y, Zeb Q, Zicola J, Fu Y, Turck F, Li wheat. Int J Mol Sci 17(10):1557 L et al (2021) Photoperiod-responsive changes in chromatin Liu J, Chen Z, Wang Z, Zhang Z, Xie X, Wang Z, Chai L, Song L, Cheng accessibility in phloem companion and epidermis cells of X, Feng M et al (2021) Ectopic expression of VRT-A2 underlies Arabidopsis leaves. Plant Cell 33:475–491 the origin of Triticum polonicum and Triticum petropavlovskyi Wan Y, Poole RL, Huttly AK, Toscano-Underwood C, Feeney K, with long outer glumes and grains. Mol Plant Welham S, Gooding MJ, Mills C, Edwards KJ, Shewry PR et al 14(9):1472–1488 (2008) Transcriptome analysis of grain development in Liu L, He ZH, Yan J, Zhang Y, Xia XC, Pena RJ (2005) Allelic hexaploid wheat. BMC Genom 9:121 variation at the Glu-1 and Glu-3 loci, presence of the 1B.1R Wang M, Li Z, Zhang Y, Zhang Y, Xie Y, Ye L, Zhuang Y, Lin K, Zhao F, translocation, and their effects on mixographic properties in Guo J et al (2021a) An atlas of wheat epigenetic regulatory Chinese bread wheats. Euphytica 142:197–204 elements reveals subgenome divergence in the regulation of Liu Y, Han C, Deng X, Liu D, Liu N, Yan Y (2018) Integrated development and stress responses. Plant Cell 33:865–881 physiology and proteome analysis of embryo and endosperm Wang X, Aguirre L, Rodriguez-Leal D, Hendelman A, Benoit M, highlights complex metabolic networks involved in seed Lippman ZB (2021b) Dissecting cis-regulatory control of germination in wheat (Triticum aestivum L.). J Plant Physiol quantitative trait variation in a plant stem cell circuit. Nat Plants 7:419–427 229:63–76 Liu Y, Hou J, Wang X, Li T, Majeed U, Hao C, Zhang X (2020) The Wang Y, Hou J, Liu H, Li T, Wang K, Hao C, Liu H, Zhang X (2019) NAC transcription factor NAC019-A1 is a negative regulator TaBT1, affecting starch synthesis and thousand kernel weight, of starch synthesis in wheat developing endosperm. J Exp Bot underwent strong selection during wheat improvement. J Exp 71:5794–5807 Bot 70:1497–1511 Lu Z, Hofmeister BT, Vollmers C, DuBois RM, Schmitz RJ (2017) Xiang DQ, Quilichini TD, Liu ZY, Gao P, Pan YL, Li Q, Nilsen KT, Combining ATAC-seq with nuclei sorting for discovery of cis- Venglat P, Esteban E, Pasha A et al (2019) The Transcriptional regulatory regions in plant genomes. Nucl Acids Res 45:e41 landscape of polyploid wheats and their diploid ancestors The Author(s) 2023 aBIOTECH during embryogenesis and grain development. Plant Cell subgenome-convergent and -divergent transcription in com- 31:2888–2911 mon wheat. Nat Commun 13:6940 Xiao J, Liu B, Yao Y, Guo Z, Jia H, Kong L, Zhang A, Ma W, Ni Z, Xu S Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al (2022) Wheat genomic study for genetic improvement of Nusbaum C, Myers RM, Brown M, Li W et al (2008) Model- traits in China. Sci China Life Sci 65:1718–1775 based analysis of ChIP-Seq (MACS). Genome Biol 9:R137 Yamamori M, Fujita S, Hayakawa K, Matsuki J, Yasui T (2000) Zhou Y, Zhao X, Li Y, Xu J, Bi A, Kang L, Xu D, Chen H, Wang Y, Wang Genetic elimination of a starch granule protein, SGP-1, of YG et al (2020) Triticum population sequencing provides wheat generates an altered starch with apparent high insights into wheat adaptation. Nat Genet 52:1412–1422 amylose. Theor Appl Genet 101:21–29 Zhu J, Fang L, Yu J, Zhao Y, Chen F, Xia G (2018) 5-Azacytidine Yang C, Zhao L, Zhang H, Yang Z, Wang H, Wen S, Zhang C, Rustgi S, treatment and TaPBF-D over-expression increases glutenin von Wettstein D, Liu B (2014) Evolution of physiological accumulation within the wheat grain by hypomethylating the responses to salt stress in hexaploid wheat. Proc Natl Acad Glu-1 promoters. Theor Appl Genet 131:735–746 Sci USA 111:11882–11887 Zhang Y, Li Z, Liu J, Zhang Y, Ye L, Peng Y, Wang H, Diao H, Ma Y, Wang M et al (2022) Transposable elements orchestrate The Author(s) 2023

Journal

aBIOTECHSpringer Journals

Published: Mar 1, 2023

Keywords: Wheat; Chromatin accessibility; Subgenome-divergence; Regulatory network; Grain development

References