A chromatin accessibility landscape during early adipogenesis of human adipose-derived stem cells.

Obesity has become a serious global public health problem; a deeper understanding of systemic change of chromatin accessibility during human adipogenesis contributes to conquering obesity and its related diseases. Here, we applied the ATAC-seq method to depict a high-quality genome-wide time-resolved accessible chromatin atlas during adipogenesis of human adipose-derived stem cells (hASCs). Our data indicated that the chromatin accessibility drastic dynamically reformed during the adipogenesis of hASCs and 8 h may be the critical transition node of adipogenesis chromatin states from commitment phase to determination phase. Moreover, upon adipogenesis, we also found that the chromatin accessibility of regions related to anti-apoptotic, angiogenic and immunoregulatory gradually increased, which is beneficial to maintaining the health of adipose tissue (AT). Finally, the chromatin accessibility changed significantly in intronic regions of peroxisome proliferator-activated receptor γ during adipogenesis, and these regions were rich in transcription factors binding motifs that were exposed for further regulation. Overall, we systematically analysed the complex change of chromatin accessibility occurring in the early stage of adipogenesis and deepened our understanding of human adipogenesis. Furthermore, we also provided a good reference data resource of genome-wide chromatin accessibility for future studies on human adipogenesis.

Introduction

Obesity has become a severe public health problem worldwide, often accompanied by metabolic disorders and chronic inflammation [1,2]. Patients with obesity are also at higher risk of other diseases, including type 2 diabetes, cardiovascular diseases and certain cancers [3]. Obesity is characterized by excessive accumulation of AT, which can be realized by expanding in adipocyte size (hypertrophy) or forming new adipocytes (hyperplasia). In physiological conditions, the AT plays not only a key regulator of lipid storage and release but also an important endocrine organ secreting many different hormones. In pathophysiological conditions of obesity, excess lipids accumulate in both the visceral depot and the ectopic sites that produce lipotoxicity in these tissues [4]. Along with the AT expansion also accompany hypoxic state, pro-inflammatory state, which become dangerous to health [2]. Clearly, a healthy and well-functional AT is of significant importance for health. The metabolic complications of both too much AT (in obesity) or too little (lipoatrophy/lipodystrophy) are quite similar deleterious. Adipogenesis, the process of adipose-derived stem cells (ASCs) converting into new adipocytes, plays an important role in AT expansion of obesity. Notably, some literature indicates that inhibited or promoted the process of adipogenesis all may contribute to reversing the metabolic disorders of obesity and also may accompany other adverse health effects [5,6]. Adipogenesis is a complex physiological process and many critical questions regarding the regulation of adipogenesis remain unresolved. Therefore, a deeper understanding of adipogenesis and whether adipogenesis confers a benefit to metabolically healthy individuals is meaningful [5].

Adipogenesis is a highly coordinated process that is accompanied by extensive transcriptional reprogramming and controlled by a series of tightly transcriptional regulations [7]. It is necessary for transcription regulation that transcription factors (TFs) interact with the transcriptional regulatory elements; thus, the open chromatin regions are a necessary condition for transcription regulation [8]. Numerous studies have shown that chromatin accessibility plays an essential role in establishing and maintaining cellular identity [9], but little is known about the systemic change of chromatin accessibility during human adipogenesis. A few studies have investigated the chromatin remodelling during adipogenesis on the murine 3T3-L1 cell line, which provides some reference for us to understand human adipogenesis [7,10,11]. Due to the considerable differences in genomic and epigenomic regulation between human and murine, the contribution of the results from the 3T3-L1 adipogenesis model to understanding the process of human adipogenesis remains controversial [12]. Therefore, it is necessary to systematically investigate the dynamic changes of chromatin accessibility during human adipogenesis.

