Reorganization of the DNA replication landscape during adipogenesis is closely linked with adipogenic gene expression

The mouse embryonic fibroblast 3T3-L1 cells can be induced to differentiate into adipocytes upon hormonal stimulation and are widely used to study obesity, which is an increasing health problem worldwide and one of the main causes of metabolic syndrome (Afshin et al., 2017; MacDougald and Lane, 1995). In vitro adipocyte differentiation (adipogenesis) is initiated by the treatment of confluent monolayer G0-phase cells with 3-isobutyl-1-methylxanthine (IBMX), dexamethasone (DEX) and insulin. Within a few hours of this hormonal treatment, the expression of the genes encoding c-Myc (Myc), c-Fos (Fos) and c-Jun (Jun) increases, followed by the expression of the genes encoding the CCAAT/enhancer-binding proteins β and δ (Cebpb and Cebpd) (Cao et al., 1991; Cornelius et al., 1994; Rosen et al., 2000; Tang and Lane, 2012). The genes encode transcription factors and activate many downstream genes that contribute to the adipocyte phenotype (MacDougald and Lane, 1995; Yeh et al., 1995). Continuous culture of these differentiation-initiated cells for up to 2 weeks gives rise to mature adipocytes that accumulate lipid droplets.

During the early stage of adipogenesis after hormonal stimulation, quiescent 3T3-L1 cells restart the cell cycle and undergo only one or two cell divisions, called mitotic clonal expansion (MCE), before differentiation into the mature adipocyte stage (Tang et al., 2003a,b). MCE is an essential step for full adipogenesis because forced inhibition of MCE by treatment with U0126 (a MEK inhibitor) or roscovitine (a CDK inhibitor) prevents cells from terminal differentiation into adipocytes, even in the presence of hormonal stimulation (Tang et al., 2003b). The indispensable role of MCE in adipogenesis is further supported by the knockdown of genes encoding proteins such as factor for adipocyte differentiation 24 (Fad24 or Noc3l) and 49 (Fad 49 or Sh3pxd2b) (Hishida et al., 2008; Tominaga et al., 2004). MCE is thought to be a window of opportunity for chromatin reorganization, although the changes that occur during this period are not well understood.

In the present study, we performed a detailed genome-wide analysis of replication timing (RT) to elucidate changes in the DNA replication program during the S phase of MCE. We found a number of RT changes, including specific changes in the first and second rounds of MCE. This analysis combined with transcriptome data revealed that these changes are closely related to changes in adipocyte gene expression.

To analyze RT regulation during MCE, we first investigated the cell cycle progression of 3T3-L1 cells upon hormonal stimulation. Based on 5-ethynyl-2′-deoxyuridine (EdU) incorporation, we estimated that cells undergoing the first and second rounds of S phase peak at approximately 18–21 h and 40–41 h after hormonal stimulation, respectively (Fig. S1, left). This is supported by the accompanied increase in total cell number (Fig. S1, right). We then measured genome-wide RT of cells undergoing MCE (for simplicity, the first and second rounds of MCE are designated as ‘MCE1’ and ‘MCE2’, respectively) by E/L Repli-seq (Hayakawa et al., 2021; Marchal et al., 2018) and compared it with that of parental control 3T3-L1 cells. In E/L Repli-seq, 5-bromo-2′-deoxyuridine (BrdU)-immunoprecipitated DNA from fluorescence-activated cell sorting (FACS)-sorted early and late replicating fractions was analyzed using next-generation sequencing (NGS) (Fig. 1A) (Hayakawa et al., 2021; Marchal et al., 2018; Takebayashi et al., 2018). The ratio of early and late read enrichment [log₂(early-S reads/late-S reads)] was determined against the chromosomal position to delineate RT domains (Fig. 1B).

