The bread wheat variety Aikang 58 (AK58) was used in this study. Seeds were sterilized in 75% alcohol and then washed three times with water. Seeds were germinated on filter paper at 25°C under 16/8 light/dark conditions. Two days later, the radicles were harvested and the remaining part was transferred to the soil. The plants were grown at 22°C in greenhouse. We harvested young spike at double-ridge stage and ovary at 2 days after pollination. Half of the above materials were stored at -80°C for RNA-seq and the other half for nuclei extraction.

Except the double-ridge stage young spike which had been stored at -80°C for more than 1 month, fresh material was collected and was immediately chopped in any of D-S-0.4 buffer (0.4 M D-Sorbitol, 20 mM MES, 20 mM KCl, 10 mM MgCl₂, 3 mM DTT, 0.4% Triton X-100 1% BSA, 1 tablet/50 ml protease inhibitor), D-Sorbitol buffer (0.6 M D-Sorbitol, 20 mM MES, 20 mM KCl, 10 mM MgCl₂, 1 mM β-mercaptoethanol, 1×PMSF buffer, 0.2% TritonX-100 0.1 mg/ml BSA, 1 tablet/50 ml protease inhibitor), Chopping Buffer (15 mM Tris-HCl, 20 mM NaCl, 80 mM KCl, 0.5 mM spermine, 5 mM 2-ME, 0.4% TritonX-100) or PP Buffer (1×PBS containing 1 tablet/50 ml protease inhibitor, 0.4% TritonX-100). After chopping, the liquid was filtered with 30 μm strainer.