ATAC-seq is an emerging high-throughput sequencing technology that defines and quantifies chromatin accessibility genome-wide by inserting the sequencing adapters into accessible regions of chromatin by Tn5 transposase [13]. Compared with traditional methods, ATAC-seq has become a more helpful tool for studying chromatin remodelling events in multiple biological processes due to its high sensitivity, high resolution and low starting cell number [14]. Here, we depicted a time series genome‐wide accessible chromatin atlas by the ATAC-seq method during adipogenesis of hASCs and discussed the complex changes of chromatin accessibility during adipogenesis of hASCs in detail. Overall, we provided a data resource of genome‐wide chromatin accessibility for future research of human adipogenesis and hoped to deepen our understanding of the chromatin remodelling process during human adipogenesis.

Results and discussion

Mapping the genome‐wide chromatin accessibility landscape during adipogenesis by ATAC-seq

To obtain the dynamic chromatin accessibility atlas during adipogenesis of hASCs, we applied ATAC-seq technology (Figure 1a). All hASCs were derived from ScienCell Research Laboratories and collected after obtaining written informed consent signed by donors. Our previous work showed that during adipogenesis of hASCs, chromatin remodelling occurred twice around the 8-hour point and two main function phases (commitment phase 0 h-8 h, determination phase 24 h-7d) were identified [15]. Therefore, we selected four time points (0 h, 2 h, 8 h and 24 h) in this study, hoping to reflect the complex chromatin accessibility changes from the lineage commitment phase to the terminal differentiation phase during adipogenesis of hASCs.

Figure 1.

Experimental methodology overview and data quality control. (a) Experimental outline of chromatin accessibility landscape construction during hASC adipogenesis. (b) The ATAC-seq signal enrichment around the TSSs for each sample. (c) The correlation of ATAC-seq peak reads between two replicates of each time point. (d) The fragment size distribution of ATAC-seq in each sample.

After ATAC-seq data collection, we preliminary evaluated our ATAC-seq data through FastQC and the results showed good data quality. We obtained an average of close to 70 million raw reads for each replicate, which was shown to be sufficient for the subsequent analysis. The observation of strong enrichment of ATAC-seq reads around transcription start sites (TSSs) in all samples (Figure 1b) further suggested the successful detection of accessible regions in our result [16]. To further evaluate the reproducibility and robustness of our result, we compared the similarity of open chromatin profiles across samples by Pearson correlation coefficients (Figure 1c). The chromatin accessibility profiles were highly correlated between two replicates (R = 0.88–0.98). Besides, we also investigated the fragment size distributions in each sample. In agreement with typical ATAC-seq profiles [17], the chromatin accessibility fragments show size periodicity corresponding to integer multiples of nucleosomes, and the fragment size distributions between biological replicates were highly correlated (Figure 1d). Overall, all these results above indicated that we obtained high-quality chromatin accessibility data by ATAC-seq and the data was robust and reproducible.

Among assays designed for detecting chromatin accessibility, ATAC-seq since first described in 2013 has gained particular popularity and becoming increasingly common which indicates its value in a wide spectrum of biological questions [18]. ATAC-seq can probe chromatin accessibility states with low cellular input and obtain a higher signal-to-noise ratio in comparison with traditional methods such as DNase, FAIRE or MNase-seq [19]. Despite the advantages, ATAC-seq also has some limitations. One of the limitations is that the classic ATAC-seq method relies on its two barcodes based library construction, which causes a part of the input tagmentation fragments to be lost after PCR amplification [20]. Several single-stranded DNA library preparation methods have already been used to address this issue [21]. A prominent example is SALP, which is overcoming the loss of chromatin accessibility information and achieving library construction simplicity and efficiency by using a novel kind of single-strand adaptor. It is worth noting that SALP is just dependent on two cheap routine enzymes, T4 DNA ligase and Taq polymerase, which greatly increases the usability of the method [22]. In addition, to the applicability limitations of different types of samples and single-cell samples, ATAC-seq-derived technologies generated include fast-ATAC, Omni-ATAC, miniATAC-seq, scATAC-seq (10 × ATAC-seq) and sci-ATAC-seq [23]. ATAC-seq and ATAC-seq-derived technologies wide applications in the future will help us gain a more comprehensive and in-depth understanding of chromatin dynamics and their biological significance.