RT data were averaged into 200-kb genomic segments (with a total of 11,798 segments/dataset) and the proportion of differences greater than those observed between biological replicates was calculated, which was used identify statistically significant RT-switching regions. These proportions were converted to false discovery rate (FDR; see details in Materials and Methods) estimates to account for multiple testing. Using this approach, we identified 1615 segments as RT-switching regions during MCE1 (control versus MCE1, 1% FDR), among which 764 and 851 segments switched toward a later (early-to-late, or EtoL) and earlier (late-to-early, or LtoE) RT, respectively (Fig. 1C,D). We also identified 1145 RT-switching segments in the period up to MCE2 (control versus MCE2, 1% FDR), among which 716 and 429 segments displayed EtoL and LtoE RT switching, respectively (Fig. 1C,D). This analysis revealed that at least 20% of the genome changed its RT during MCE.

We performed a K-means clustering of these significant RT switching segments (2298 segments), revealing seven clusters with distinct patterns of RT changes in each cell state (Fig. 1E,F). Among the EtoL switching segments, those in clusters 1 and 2 tended to show relatively large RT shifts only during MCE1 and MCE2, respectively. Similarly, among the LtoE switching segments, those in clusters 5 and 6 tended to show relatively large RT shifts only during MCE1 and MCE2, respectively. These results demonstrate that each MCE cycle is required to induce specific RT changes in a set of genomic regions.

RT data obtained from E/L Repli-seq are the average of many cells; therefore, we next examined genome-wide RT changes by single-cell DNA-replication sequencing (scRepli-seq). scRepli-seq consists of sequencing amplified genomic DNA from single S phase cells and detection of the copy number difference between replicated and unreplicated regions of the genome (Fig. 2A) (Miura et al., 2020; Takahashi et al., 2019). We collected single mid-S phase cells from MCE1 and MCE2 and analyzed copy number differences to measure RT at the single cell level. We observed that RT reorganization in the population data were well recaptured in individual cells undergoing MCE (Fig. 2B,C). The existence of MCE stage-specific RT changes at the single-cell level clearly demonstrates that MCE2 is not a mere repetition of cell division, but that it is distinct from MCE1 in terms of the genome-wide RT landscape.

