Genome-Wide Characterization of the Accessible Chromatin Regions in Chickpea

1. Introduction

Chromatin, a complex structure composed of DNA and proteins, is formed by the binding of DNA to histones to form nucleosomes, which are further folded and compacted [1-3]. These structures undergo dynamic changes during DNA replication or transcription, allowing the tight structure of chromosomes to become open [4, 5]. The accessible chromatin regions (ACRs) in the genome facilitate the binding of trans-acting factors to the cis-regulatory elements present in promoters, enhancers, insulators, and silencers. This property, which allows for the binding of transcription factors (TFs), is defined as chromatin accessibility [6-8]. In addition to being regulated by TFs, the chromatin is also subject to various histone modifications, which affect the chromatin structure and the subsequent gene expression [9, 10]. Currently, ATAC-seq (Assay for Transposase-Accessible Chromatin with high throughput sequencing) [4, 11], as a common method for ACR analysis, uses the transposase Tn5 [12], an enzyme capable of labeling ACRs, to fragment these sequences for library construction to improve the understanding of the regulatory mechanisms of gene expression across plant genomes.

Cicer arietinum (chickpea) is an herbaceous plant, belonging to the genus Cicer in the family Fabaceae. It is an important food crop in China, being rich in protein, various amino acids, vitamins, crude fiber, calcium, magnesium, iron, and other ingredients, and has a high economic value and extensive application prospects [13, 14]. However, studies on regulatory sequences in chickpea have remained limited. In this study, we have adopted ATAC-seq to identify the ACRs across the chickpea genome. Integrated with other data, our results provided a comprehensive understanding of the genomic features of these important cis-regulatory sequences and their impact on downstream gene expression. Furthermore, our data provided genomic resources for further exploration of non-coding regulatory ACRs for genetic diversity and crop improvement in chickpea.

2. Materials and Methods

2.1. Plant Materials

Chickpea (Cicer arietinum) were grown naturally in a greenhouse, and leaf samples were collected 15 days after germination for subsequent exprements.

2.2. RNA-sequencing (RNA-seq) data analysis