The macro-dynamics changes of chromatin accessibility during adipogenesis of hASCs

Next, we globally analysed the dynamic changes of chromatin accessibility during adipogenesis of hASCs. Statistical analysis of accessible chromatin peaks revealed that the number of chromatin accessible peaks gradually increased from 0 h to 8 h and then declined from 8 h to 24 h (Figure 2a). Consistent with this result that our previous work demonstrated, hASCs experienced two chromatin remodelling before and after 8 h point during adipogenesis [15]. Our prior work also showed an obvious phase mode, commitment phase (0 h ~ 8 h) and determination phase (8 h-7d). Thus, we argue that chromatin accessible peaks gradually increased during the commitment phase (0 ~ 8 h) and then declined during the determination phase (after 8 h) at a genome level. Adipogenesis is the process that hASCs lose their developmental plasticity gradually and then take on specialized functions of adipocytes. Compared with hASCs, it seems a smaller number of open chromatin regions were sufficient to maintain adipocyte identity. The same phenomenon was observed in the differentiation of human embryonic stem cells into dopaminergic neurons that open chromatin were regions decreasing [24]. Based on the above results, 8 h may be the critical transition node of adipogenesis chromatin states from the commitment phase to the determination phase.

Figure 2.

Overview of chromatin accessibility variation during hASC adipogenesis.

We further merged the fragment size distribution curve between each sample for comparison (Figure 2b). Data showed that the mode of fragment size distribution at 24 h was significantly different from the others and the higher proportion of mono-nucleosome-bound fragment peak means that the chromatin accessibility reduced at 24 h. Principal components analysis (PCA) also exhibited significant differences among different time points (Figure 2c). Moreover, we also visualized genome-wide chromatin accessible peaks among all samples. This circle diagram clearly displayed the radical and global difference of chromatin accessibility between time points (Figure 2d). To annotate the genomic function of accessible chromatin peaks, we used ChIPseeker, a suite of tools for sequencing data analysis. The proportion of intron regions in chromatin accessible regions obviously increased from 0 h to 8 h and reduced from 8 h to 24 h. In contrast, the proportion of promoter regions in chromatin accessible regions reduced from 0 h to 8 h and increased from 8 h to 24 h (Figure 2e). All results above demonstrated that the chromatin accessibility and the proportion of gene function element distribution were changed significantly and dynamically during the early stage of hASC adipogenesis. In addition, the dramatic chromatin accessible change of promoter regions and intron regions suggests the important role of these regions during adipogenesis of hASCs.

Cluster analysis revealed five distinct chromatin accessibility change modes during adipogenesis of hASCs

To identify the specific dynamic change pattern in the accessible chromatin peaks landscape during adipogenesis of hASCs, we employed TCseq [25]. Within the TCseq framework, we chose C-means Clustering. We identified five specific Clusters with different accessibility change behaviours across the four-point time course (Figure 3a). The variation of peak amplitudes in each cluster reflected transient opening response (Cluster I), late opening response (Cluster II), gradual opening response (Cluster III), gradual closing response (Cluster IV) and late closing response (Cluster V) (Figure 3a). It is noteworthy that Siersbæk et al. mapped extensive chromatin remodelling during early adipogenesis of 3T3L-1 by DNase I Hypersensitive site analysis but obtained an obviously different result as ours. An obvious difference, there were only opening response clusters but no closing response clusters during early adipogenesis (0 h-24 h) in their experiment [11]. This result may reflect the difference in chromatin remodelling between human and mouse adipogenesis models.

Figure 3.