Crude nuclei were stained with DAPI and loaded into a flow cytometer (BD FacsAria II SORP). 50,000 nuclei were sorted into 400 μl buffer. The nuclei were collected by centrifugation at 1,000 g for 10 min at 4°C. After removing the supernatant but 20 μl liquid remained, premixed 1×TD buffer (10 mM Tris-HCl, 5 mM MgCl₂, 10% v/v DMF) and 1.5 μl TTEMix V50 (Vazyme, #TD501) were added to the nuclei solution and incubated at 37°C for 30 min. The reaction was terminated and DNA was purified using the QIAGEN MinElute PCR Purification Kit (QIAGEN, #28006). The eluted DNA was used to prepare library with either TruePrep DNA Library Prep Kit V2 (Vazyme, #TD501) or NEBNext High-Fidelity2X PCR Master Mix (NEB, M0541). The primer sequences are listed in Table S1 . The precise number of PCR cycles was determined by qPCR (Bajic et al., 2018). The libraries were purified with AMPure beads. The concentration of ATAC-seq libraries were measured using Qubit. The fragment sizes of the library were calculated by Agilent 2100 Bioanalyzer. The libraries were sequenced by the NovaSeq 6000 system with the 150-bp paired-end mode.

Tissue was snap-frozen in liquid nitrogen and stored at -80°C. Total RNA from young spike at double-ridge stage and ovary at 2 days after pollination was extracted using Trizol (Ambion, #15596018). The concentration and purity of RNA were detected by Nanodrop (Thermo Fisher Scientific) and electrophoresis. RNA-seq libraries from the two tissues were prepared using the VAHTS^® Universal V8 RNA-seq Library Prep Kit for MGI (Vazyme, #NRM605). RNA-seq libraries were sequenced by DNBSEQ-T7 (BGI, Shenzhen) with the 150-bp paired-end mode. RNA-seq was performed in triplicates.

The quality of raw sequencing data was estimated using Trimmomatic v0.32 (https://github.com/usadellab/Trimmomatic) with parameters of “LEADING: 20 TRAILING: 20 SLIDINGWINDOW: 4: 15 MINLEN: 36”. The clean data was mapped to IWGSCv1.0 using bowtie2 v2.4.4 with parameters “–sensitive -k 3” (Langmead et al., 2009). After removing PCR duplicates using Picard v2.23.9 and low-quality reads (Filter unmatched & MAPQ < 5 reads) using SAMtools v1.9 (Li et al., 2009), the high quality clean reads were analyzed using MACS2 v2.2.7.1 (Zhang et al., 2008) with the following parameters: –nomodel –shift -75 –extsize 150 –to-large -B –SPMR. The shift length is our empirical value after hundreds of predictions using phantompeakqualtools (https://github.com/kundajelab/phantompeakqualtools). In order to remove false positive peaks, we used IDR to get the consistent peaks, which was recommended by ENCODE. The FRiP and the SPOT coefficients were used to evaluate the signal-to-noise ratio of the libraries. The final bam files produced from reads alignment were converted to bigwig format for visualization on genome browser using bamCoverage of deepTools v3.5.0 (Ramirez et al., 2016). The correlation between samples were calculated and profiled by the multiBigwigSummary and plotCorrelation command of deepTools toolkit, respectively. Generally, the Pearson Correlation Coefficient greater than 0.8 is regarded as better correlation. Further, the deepTools computeMatrix command was used to analyze the signal density around the TSS and TES with the following parameters: computeMatrix scale-regions -p 10 -b 3000 -a 3000 -m 5000 –skipZeros. The distance from the peak summit to gene TSS was calculated by BEDTools v2.27, and then the locations of peaks were divided into genebody, promoter, intergenic upstream and downstream region using our developed R script. Genes with at least one specific peak within 4 kb upstream and 1.5 kb downstream of TSS were labeled as chromatin accessibility-related genes. ClusterProfiler was used to perform GO analysis with the obtained OCRs related genes (q-value cutoff = 0.05) (Wu et al., 2021). We developed a pipeline, called wheatATAC to analyze ATAC-seq data in wheat (https://github.com/hcph/wheatATAC.git).

We randomly extracted 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90%, 100% and virtual 110%, 120%, 130%, 140%, 150%, 160%, 170%, 180%, 190% and 200% of the high quality unique mapped clean reads via SAMtools toolkit combined with Linux split command and performed peak calling via MACS2 (Zhang et al., 2008; Li et al., 2009). The saturated sequencing depth was determined via the saturationPlot function of the R-package, ATACseqQC (Ou et al., 2018).

Since grinding in liquid nitrogen is not suitable for nuclei extraction of small tissues, we chopped the tissue with blade for homogenization. In order to fix the suitable buffer condition to isolate nuclei in wheat, we tested three frequently used plant nuclei extraction buffers, PP buffer and D-Sorbitol buffer that worked nicely for maize (Tu et al., 2020; Sun et al., 2020), chopping buffer which had been widely used to extract nuclei from various plants (Lu et al., 2017; Maher et al., 2018; Lu et al., 2019; Zhu et al., 2020). We found that PP buffer was of low efficiency and generated cell debris ( Figure 1A ). Chopping Buffer produced more nuclei but the nuclei tended to aggregate with lots of cell debris ( Figure 1B ). The highest quality wheat nuclei were obtained using D-Sorbitol Buffer, not only regarding the number of nuclei but also the cleanest background as compared with the other two buffers ( Figure 1C ). Considering that high concentration of D-Sorbitol may affect later Tn5 tagmentation, we decreased the concentration of D-Sorbitol to 0.4 M (D-S-0.4). We then verified the applicability of D-S-0.4 buffer in different tissues, including aboveground part of 10-day old seedling, radical at 2 days after germination, double-ridge stage of spike (about 2 mm), and ovary at 4 days after pollination. As shown, the nuclei from different tissues are intact, which indicated that the D-S-0.4 buffer was suitable for nuclei extracting from different wheat tissues. However, different cellular structures and tissue specific components can be found in the suspension, such as root hairs from root tissue and starch grains from the ovary, but in general, cleaner nuclei could be obtained from tissues with meristematic activity, such as young inflorescence ( Figures 1D-G ).

In plant cells, besides the DNA fragments with open chromatin in the nuclei, there are exposed organelle DNA in mitochondria and chloroplasts, which can be cleaved by Tn5. Fluorescence activated nuclei sorting (FANS) sorts stained nuclei to get rid of other organelles. Combination of FANS with ATAC-seq (FANS-ATAC-seq) reduced the alignment rate of chloroplasts and mitochondria to 30% in Arabidopsis (Lu et al., 2017). To develop ATAC-seq protocol in wheat, we used simple chopping method followed by flow cytometry sorting to purify nuclei for Tn5 tagmentation ( Figure 2A ). In this method, the nuclei were stained with 4, 6-Diamidino-2phenylindole (DAPI) and then submitted to FANS. The nuclei can be divided into peaks of 2C and 4C. 2C represents the nuclei from somatic cells and 4C is for the dividing nuclei. We collected all the nuclei of 2C and 4C ( Figure 2B ). A total of 50,000 sorted nuclei were obtained for Tn5 tagmentation at 37°C for 30 min. The sequencing tag was inserted into the open chromatin through Tn5 and then the DNA was purified for library preparation. After tagmentation and PCR amplification, clear nucleosome pattern with different fragment sizes can be seen in the ATAC-seq libraries of aboveground of 10-day old seedling, ovary at two days after pollination and young spike at double-ridge stage. ( Figures 2C–D ), indicating proper digestion of Tn5. In order to evaluate the feasibility of FANS-ATAC-seq in wheat, especially with low amount of tissue input, we selected young spike at double-ridge stage and ovary at two days after pollination to generate ATAC-seq data ( Figure S1 ). Only with ten spikes and ten 2-day old ovaries, we successfully got more than 300,000 nuclei and 50,000 were used for subsequent experiment.

How many reads are needed is always a question for next-generation sequencing based methods. The sequencing depth is quite important for subsequent qualitative and quantitative analysis. Insufficient sequencing data may cause false positive results and lost positive signals due to the interference from background noise. Exceeded sequencing not only wastes money, but also reduces the signal-to-noise ratio. To explore the sequencing saturation of ATAC-seq in wheat, we sequenced one library from each of the young spike and the ovary for 16-17 Gb high quality clean data, which is similar to the size of wheat genome. Then, we randomly extracted from 10% to 200% of the reads from the clean data to call peaks, respectively. The results show that the saturation of sequencing depth is sample dependent, but in general, 16-17 Gb data is optimal for ATAC-seq in wheat ( Figure 3 ). Thus, we aimed for 16 Gb data for each of the four libraries.

ATAC-seq data analysis is challenging in wheat due to the presence of more than 80% repetitive sequences in the genome. We established a workflow and a pipeline to analyze ATAC-seq datasets in wheat and it can be easily extended to other plant species. To solve the problem of high ratio of repetitive sequences in wheat, we set the number of allowed mismatches to 0 when aligning the clean reads to the reference genome of IWGSCv1.0. The imperfectly aligned reads and reads from PCR duplication were removed. By using Bowtie2 for alignment and the “SAMtools index -c” command to build a csi index for bam files, the reads can be aligned directly to the reference genome without splitting the chromosomes, thus avoiding the hassle of position conversion. To reduce the false positive rate and obtain the high-quality, consistent, and reproducible peaks, the Irreproducible Discovery Rate (IDR) values between original replicates, the pseudo-replicates which respectively and randomly split from each of the original replicate, and the pseudo-replicates that randomly split after pooling the replicates were calculated by IDR v2.0.4.2 (https://github.com/nboley/idr), respectively ( Figure 4A ). The IDR value of 0.05, 0.02 and 0.01 are the recommended cutoffs for the original replicates, the pseudo-replicates split from original replicates and pooled pseudo-replicates, respectively for wheat ATAC-seq analysis. The whole procedure can be divided into five steps, including i) removing the adapters to obtain clean reads, ii) generation of unique mapped reads, iii) preparing files for IDR analysis, iv) calling peaks and v) combination of consistent peaks ( Figure 4A ).

A successful ATAC-seq should capture OCRs with high signal-to-noise ratio. In general, the most direct way for data estimation is to visualize peaks in genome browser like Integrative Genomics Viewer (IGV). We integrated the quality evaluation process in the pipeline by including the production of visualization files, calculation for library complexity, insert fragment distribution, sample cross correlation, TSS enrichment score, FRiP and the Signal Portion Of Tags (SPOT) values ( Figure 4B ). In addition, the functional feature of the peak region and the overlapping ratio of the peaks with other tissues are also indicators of data quality ( Figure 4B ). Both the data processing and date quality evaluation pipelines are available in Github (https://github.com/hcph/wheatATAC.git). The users can take raw sequence data as input to obtain OCRs and all the results for quality evaluation ( Figure 4B ). Since the parameters are easy to be adjusted, the pipeline is generally applicable to other plant species. Importantly, as the finally obtained signals would have been filtered by multiple comparisons, the open chromatin signals are reproducible and comparable across different libraries. In addition, the pipeline is based on Linux system and is highly flexible.

Using the data processing strategy we proposed, we found that at least 90% of the four datasets we obtained here from young spike and ovary can be aligned to the wheat genome with the highest genome alignment rate of 95% ( Table S2 ). We then used these data to predict OCRs in young spike and ovary of wheat. After preformed IDR analysis, we finally obtained 70,766 and 86,041 OCRs in young spike and ovary, respectively. These OCRs are consistent among replications. The number of OCRs in ovary is higher than that in young spike. The signal-to-noise ratios of the ATAC-seq data from both tissues were high ( Figure 5A ). The correlation co-efficiency between the replication of tissues is higher than 0.97, especially in young spike, with a co-efficiency of almost 1.00 (R score), indicating the reproducibility of our ATAC-seq protocol ( Figure 5B ). Not surprisingly, chromatin accessibility between young spike and ovary is lowly correlated.

We defined the 4 kb upstream to 1.5 kb downstream of TSS as promoter region and the regions beyond 4 kb upstream of TSS or downstream of transcription end site (TES) were treated as intergenic regions where distal regulatory elements like enhancers may exist. The distribution of the predicted OCRs in TSS and TES showed the same trend in young spike and ovary ( Figure 5C ). This is consistent with the result in Arabidopsis thaliana (Lu et al., 2017). However, different from the situation in Arabidopsis, more than 50% of the open chromatin can be found in intergenic region in wheat, which might indicate the presence of distal regulatory elements in these regions or it might be due to the super large genome size with large number of repetitive sequences and more complex spatial genome organization in wheat ( Figure 5D ).

A high correlation between chromatin accessibility and gene activity is expected since OCRs are more preferable to be bound by regulatory proteins to alter target gene transcription. In total, 29,909 genes in young spike and 33,452 genes in ovary were identified with ATAC-seq peaks in the promoter regions. We profiled the gene transcription of young spike and ovary. The transcription level of genes with ATAC-seq peaks is significantly higher than those without any ATAC-seq peaks ( Figure 6A ). Strikingly, the expression level is directly proportional to the number of peaks ( Figure 6B ).

We performed GO analysis on the relevant genes obtained from the tissue-specific ATAC-seq peaks of young spike and ovary, respectively. Genes with OCRs in young spike are mainly clustered in the pathways of protein-containing complex assembly, chromosome organization, protein folding and DNA conformation change ( Figure 6C ), which suggests that open chromatin mediates rapid cell division and differentiation in young spike by affecting protein complex assembly and DNA conformation. Most of the open chromatin related genes in ovary are in the coenzymes synthesis and nucleotide metabolism pathways ( Figure 6D ), indicating active metabolic reactions in the ovary. All these results above again approved the high quality of our ATAC-seq data from the aspect of biological sense.