The RNA-seq raw data of chickpea leaves were downloaded from NCBI under accession number (PRJNA742483). Adapter trimming of the raw data was performed using TrimGalore (http://github.com/FelixKrueger/TrimGalore, default parameters), and then mapped to the chickpea reference genome [15] using the HISAT-StringTie pipeline [16]. Gene expression levels were calculated using the default parameters and expressed as fragments per kilobase per million mapped reads (FPKM).

2.3. Assay for Transposase-Accessible Chromatin Using Sequencing (ATAC-seq)

Leaf ATAC-seq library construction was performed as previously reported [17, 18]. Briefly, approximately 0.3 g of chickpea leaf sample was flash-frozen in liquid nitrogen, and then rapidly ground to powder. Nuclei were extracted by adding 10 mL of lysis buffer, followed by a washing step with 1 mL of Wash Buffer I and Wash Buffer II, respectively. Subsequently, about 105 nuclei were then resuspended in 50 μL of Tn5 reaction solution (TruePrep DNA Library Prep Kit V2 for Illumina, Vazyme, TD501) and incubated at 37 °C for 30 min. The tagged DNA was purified using the DNA Clean & Concentrator Kit (Zymo Research, D4014) and amplified using Q5 DNA polymerase (NEB, M0491) with indexed primers from the TruePrep Index Kit V2 for Illumina (Vazyme, TD202). The library was sent for sequencing on the Illumina NovaSeq PE150 platform.

Raw ATAC-seq data were filtered by TrimGalore and mapped to the chickpea genome by Bowtie2 [19] with the parameter set as: -X 1000. Samtools [20] was used to extract the reads with the MAPQ value above 30, which were further used for peak calling by MACS3 (parameters set as “--nomodel --shift 100 --extsize 200”) [21]. Chickpea genomic DNA, processed by the Tn5 enzyme, was used as an input control. The reproducible peaks between two replicates were defined as ACRs. Annotation of the ACRs was performed using HOMER (http://homer.ucsd.edu/homer). Visualization of ATAC-seq signals was performed using Deeptools [22] and IGV [23]. MEME Suite [24] was used for motif analysis.

2.4. Histone Chromatin Immunoprecipitation Sequencing (ChIP-seq)

Chickpea histone ChIP-seq was performed as described in the previous study with minor modifications [25-27], using ChIP-seq grade anti-H3K27ac (PTMA, PTM-116) and anti-H3K4me3 (PTMA, PTM-613) antibodies. Approximately 5 ng of ChIP-seq DNA and input DNA were utilized for library construction using the TruePrep DNA Library Prep Kit V2 for Illumina (Vazyme, TD501). Aptamer trimming and removal of low-quality reads was performed using Trim_Galore, and mapping of reads to the genome was performed using Bowtie2 [19]. The input library was used as a control. Mapped reads (MAPQ >30) were peaked using MACS3 with the narrow peak model [28].

3. Results

3.1. Identification of ACRs in the Chickpea Genome by ATAC-seq

In this study, ATAC-seq was employed to identify and characterize ACRs across the chickpea genome. Leaf ATAC-seq libraries from two biological replicates of chickpea generated approximately 99.6 million (M) and 100.8 M reads, respectively (Figure S1, Table S1). As expected, ATAC-seq signals were mainly enriched around the transcription start site (TSS) and exhibited a positive correlation with the expression levels of coding genes (Figure 1A and 1B), indicating that the chromatin accessibility is associated with gene activation, which is consistent with previous reports. A total of 11,555 ATAC-seq enriched regions were identified in two replicates and defined as ACRs, which were found to be associated with 7,609 coding genes (Figure 1C, Table S2). Genomic annotation of these ACRs revealed that the majority of ACRs were in the promoter (35%) and the distal intergenic regions (42%), while 12% of ACRs were observed in the downstream region and 11% in the exonic region (Figure 1D, Table S2).

3.2. Enrichment Analysis of Conserved Transcription Factor Motifs Enrichment in the ACRs

TFs can regulate gene expression by recognizing specific DNA sequences. During chromatin opening, TFs bind to DNA sequences (e.g. Motifs) in the open chromatin regions. Therefore, ACRs tend to enrich for conserved TF binding motifs [29, 30]. To assess the enrichment of TF binding motifs on chickpea ACRs, the TF binding motifs reported in the PlantPAN3.0 database were mapped to these ACRs [31]. The result showed that conserved TF binding motifs were significantly enriched on ACRs, with the top 10 enriched motifs being AP2/ERF, bHLH, bZIP, C2H2, Dof, MADS_box, MYB, NAC, TCP and WRKY (Figure 2A and 2D, Table S3). Subsequently, a potential TF regulatory network was then constructed using TF binding motifs enriched in ACRs and genes related to the downstream of ACRs (Figure 2B). In addition, a total of 5,881 (77.5%) genes were observed to be co-targeted by the MYB, NAC, and AP2/ERF families (Figure 2C). Taken together, these results suggest that there are several conserved TF-binding motifs in chickpea ACRs that could synergistically regulate the expression of downstream genes.

3.3. histone Modifications in ACRs

It has been reported that histone modifications are usually associated with ACRs to regulate gene expression [32-35]. In this study, ChIP-seq was adopted to investigate the distribution of H3K27ac and H3K4me3 in chickpea leaf ACRs and the relationship between these histone modifications and ACRs. A total of 32,260 H3K27ac-enriched peaks and 22,164 H3K4me3-enriched peaks were identified in chickpea (Tables S4 and S5). H3K4me3- and H3K27ac-modified ACRs were mainly enriched in the promoter region, and positively correlated with gene expression as expected (Figure 3A and S3, Table S4 and S5).

Using the K-mean algorithm, the histone-modified ACRs could be classified into two major clusters (Figure 3B). The results indicated that the majority of these ACRs were modified by H3K4me3, with a smaller number modified by H3K27ac in Cluster I, whereas most of the modified ACRs in Cluster II exhibited minimal modified signal. Furthermore, it was observed that 743 and 1,736 ACRs were modified only by H3K4me3 and H3K27ac, respectively, with 3,565 ACRs being dual-modified by these two histone marks (Figure 3C and 3D, Table S6). Interestingly, in the dual-modified ACRs regions, the overall signal coverage of H3K4me3 was higher than that of H3K27ac (Figure 3E). The expression of genes associated with dual-modified ACRs was higher than that of those with only a single modification (Figure 3F, Table S6). Taken together, these data suggest that dual modifications of H3K27ac and H3K4me3 at ACRs may play an important role in chromatin opening and gene activation.

Distal ACRs (dACRs), located in intergenic regions, might act as enhancers and also play an important role in regulating gene expression [36, 37]. In this study, a total of 3,639 dACRs were identified and mainly distributed in the 2-8 kb intergenic region in chickpea. In addition, 688 dACRs were detected in regions more than 20 kb away from the nearest TSS (Figure 4A, Table S2). In contrast to the proximal ACRs (pACRs), it was found that the majority of the dACRs (77.08%) were found to be unmodified by either H3K4me3 or H3K27ac ((Figure 4B and 4E, Table S4 and S5). The dACRs with H3K4me3 or H3K27ac modification were located in the regions close to the TSS (<10 kb) (Figure 4C). As previously reported, the expressions of genes associated with dACRs and pACRs were more significant than those with dACRs or pACRs only, indicating an enhancer role of dACRs in gene expression regulation (Figure 4D).

3.4. Effects of ACRs Absent on the Homologous Gene Expression in Chickpea

Gene duplications are prevalent during whole genome duplication, resulting in the divergence of the homologous genes (hGenes). Indeed, during the duplication of hGenes, the gain and loss of their corresponding cis-regulatory sequences, such as the promoter, could also affect their neofunctionalization [38, 39]. The hGene pairs were derived from the ACR-associated genes (E-value < 1e-5) (Table S7). Both hGenes associated with ACRs showed no significant difference in their gene expression, while the expression of the hGene without ACRs was significantly lower than that of their hGene with ACRs (Figure 5). Furt hermore, the absence of ACRs in hGenes also affected the histone modification status (Table S7). This observation suggests an important role for ACRs in maintaining hGene expression during the progress of hGene duplication.

4. Discussion

In this study, ATAC-seq has been employed for genome-wide characterization of ACRs in chickpea. The results revealed that ACRs were positively correlated with their associated gene expression, and were also enriched with many conserved TF motifs, suggesting that ACRs are important for gene regulation by acting as cis-regulatory sequences. In addition, ACRs were shown to be associated with active histone marks (e.g. H3K4me3 and H3K27ac) in the regulation of gene transcription. Finally, the presence or absence of ACRs could influence the expression levels of hGenes. Taken together, these data could be considered valuable genomic resources for future studies on gene regulation in chickpea.

ATAC-seq has been widely used as an efficient tool for gene regulation studies, such as TF-DNA interaction analysis and promoter-enhancer prediction in plants. For example, HSFA1a played an important role in the dynamic formation of promoter-enhancer contacts near the ACRs in response to heat stress in tomato, by ATAC-seq and other methods [40]. A total of 446 regeneration-related core TFs are associated with the open chromatin regions in wheat [7]. Moreover, ATAC-seq is currently being used to study a variety of plants, such as Arabidopsis thaliana [11, 41], maize [42-44], cultivated soybean [17, 43, 45-48], rice [49], grape [50], tomato [51] and the medicinal plant, Andrographis paniculata [52].

The function of distal ACRs has received much attention in recent investigations, as they may serve as enhancer-like elements that could mediate contact with proximal cis-regulatory sequences to promote gene expression [53, 54]. Similar to the proximal ACRs, the distal regions also contained conserved TF-binding motifs (e.g. AP2/ERF, bHLH, bZIP, etc.), through which the expression levels of genes could be precisely controlled, thereby affecting various biological processes in plants [55]. Therefore, several reports have indicated that variants in these non-coding regions, such as patch chromatin variants, single nucleotide polymorphisms, indels, etc., would affect the domestication of the plants and their final phenotype [39, 45], which has been reported in plants, such as wheat [21], tomato [38], maize [42], etc. In this study, several cases of gain and loss of hGenes were observed, which significantly affected gene expression. It was suggested that the non-coding ACRs may also play an important role in the adaptation and domestication of chickpea, but further investigation is required to confirm this hypothesis.

5. Conclusions

In this study, we have provided a comprehensive overview of ACRs and their genomic features across the chickpea genome, using ATAC-seq and histone ChIP-seq methods. Furthermore, this study has shown that ACRs play important roles in gene regulation and potential roles in plant domestication, and our data provide an important resource for further exploration of the chickpea genome.