Dynamic temporal variation mode of chromatin accessibility during hASC adipogenesis. (a) Five dynamic accessibility behaviour modes of chromatin accessibility data analysis by c-means cluster during adipogenesis. (b) GO analysis of the regions in each cluster by GREAT.

To further identify the characteristics of these clusters, we employed the Gene Ontology (GO) analysis to each cluster by GREAT programme [26] (Figure 3b). In general, in eukaryotic cells, at gene cis-acting elements such as promoters and enhancers, chromatin needs to be in an accessible state for TF binding to activate target gene [27]. However, condensed chromatin always results in gene silencing [8]. Therefore, higher chromatin accessibility means higher transcription activity [28]. As we expected, GO analysis showed that many developments and differentiation-related regions (Cluster II and Cluster III) gradually opened during adipogenesis. Compared to Cluster III, regions related to the terminal differentiation of adipocyte (including signalling pathway of peroxisome proliferator-activated receptor, mitochondrion distribution and cholesterol metabolic [29]) were only enriched in Cluster II. The observation above suggests that the regions related to promoting adipogenic differentiation gradually opened after induced and the region mainly related to terminal differentiation opened after 8 h, which was also consistent with our previous work [15]. Interestingly, many regions related to the regulation of pri-miRNA transcription were also significantly enriched in Cluster II. Many miRNAs have been proved to play regulatory roles in adipogenesis [30]. Our results pointed out that the accessible regions related to the regulation of pri-miRNA may be important in the terminal differentiation of adipogenesis.

Regions in Cluster IV and Cluster V gradually lost accessibility during adipogenesis, and the lower chromatin accessibility means a lower transcription activity. GO analysis showed that many anti-adipogenesis-related regions gradually closed during adipogenesis (Cluster II and Cluster III), as we expected [7]. A typical example was the PKC-related process, which has been widely proven in human and murine [31-33]. The regions related to proliferation, cell differentiation, systems development and cell fate commitment also lost accessibility from 8 h to 24 h and showed that chromatin state trans to terminal differentiation phase after 8 h.29 Interestingly, we found that regions related to negative regulation of the phospholipid metabolic process gradually closed (Cluster IV) and the regions related to positive regulation of organophosphate catabolic process gradually opened (Cluster III). This hinted that ensuring the up-regulation of the phospholipid metabolic process might be important in the early process of adipogenesis.

Different from other clusters, Cluster I was temporarily open. GO analysis showed that the regions in Cluster I were highly enriched in the biological process of cell differentiation regulation. Our result implied that the temporarily opening of the regions in Cluster I might be necessary for adipogenesis of hASCs. For example, the regions related to the regulation of cellular response to platelet-derived growth factor (PDGF) stimulus cytokine secretion were enriched in Cluster I. The PDGF-related process is traditionally recognized as anti-adipogenic. Artemenko et al. found that PDGF directly inhibited murine 3T3‐L1 preadipocyte and human preadipocyte differentiation [34]. Studies also showed that members in PDGFRs were also inhibitors of adipogenesis [35,36]. Interestingly, the stimulus of PDGF also may promote adipogenesis in recent research. PDGF -BB enhances adipogenesis in orbital fibroblasts [37]. Moreover, de novo adipogenesis in Matrigel also demonstrated that PDGF was a potent inducer as important as BFGF in de novo adipogenesis during many candidates [38]. These seemingly contradictory results implied that the short activation of the regions that respond to PDGF stimulus may be necessary for adipogenesis and also supported the potential importance for the temporarily opening of the regions in Cluster I during adipogenesis.