Early and late RT are associated with active and inactive chromatin structures, respectively (Hiratani et al., 2008; Ryba et al., 2010). EtoL and LtoE RT switching frequently coincide with transcriptional downregulation and upregulation, respectively (Hiratani et al., 2010; Rivera-Mulia et al., 2015; Takebayashi et al., 2021); therefore, we next examined the relationship between RT changes and transcriptional changes during adipogenesis. We performed transcriptome analysis in each cell state (control, MCE1 and MCE2) and identified a number of transcriptionally upregulated and downregulated genes during adipogenesis. Globally, no clear coordination between RT switching and transcriptional changes was observed (Fig. 3). This was somewhat expected because among the genes in the RT switching regions, only a subset of genes showed significant changes in transcription levels and the majority of genes remained unchanged or were not even expressed. We next focused on adipogenesis-associated genes (196 genes in the mouse gene set: HALLMARK_ADIPOGENESIS – GSEA; https://www.gsea-msigdb.org/gsea/msigdb/mouse/geneset/HALLMARK_ADIPOGENESIS.html) to examine the coordination between RT and transcription (Fig. 4A–C). As expected, these genes are upregulated during adipogenesis; in many cases, induction of low-level transcription occurred during MCE followed by remarkable upregulation after MCE (days 4–8, late adipogenesis stage). We found a clear tendency of these adipogenic genes being associated with shifts towards earlier RT. Of particular interest, RT shifts mainly occurred during MCE2 (Fig. 4B). Approximately 84% of adipocyte-associated genes were in the early-replicating regions in control cells (before differentiation), and these underwent an RT shift even earlier than control cells during MCE (Fig. 4B,D). Genes involved in fat metabolism, such as Cd36 and Fabp4, are late replicating before differentiation but switch to early replicating during MCE and before strong transcriptional activation occurs (Fig. 4E).

Close examination of the order of RT shift and transcriptional induction during MCE revealed that some adipogenic genes undergo low-level transcriptional induction prior to the RT shift, whereas others undergo low-level transcriptional induction that coincides with or follows the RT shift (Fig. 5). Thus, there is no specific order in the RT shift and transcriptional induction during MCE, whereas strong transcriptional activation generally occurs after MCE.

The observation that changes to the timing of adipogenic gene locus replication mainly occur during MCE2 (Fig. 4B) raises the possibility that the number of DNA replication cycles during MCE might be linked with cellular potential for adipogenesis. To test this, we performed double S-phase labeling with EdU and 5-chloro-2′-deoxyuridine (CldU), which distinguishes cells undergoing two rounds of replication during MCE (double positive for EdU and CldU) from those undergoing only one round of replication (only positive for either EdU or CldU) and non-replicating cells (negative for both EdU and CldU) (Fig. 6A,B; see details in the Materials and Methods). We then examined their differentiation into mature adipocytes based on the accumulation of lipid droplets (Oil Red O staining in Fig. 6B). We found that more than 80% of cells undergo one or two rounds of replication during MCE under these experimental conditions (Fig. 6C). However, the percentage of lipid droplet-positive cells was significantly higher in cells with two rounds of replication compared with those with one round of replication and those without replication (Fig. 6D). These results clearly show that the difference in the number of replication cycles during MCE affects adipocyte differentiation potential and further strengthen the significance of RT changes during MCE2.

In the present study, we carefully examined the genome-wide temporal order of DNA replication during MCE by E/L Repli-seq and scRepli-seq approaches and identified a number of RT changes. We found that (1) each round of MCE is qualitatively different in terms of the replication landscape; (2) RT switching from late to early in the second round of MCE is linked to the transcriptional induction of adipogenic genes in early adipogenesis and with full transcriptional activation in late adipogenesis; and (3) the differentiation potential of cells is linked with how many DNA replication cycles cells undergo during MCE. Thus, our study revealed a previously unrecognized dynamic RT regulation during MCE and provides a basis for understanding the biological significance of MCE.

Accumulating evidence indicates that RT is regulated in the context of cell differentiation and disease; however, understanding of the biological significance of these changes is far from complete. One can imagine a scenario in which early RT simply leads to early doubling of gene copy number during S phase and, as a result, leads to an increase in gene expression. This seems to be the case for some essential genes, such as those encoding histones and cell cycle factors (Gómez, 2017). Studies in Saccharomyces cerevisiae show that delayed RT of the histone HTA1-HTB1 locus by inactivating nearby replication origins results in failure to maintain timely and appropriate expression of these genes (Müller and Nieduszynski, 2017). S. cerevisiae has a mechanism to compensate replication-coupled gene-dosage imbalance through H3K56 acetylation by Rett109; however, some genes, such as histone and cell cycle factor genes, escape from this dosage compensation (Müller and Nieduszynski, 2017; Voichek et al., 2016). A similar dosage compensation mechanism seems to exist in mammalian cells (Padovan-Merhar et al., 2015), although the details of this phenomenon are not entirely clear. We found that many adipogenic genes undergo a shift in RT toward earlier replication in MCE (Fig. 4). It is possible that earlier replication of these genes might be advantageous for adipogenesis through earlier doubling of copy number and escape from dosage compensation, but further studies are necessary to address this point.

Different types of chromatin are assembled in different stages of S phase; DNA replicated in early and late S phase preferentially forms active and inactive chromatin, respectively, via an unknown mechanism (Zhang et al., 2002). Indeed, disruption of the RT program by inactivating RIF1 causes disruption of genome-wide histone modification patterns in a replication-dependent manner (Klein et al., 2021). Taking these observations into consideration, it is possible that LtoE switching of adipogenic genes might significantly impact chromatin structure and transcriptional activation of these genes. Our measurement of RT and transcriptome analysis at multiple stages during adipogenesis revealed that some genes undergo transcriptional activation prior to RT switching (Fig. 5). In this case, RT switching might be regarded as a consequence of transcriptional activation. However, once RT changes are induced as a downstream event of transcription, early replication could in turn be a driving force to assemble and maintain a transcriptionally competent active chromatin state via the mechanism described above. Such a positive feedback loop, regardless of which changes occur first, might play an instructive role in setting up the adipocyte-specific transcriptional program.

The genome is divided into two nuclear compartments called A/B compartments and each compartment consists of multiple topologically associating domains, as revealed by the chromosome conformation capture technique, Hi-C (Lieberman-Aiden et al., 2009). The A and B compartments are associated with active and inactive chromatin, respectively. The B compartments are characterized as preferentially interacting with the nuclear periphery, which favors transcriptionally inactive chromatin formation. This is supported by the observation that forced tethering of active genes to the nuclear periphery induces transcriptional silencing (Finlan et al., 2008; Kumaran and Spector, 2008; Reddy et al., 2008). RT profiles can be used to predict three-dimensional genome organization in the nuclear space because early and late RT domains are largely consistent with the A and B compartments (Ryba et al., 2010). Reorganization of RT domains during cell differentiation is accompanied by changes in the A/B compartments (Miura et al., 2019; Takebayashi et al., 2012). It is possible that LtoE switching of adipogenic genes is accompanied by spatial repositioning from the B to A compartment, thereby facilitating a switch from inactive to active chromatin. Replication-coupled positional changes of adipogenic genes, if occurring, might impact maintenance of the chromatin state even after cells exit from MCE (late stage of adipogenesis).

During MCE, up to two rounds of cell division occur, but little is known about exact differences between cells undergoing only one cell division and those undergoing two cell divisions. Among the RT switching regions we identified during MCE were regions undergoing gradual RT changes through two rounds of MCE, or significantly only in MCE2. Interestingly, many adipogenic genes, including those directly involved in lipid metabolism, show RT changes during MCE2 (Fig. 4B). Based on the assumption that RT is an epigenetic mark that defines a specific cell state (Hiratani and Gilbert, 2009), cells are thought to acquire different identities depending on how many MCE cycles they undergo. Indeed, the difference in the number of MCE cycles affected differentiation potential (Fig. 6).

Collectively, our study uncovered dynamic changes in the temporal order of chromosomal replication in the early step of adipogenesis, which indicates its possible involvement in adipogenic gene expression. This further strengthens the hypothesis that in adipogenesis, MCE is a window of opportunity for chromatin reorganization.

Mouse 3T3-L1 cells (Japanese Collection of Research Bioresources, lot #07232012) were cultured in Dulbecco's modified Eagle's medium (DMEM, Nacalai tesque), supplemented with 10% fetal bovine serum (FBS, lot #S05692S1820, Biowest), 2 mM L-glutamine (Nacalai tesque), 168 µM penicillin (Nacalai tesque) and 137.5 µM streptomycin (Sigma-Aldrich) in a humidified 5% CO₂ incubator at 37°C. Cells were confirmed to be negative for mycoplasma contamination. For differentiation experiments, cells were seeded in six-well plates and grown until confluency. At 2 days post confluency (defined as day 0), cells were treated with DMEM containing 10% FBS supplemented with 10 µg/ml insulin (Sigma-Aldrich), 0.25 µM dexamethasone (DEX, Sigma-Aldrich) and 0.5 mM 3-isobutyl-1methylxanthine (IBMX, Sigma-Aldrich) for 48 h to induce adipogenesis. These cells were then treated with DMEM containing 5 µg/ml insulin and the medium was replaced every 2 days up to 8 days (day 8) for full differentiation into mature adipocytes.

To monitor cell cycle progression during early stage adipogenesis, pulse-replication labeling with 100 µM EdU (Thermo Fisher Scientific) was performed every hour for the first 48 h of adipogenesis. The labeled cells were fixed with 4% paraformaldehyde for 10 min at room temperature and permeabilized with 0.5% Triton X-100 (Nacalai tesque) in 1×PBS for 15 min at room temperature. A Click-iTTM EdU Alexa Fluor 555 Imaging Kit (Thermo Fisher Scientific) was used to detect incorporated EdU following the manufacturer's protocol. Cells were stained with 200 ng/ml 4′,6-diamidino-2-phenylindole (DAPI) in 4× saline-sodium citrate buffer for 5 min before imaging. Cell images were captured using a ZEISS Axioplan 2 fluorescence microscope equipped with a ×40 NA 1.4 oil immersion objective and a CCD camera (ORCA-R2, Hamamatsu Photonics, Shizuoka, Japan) and Metamorph version 7.8.13.0 software (MDS Analytical Technologies). The percentage of S phase cells was calculated by dividing the number of EdU-positive cells by the total number of DAPI-stained cells.

Cells were labeled for 2 h with 50 µM BrdU (Sigma-Aldrich) before being fixed with 70% ethanol. We used a Sony SH800 cell sorter in ultrapurity mode to fractionate the early- and late-S phase populations based on DNA content (propidium iodide intensity). We followed our routine BrdU-immunoprecipitation (IP)-based protocol (Ryba et al., 2011; Takebayashi et al., 2018) except that we used a Bioruptor (Cosmo Bio, Tokyo, Japan) for genomic DNA sonication to obtain DNA fragments of ∼100–700 bp. After BrdU-IP, the purity of early- and late-S-phase cell fractions was validated by PCR using specific primer sets that detect early- and late-replicating regions (Ryba et al., 2011). The immunoprecipitated DNA samples were then subjected to whole genome amplification (WGA) using a SeqPlex kit (Sigma-Aldrich, SEQXE). NGS libraries were constructed from early- and late-replicating DNA after WGA using an NEBNext Ultra II DNA Library Prep Kit for Illumina (New England Biolabs) according to the manufacturer's instructions, followed by NGS on an Illumina HiSeq X system. After NGS, the raw FASTQ files were trimmed to remove adapter sequences using the Cutadapt program (https://github.com/marcelm/cutadapt) before mapping. We performed two-step adapter trimming, first removing the Illumina adapter according to the index of each NGS library and then removing the SEQXE adapter (Hayakawa et al., 2021). Mapping to the mouse mm9 reference genome (chr1–19, chrX, chrY) was performed using the Burrows–Wheeler alignment tool (https://bio-bwa.sourceforge.net/) and duplicated reads were removed using the Picard tool (http://broadinstitute.github.io/picard/). We defined MAPQ>10 as uniquely mapped reads. We also filtered out the reads that overlapped with the mm9 blacklists. After blacklist filtering, we counted the reads from the early- and late-S phase BrdU-IP samples in sliding windows of 300 kb at 40 kb intervals, and performed reads per million (rpm) normalization. The ratio of early-S phase to late-S phase read counts [log₂(early-S reads/late-S reads)] was calculated for each bin.

For some analyses, datasets were averaged into 200-kb windows (fixed position). To determine the significant RT switching regions that were independent of changes between replicates, we calculated the Euclidian distance at 11,798 200-kb genomic segments between groups (i.e. control versus MCE1) and within groups (i.e. control replicate 1 versus control replicate 2), which was used to calculate the P-values at each 200 kb genomic segment (Takebayashi et al., 2013, 2021). Statistical significance was then calculated using the q-value package in R/Bioconductor (https://cran.r-project.org/), which yields a q-value for each segment that reflects the proportion of false positives (i.e. the FDR) among segments deemed to have significant RT changes. High-confidence RT-switching domains were selected with a q-value cutoff of 0.01, corresponding to an overall FDR of 1%. Data from sex chromosomes were excluded from the analysis because the sex karyotype of 3T3-L1 cells is uncertain.

scRepli-seq was performed as described previously (Miura et al., 2020; Takahashi et al., 2019). Single mid-S-phase cells were sorted with a Sony SH800 cell sorter in single-cell mode based on DNA content (propidium iodide intensity). Single cells were sorted directly into a 96-well plate with 6 µl of cell lysis buffer [447 µl H₂O, 3 µl Proteinase K (Sigma-Aldrich, P4850), 50 µl 10× single-cell lysis and fragmentation buffer (Sigma-Aldrich, L1043)], incubated at 55°C for 1 h and then at 99°C for 4 min for genomic DNA isolation and fragmentation. Then, 6 µl of the genomic DNA solution was subjected to WGA using a SeqPlex kit (Sigma-Aldrich, SEQXE). NGS libraries of amplified genomic DNA were constructed using an NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs), according to the manufacturer's instructions, followed by NGS (Illumina HiSeq X system). After NGS, the raw FASTQ files were subjected to two-step adapter trimming, mapping to the mm9 reference genome. We filtered out duplicated reads that had an identical chromosome start position and strand information relating to an existing read. We also filtered out reads that overlapped the mm9 blacklists. We then counted the reads in sliding windows of 200 kb at 40-kb intervals and used the AneuFinder's correctMappability command (http://bioconductor.org/packages/release/bioc/html/AneuFinder.html) for normalizing mid-S-phase data based on a G1 control. From the mappability-corrected mid-S-phase read counts, the genome-wide median was obtained and was used to generate log₂[(mappability-corrected mid-S reads)/median] scores, which we defined as the single-cell mid-S-phase RT score. Positive and negative RT scores were classified as early and late RT. The binarization (replicated or unreplicated) for each genomic bin was performed using the findCNVs command in AneuFinder as described previously (Miura et al., 2020).

RNA was extracted according to the standard protocol using TRIzol (Invitrogen). The concentration of extracted total RNA and the RNA integrity number were measured using the 4200 TapeStation system (Agilent Technologies). mRNA was enriched from total RNA using the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs, E7490). RNA sequencing library preparation was performed using the NEBNext Ultra II RNA Library Prep Kit for Illumina (New England Biolabs, E7770) and fragment sizes of the obtained library were measured using the TapeStation system. NGS was performed on a Hiseq X Ten (Illumina). Index sequences added during library preparation were removed using Cutadapt software, then NGS reads were mapped to the mm9 mouse reference genome using Hisat2 (http://daehwankimlab.github.io/hisat2/). Mapped NGS read counts were used to calculate transcripts per million (TPM) values using Stringtie software (Pertea et al., 2015). The average of three biological replicates was used for the analysis.

Cells seeded on coverslips (Thermo Fisher Scientific) were induced to differentiate as described above. Cells were first labeled with 100 µM EdU for 10 min at 18 h after differentiation and subsequently labeled with 100 µM CldU (Sigma-Aldrich) during the period between 48 and 96 h after differentiation. It should be noted that cells on coverslips appeared to take longer to undergo adipogenesis, probably because of less adhesive culture conditions. After the double labeling, cells were cultured for a further 4 days. Cells were then fixed in 4% paraformaldehyde in PBS for 10 min at room temperature and then permeabilized with 0.5% Triton X-100 in PBS. For immunostaining of CldU-labeled DNA, cells were treated with 2 M HCl for 20 min at 37°C to denature cellular DNA, followed by neutralization with 0.1 M borate buffer for 15 min. Cells were then blocked with 3% Blockace (Dainippon Pharmaceutical) in PBS with 0.1% Tween 20 (PBS-T) for 30 min at 37°C, incubated with a mouse anti-BrdU monoclonal antibody (1:100; clone MoBU-1, mouse IgG, B35141; Thermo Fisher Scientific) in 1% Blockace/PBS-T for 1 h at 37°C, washed in PBS-T for 5 min, and incubated with Alexa Fluor 555-conjugated goat anti-mouse IgG (1:200; A-21425; Thermo Fisher Scientific) in 1% Blockace/PBS-T for 30 min at 37°C. After washing with PBS-T, EdU-labeled DNA was detected using Click-iT EdU Alexa Fluor 488 (Thermo Fisher Scientific) according to the manufacturer's protocol. Nuclei were stained with DAPI. To visualize lipid droplets, Oil Red O (Sigma-Aldrich) staining was performed after detection of CldU and EdU. Cells were washed with PBS and subsequently with water. These cells were incubated in 60% isopropanol for 1 min at room temperature, stained with 0.5% Oil Red O in 60% isopropanol for 30 min at room temperature, and washed again with 60% isopropanol. After washing in water, lipid droplets were visualized by light microscopy.

We thank the members of the Laboratory of Molecular and Cellular Biology at Mie University for their help and advice throughout the course of this work.

The peer review history is available online at https://journals.biologists.com/jcs/article-lookup/doi/10.1242/jcs.260778.reviewer-comments.pdf