We also found that regions related to the anti-apoptotic process(including negative regulation of extrinsic apoptotic signalling pathway, negative regulation of extrinsic apoptotic signalling pathway in absence of ligand) gradually opened (Cluster III). Apoptosis is an important homoeostatic mechanism controlling adipocytes number, and excessive adipocyte apoptosis appears to be of some concern related to detrimental metabolic effects of fat mass loss that will increase blood lipid concentrations, ectopic lipid storage and other detrimental metabolic effects [39]. Thus, tight control of apoptosis induction is necessary to maintain a healthy AT. The study by Sorisky et al. confirmed that adipocytes acquire a relative resistance to apoptosis as the process of adipogenesis which supports our result [40]. Although several pieces of evidence have shown that the activation extrinsic apoptotic signalling pathway regulates adipogenesis [39], the function of an extrinsic apoptotic signalling pathway in absence of ligand in adipogenesis has not been elucidated. Extrinsic apoptotic in absence of ligand involves a special variety of receptors, called dependence receptors or addiction receptors. Dependence receptors activate apoptotic pathways following the withdrawal of trophic factors and other supportive stimuli [40]. It has been shown that dependence receptors play a major role in organogenesis, but there is no evidence that dependence receptors play a role in adipogenesis. Our result also hinted that the apoptotic process with the participation of dependence receptors might regulate adipogenesis.

It is also worth noting that a large amount regions related to the T cell differentiation and activation were significantly enriched in Cluster I, Cluster II and Cluster III. T cells are the second-largest cell population in AT and are also considered to play an important role in obesity-induced inflammation and metabolic homoeostasis disorder [41]. In the lean state, AT is dominated by anti-inflammatory Th2 and Treg cells, which help to maintain an anti-inflammatory milieu and metabolic homoeostasis. With an excessive nutrient intake and/or reduced energy expenditure, the total number of CD4 + T cells is significantly increased, and the populations of Th2 and Treg cells are decreased, concurrently with increased Th1 and Th17 cells, which led to a state of chronic low-grade inflammation [42]. Hence, strict control of T cell homoeostasis is required to maintain a modest immune response in AT. In our study, pro-inflammation regulation processes (such as positive regulation pro-inflammation T cell differentiation and activation, interleukin-1-mediated signalling pathway [43]) related regions gradually opened (Cluster II and Cluster III). At the same time, anti-inflammation and the immunomodulatoryprocess (such as NK T cell differentiation, interleukin-2 production, transforming growth factor beta production, negative regulation of immune response) related regions gradually opened (Cluster II and Cluster III) [44-46]. The simultaneous opening of anti-inflammatory and pro-inflammatory-related regions suggests that adipogenesis activated the immunoregulatory functions of adipocytes. This might partly explain why new adipogenesis can reverse obesity-associated inflammation [5]. Traditionally, as an endocrine cell, adipocyte can greatly modulate an inflammatory response by secretion of adipokines, cytokines and chemokines [47]. Recent evidence has shown that extensive crosstalk exists between adipocytes and immune cells [42]. Our results support that adipocytes might serve as immune regulatory cells in AT inflammation through modulation of conventional immune regulatory cell types [48]. Furthermore, many regions that were related to the negative regulation of the differentiation and activation of T cell (T-helper cell, CD4-positive T cell and alpha-beta T cell) temporarily open (Cluster I) during the early stage of adipogenesis. Evidence suggests that pro-inflammatory cytokines secreted by T cells inhibit adipogenesis through various pathways [49] which showed that the temporary activation of the regions that negative regulation of T cells differentiation and activation may ensure the smooth adipogenesis process during the early stage of adipogenesis.

The growth of AT requires continuous remodelling of the vascular network which ensures close vascular endothelial cell-adipocyte crosstalk and the dynamic metabolic balance of adipocyte [50]. In our result, many regions related to the angiogenic were significantly enriched in Cluster II, Cluster III and Cluster V. Insufficient vascularization during AT expansion promotes hypoxia causing the release of inflammatory cytokines and chemokines, which may trigger further complications including inflammation, fibrosis and apoptosis, contributing to AT dysfunction [51]. Our result showed that many promotion angiogenesis related regions (Cluster II and Cluster III) gradually opened during adipogenesis, which includes vascular endothelial growth factor and fibroblast growth factor 2 (FGF2) that the most important factor coordinated vessel growth and increased [50]. In addition, many regions related to endothelial cell differentiation gradually lost accessibility during adipogenesis (Cluster V). Evidence that hASCs and pre-adipocytes could direct transdifferentiate into endothelium in the AT [52]; therefore, inhibition of auto-endothelial differentiation is necessary to promote peripheral angiogenesis while ensuring auto-adipogenesis.

Rich chromatin accessibility variation at PPARG region during adipogenesis of hASCs and abundant TF binding motifs at these variation regions

Adipogenesis is a highly orchestrated process tightly controlled by a complex transcription regulation network. Even though hormones, miRNAs and many other factors can regulate adipogenesis, there is no doubt that TFs are the most important regulator of adipogenesis [7]. TFs need to recognize and bind specific DNA sequences (motif) to control chromatin remodelling and transcription. Open, partially nucleosome-free chromatin regions are much more accessible to be bound by TFs than closed nucleosome-dense regions [53]. Therefore, the accessibility change of the TF binding motifs is important for the regulation of TFs in adipogenesis. We widely scanned the TF binding motifs in the regions of each Cluster using HOMER. Then, we sorted and displayed the motifs of TFs according to the P-value. Many TF binding motifs were located in the variable accessible chromatin regions and almost all the TF binding to these motifs have been reported to play roles in adipogenesis [7] (Figure 4a).

Figure 4.

Rich TF binding motifs distribution at chromatin accessibility variation regions of PPARG locus during adipogenesis of hASCs. (a) Enrichment of TF binding motifs for the dynamical variation regions of chromatin during adipogenesis of hASCs. (b) TF binding motifs distribution of chromatin accessibility variation regions of PPARG locus.

Peroxisome proliferator‐activated receptor γ (PPARG) is widely considered as the most central regulator of adipogenesis that functions have been widely confirmed in in vitro or in vivo investigations [12,54]. Many other TFs regulating adipogenesis also play their roles by directly regulating PPARG [55]. Therefore, we focused on the motif distribution in the various regions of PPARG (Figure 4b). Regions located at PPARG enriched in abundant variable accessible chromatin regions and the change mode of various chromatin accessibility regions at PPARG mainly belong to Cluster II and Cluster III. Abundant TF binding motifs were distributed at these various regions, and many of these TFs had been verified that can play roles in adipogenesis through PPARG [56]. The chromatin accessibility change mode of binding motifs of TFs (including JUNB, FRA1, FRA2, CEBPB, BATF, ATF3 and PU.1) at PPARG belonged to Cluster II and Cluster III. The chromatin accessibility change mode of binding motifs of CJUN and BACH2 at PPARG only belonged to Cluster II. It is worth noting that most of these TF binding motifs were only located in intron regions rather than in exons or promoters of PPARG, which hinted that many TFs might subtly regulate the intron regions of PPARG during the early stage of adipogenesis. Although researchers always used to pay more attention to the TF binding sites in promoter regions. Nevertheless, increasing data show that TF binding sites in the intron also play critical roles in gene regulation [56-59]. The above results indicated that during the early stage of adipogenesis, many nascent accessible chromatin regions gradually appeared in the intron region of PPARG and numerous potential TF binding sites were exposed for further possible regulation.

Besides, TFs in CCAAT-enhancer-binding proteins (CEBPs) family (including CEBPA, CEBPB and CEBPD) also play key roles in adipogenesis [29]. Interestingly, unlike PPARG, there was no significant change in the chromatin accessibility state of CEBPs family. Among them, the regions of CEBPB and CEBPD maintained a high degree of accessibility, but the regions of CEBPA maintained a low degree of accessibility. Maintaining an open state is possibly conducive to the rapid response of CEBPB and CEBPD to the stimulus of adipogenesis [60]. Although CEBPA and PPARG are all playing the most prominent roles in the terminal differentiation phase of adipogenesis [61], the chromatin accessibility variation of PPARG regions is obviously apparent earlier than CEBPA after the stimulus of adipogenesis. This also further reminds the core role of PPARG in the process of human adipogenesis.

In general, we depicted a high-quality chromatin accessibility atlas by ATAC-seq method during adipogenesis of hASCs and systematically illustrated a complex change of chromatin accessibility occurring in the early stage of adipogenesis (0 h-24 h) in this study. We also found that, upon adipogenesis, the chromatin accessibility of regions related to anti-apoptotic, angiogenic and immunoregulatory gradually increased, which is beneficial to maintaining the health of AT and also provides new ideas for the treatment of human pathological obesity. It is worth noting that the chromatin accessibility in many intronic regions of PPARG changed significantly and quickly after adipogenesis induction and these regions were enriched in numerous TF binding motifs. These TF binding motifs were exposed for further regulation in the early stage of adipogenesis and may become novel targets to regulate adipogenesis. Finally, our ATAC-seq data provided a good reference database for future studies on human adipogenesis.

Materials and methods

Cell culture and differentiation

The hASCs were proliferated in Mesenchymal Stem Cell Growth Media (#HUXMD-90011, Cyagen) and induced differentiation to adipocytes by Adipogenesis Basal Medium (#HUXMD-90031, Cyagen) as described previously [15]. Cells were respectively collected at 0 h, 2 h, 8 h and 24 h after induced adipogenesis. Two biological replicates were conducted for each time point.

ATAC-seq sample preparation and data collection

Briefly, we used the Nextera DNA Library Preparation Kit (#FC-121-1031, Illumina) to execute the transposition after cells were harvested and lysed in lysis buffer. The transposed DNA fragments were purified by the MinElute PCR Purification Kit (#28,004, Qiagen). Then, samples were amplified by using 1X NEBNext High-Fidelity PCR Master Mix (#M0541L, New England Biolabs), the subsequent library was purified by the MinElute PCR Purification Kit (#28,004, Qiagen) and subjected to sequencing on Illumina Novaseq 6000 using PE150.

Data quality control and preprocessing

Raw data (Fastq files) was processed by FastQC v0.11.8, including the GC content, sequence length distribution, sequence duplication levels, contamination of primers and adapters [18]. The sequencing adaptors and poor-quality reads were detached by Trimmomatic v0.36 with the parameter ‘ILLUMINACLIP:adapter.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36’ [62], clean reads of high quality were obtained and the clean reads were aligned to hg19 reference genome sequences by bwa v0.7.17-r1188 with parameter ‘-v 3’ [63]. Finally, We used MACS2 with parameter ‘macs2 callpeak – nomodel -f BAMPE – keep-dup 1 -q 0.05 -B – SPMR’ [64] for peaks calling (the initial threshold q-value 0.05 as the cut-off value) and the peaks were used in the downstream analyses. The Pearson correlation coefficient (R) between two replicates was calculated by R package and the ATAC-seq signal enriched around the TSSs was handled by Deeptools v 3.5.0 [65].

Data analysis

ATAC-seq signal tracks were visualized through Integrative Genomics Viewer (IGV v2.5.0) [66]. The normalized read counts against all merge peaks were used for PCA, and the PCA plot was output by R package. Peaks were annotated by ChIPseeker [67]. The unsupervised Clustering R package, TCseq, was used to investigate temporal patterns of the time course sequencing data, the parameter is algo = ‘cm’, k = 5,standardize = TRUE [25]. GO analysis was performed by GREAT programme [68]. The HOMER’s findMotifsGenome.pl tool was used for motif analysis [68].

Circos analysis

We displayed the radical and global difference of chromatin accessibility between time points by circos plot with R package circlize(circos.initializeWithIdeogram(species=‘hg19’), circos.trackPlotRegion and circos.trackLines with default parameter).