Abstract
E/L Repli-seq is a powerful tool for detecting cell type-specific replication landscapes in mammalian cells, but its potential to monitor DNA replication under replication stress awaits better understanding. Here, we used E/L Repli-seq to examine the temporal order of DNA replication in human retinal pigment epithelium cells treated with the topoisomerase I inhibitor camptothecin. We found that the replication profiles by E/L Repli-seq exhibit characteristic patterns after replication-stress induction, including the loss of specific initiation zones within individual early replication timing domains. We also observed global disappearance of the replication timing domain structures in the profiles, which can be explained by checkpoint-dependent suppression of replication initiation. Thus, our results demonstrate the effectiveness of E/L Repli-seq at identifying cells with replication-stress-induced altered DNA replication programs.
Introduction
In mammalian cells, DNA replication is regulated in Mb-sized chromosomal units in which clustered replication origins are synchronously activated [Rhind and Gilbert, 2013; Rivera-Mulia and Gilbert, 2016a; Takebayashi et al., 2017]. Each unit is defined as an early or late replication timing domain according to replication timing within the S phase. Genome-wide studies to map replication timing domains have revealed that chromosomes consist of mosaic patterns of early and late replication timing domains, the patterns of which are highly cell type-specific [Hiratani et al., 2008, 2010; Rivera-Mulia et al., 2015, 2018; Dixon et al., 2018]. For instance, the domain organization in embryonic stem (ES) cells is almost indistinguishable in several independent ES cell lines but is clearly distinct from that in differentiated cell types [Hiratani et al., 2008; Ryba et al., 2010]. Because early and late replication timing domains are generally associated with active and inactive chromatin modifications [Ryba et al., 2010; Rhind and Gilbert, 2013; Dileep et al., 2015b], reorganization of the replication timing domain, which occurs during cell differentiation, is thought to reflect changes in transcriptional competence.
E/L Repli-seq is one of the methodologies that have been used to map genome-wide replication timing domains. The method employs immunoprecipitation of 5-bromo-2′-deoxyuridine (BrdU)-labeled replicating DNA from FACS-sorted early and late S-phase cells and determines quantitative the enrichment ratio between early and late replicating sequences in each genomic region using next-generation sequencing (NGS) technologies [Marchal et al., 2018; Takebayashi et al., 2018]. Several studies have shown that the distribution of early and late replication timing domains along the chromosome is largely in line with those of A (active) and B (inactive) nuclear compartments based on chromatin interaction patterns [Ryba et al., 2010; Pope et al., 2012; Dileep et al., 2015a]. This suggests that replication timing domains are fundamental structural units on mammalian chromosomes, and also that E/L Repli-seq has the potential to predict genome organization. Other applications for E/L Repli-seq are possible and much remains to be examined.
In addition to developmental cues, replication stress can influence the replication landscape in mammalian cells. When cells experience replication stress, significant changes in origin firing regulation occur. It has been shown that DNA damage during S phase induces 2 major (seemingly paradoxical) changes in origin firing regulation [Ge and Blow, 2010; Blow et al., 2011; Yekezare et al., 2013]. First, replication fork slowing induced by replication stress activates dormant replication origin firing (hyperactivation of origins), which is thought to be a compensatory mechanism for fork slowing to allow the completion of replication within already actively replicating genomic regions. Second, origin activation is suppressed in a checkpoint-dependent manner, which is thought to prevent cells from further damage and/or allows cells to repair damaged sites. These observations, however, have been made through single-molecule analysis of origin firing events on extended DNA fibers and visualization of spatiotemporal regulation of replication factories in the nucleus [Dimitrova and Gilbert, 2000; Feijoo et al., 2001; Ge et al., 2007]. It remains unclear how such complex replication regulation is linked to the domain-level regulation detected by E/L Repli-seq. Therefore, the rationale for this study was to directly examine the effect of dormant origin activation and checkpoint-dependent origin suppression on E/L Repli-seq profiling using camptothecin (CPT)-treated human retinal pigment epithelium (RPE) cells as a model.
Materials and Methods
Cell Culture and Drug Treatment
Human RPE cells were grown in MEMα (Wako, Osaka, Japan) supplemented with 10% FBS at 37°C with 5% CO2. DNA replication stress was induced by exposure to 1 μM CPT (TopoGEN, Buena Vista, VA, USA) with or without 1 mM caffeine (Nacalai Tesque, Kyoto, Japan).
Immunofluorescence Staining of γH2AX
RPE cells were fixed with 4% paraformaldehyde (Nacalai Tesque)/phosphate-buffered saline (PBS) for 10 min. After fixation, samples were washed 3 times with PBS and permeabilized with 0.5% Triton X-100 (Nacalai Tesque)/PBS for 15 min at room temperature (RT). After washing with PBS, blocking was done in 3% Blockace (Yukijirushi, Hokkaido, Japan)/PBS for 1 h at 37°C. Anti-phospho-histone H2A.X (Ser139) (1:500; Merck, Darmstadt, Germany) was the primary antibody, followed by Alexa Fluor® 488 donkey anti-mouse IgG (H + L) secondary antibody (1:500; Invitrogen, Waltham, MA, USA, cat. #21202). Both primary and secondary antibody reactions were performed at 37°C for 1 h.
Replication Labeling with EdU
Cells were labeled for 2 h with 50 μM 5-ethynyl-2′-deoxyuridine (EdU; Thermo Fisher Scientific, Waltham, MA, USA) at 37°C with 5% CO2. The labeled cells were fixed with 4% paraformaldehyde (Nacalai Tesque)/PBS for 10 min and permeabilized with 0.5% Triton X-100 (Nacalai Tesque)/PBS for 15 min at RT. After washing with PBS, Click-iTTM EdU Alexa Fluor® 555 Imaging Kit (Thermo Fisher Scientific) was used to detect EdU following the manufacturer’s protocol. After EdU detection, nuclei were stained with 200 ng/mL 4′,6-diamidino-2-phenylindole (DAPI; Roche, Indianapolis, IN, USA)/4× saline-sodium citrate (SSC) for 3 min and 3 PBS washes before mounting in Vectashield (Vector Laboratories, Burlingame, CA, USA).
Molecular Combing and Immunofluorescent Detection
Cultured cells were labeled for 20 min with 100 μM 5-iodo-2′-deoxyuridine (IdU; Sigma-Aldrich, St. Louis, MO, USA) in culture medium. After washing with PBS, cells were labeled for 20 min with 100 μM 5-chloro-2′-deoxyuridine (CldU; Sigma-Aldrich). DNA replication labeling was stopped with PBS at 4°C, and cells were centrifuged (1,500 rpm, 4°C, 3 min). After removing supernatant, buffer A (250 mM sucrose, 20 mM HEPES pH 7.0, 10 mM KCl, 1.5 mM MgCl2, 1 mM EDTA pH 8.0, 1 mM EGTA) was used to collect genomic DNA. Then, 1.4% agarose LM-MP (Roche) were mixed, and the solution was run into plug mold. After 2 h on ice, the gel was incubated in ESP solution (10 mM N-lauroylsarcosine sodium salt, 10 mM Tris-HCl pH 8.0, 0.5 M EDTA pH 8.0, and 100 μM Proteinase K) for 16 h to break down proteins. Next, the gel was washed with 0.5 M EDTA for 3 times and twice in 1× TE (10 mM Tris-HCl, 1 mM EDTA pH 8.0). To stain genomic DNA, T40E2 (40 mM Tris-HCl, 2 mM EDTA pH 8.0) with 4 μM YOYO-1 (Life Technologies, Carlsbad, CA, USA) was added. After 1 h at RT, the gel was washed with 1× TE. Then, the gel was melted at 70°C, and β-agarase I (NIPPON GENE, Tokyo, Japan) was added to complete gel digestion at 65°C. Finally, the DNA solution was diluted with 0.5 M MES and kept at 4°C. Genomic DNA was prepared and combed onto the silanized coverslips (Matsunami Glass, Osaka, Japan). Combed DNA molecules were denatured in 50% formamide (Roche) and 2× SSC at 72°C for 12 min. For immunodetection of IdU- and CldU-labeled DNA, denatured DNA molecules were incubated with mouse anti-BrdU monoclonal antibody (1:5; Becton Dickinson, San Jose, CA, USA) and rat anti-BrdU monoclonal antibody (1:50; Abcam, Cambridge, UK) for 1 h at 37°C. After washing with 0.05% Tween 20/PBS, coverslips were incubated with Alexa Fluor 555-conjugated goat anti-mouse IgG (1:100; Life Technologies) and Alexa Fluor 488-conjugated rabbit anti-rat IgG (1:100; Life Technologies) for 45 min at 37°C. After washing with 0.05% Tween 20/PBS, coverslips were mounted in Vectashield.
Statistical Analysis
Statistical analyses of each CldU replication fork rate and ori-to-ori distance were performed by Mann-Whitney U test. p values <0.05 were considered statistically significant.
Imaging
DNA fibers and cell images were captured by ZEISS Axioplan 2 imaging fluorescence microscopy (×40 with an NA 1.4 oil immersion objective equipped with a CCD camera; ORCAR2, Hamamatsu Photonics, Shizuoka, Japan) and Meta Morph version 7.8.13.0 software (MDS Analytical Technologies). DNA molecule extension was estimated on coverslips prepared with λ-DNA, and the DNA molecules were stained with YOYO-1 (Invitrogen) for 1 h at RT. The DNA molecules measured 19 ± 1.2 μm. Because the length of the λ genome is 48.5 kb, the extension of DNA molecules is 2.55 ± 0.16 kb/μm.
Genome-wide Replication Timing Analysis by E/L Repli-Seq
We followed our routine BrdU-immunoprecipitation (IP)-based protocol [Takebayashi et al., 2013; Takahashi et al., 2019] except that we used a Bioruptor (Cosmo Bio, Tokyo, Japan) for gDNA sonication to obtain DNA fragments of 100–700 bp. We used a Sony SH800 cell sorter in ultrapurity mode to fractionate the early and late S-phase populations based on DNA content (PI intensity). Considering that the length of S phase is about 10 h in exponentially growing RPE cells and the duration of BrdU treatment is 2 h, the percentage of BrdU-labeled DNA in each sorted early or late S-phase fraction is estimated to be up to 20% of total DNA. After BrdU-IP, the purity of early and late S-phase cell fractions was validated by PCR using specific primer sets that detect control early and late replicating regions [Ryba et al., 2011]. These immunoprecipitated DNA samples were subjected to whole-genome amplification (WGA) with a SeqPlex kit (Sigma, SEQXE). NGS libraries were constructed from early and late replicating DNA after WGA with a TruSeq DNA PCR-Free High Throughput Library Prep Kit (Illumina, San Diego, CA, USA) and NEXTflex® DNA-Seq Kit, both according to the manufacturers’ instructions, followed by NGS (Illumina Hiseq X system). After NGS, the raw Fastq files were trimmed to remove adapter sequences using the cutadapt program before mapping. We performed 2-step adapter trimming, first removing the Illumina adapter according to the index of each NGS library and then removing the SEQXE adapter [Martin, 2011; Miura et al., 2020]. hg19 reference genome (chr1–22, chrX, chrY) assemblies were used. Mapping to the hg19 reference genome was done using the BWA program. We used the Picard tool to remove duplicated reads and defining MAPQ >10 as uniquely mapped reads. We also filtered out the reads that overlapped with the hg19 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 normalization. The ratio of early S-phase to late S-phase read counts log2 (Early S reads/Late S reads) was calculated for each bin.
Results and Discussion
To evaluate the extent to which replication stress-induced fork slowing and subsequent altered origin firing affect the replication profiles detected by E/L Repli-seq, we used human RPE cells exposed to 1 μM CPT, a potent topoisomerase I inhibitor. At this concentration, studies have shown that CPT induces DNA damage at replication forks and associated fork slowing [Sugimura et al., 2008; Regairaz et al., 2011]. Consistent with this, molecular combing analysis of replicating DNA showed that replication fork progression slowed in the presence of 1 μM CPT (Fig. 1a). This was accompanied by decreased ori-to-ori distances (Fig. 1b) via dormant origin activation, confirming that the compensatory mechanism for fork slowing is maintained in this cell line.
Replication stress induces characteristic changes in the E/L Repli-seq profile. a CPT-induced replication fork slowing revealed by DNA molecular combing. A representative image of replication fork progression visualized by consecutive pulse labeling with 100 μM IdU for 20 min (red) and CldU for 20 min (green). Labeled DNA combed on the coverslip was detected by immunofluorescence staining (see Materials and Methods). Scale bar, 10 kb. The histogram shows the measured replication fork speed (kb/min) in the control and in 1 μM CPT-treated cells (control n = 257, CPT-treated n = 195). p < 0.01 versus control. Replication labeling with IdU and CldU was performed in the presence or absence of CPT. b CPT-induced dormant origin firing revealed by DNA molecular combing. A representative image of replication origin firing visualized by the same method as in a. Scale bar, 10 kb. The histogram shows the measured ori-to-ori distances (kb) in the control and in 1 μM CPT-treated cells (control n = 128, CPT-treated n = 99). p < 0.01 versus control. c Experimental design for inducing different degrees of replication stress. d CPT-induced replication stress resulted in DNA damage and checkpoint-dependent replication inhibition. Human RPE cells were replication-labeled with 50 μM EdU for 30 min without (control) or with 1 μM CPT. γH2AX immunofluorescence staining (green) was performed simultaneously with the click chemistry-mediated detection of EdU (red). Nuclei were stained with DAPI (blue). Caffeine was also added during CPT treatment in some cases to see the effect of checkpoint inhibition. Checkpoint-dependent replication inhibition is more prominent in the longer (7 h) CPT treatment. Arrowheads indicate that γH2AX-positive cells are also positive for EdU. Scale bar, 50 μm. e Experimental overview of E/L Repli-seq analysis. Cells treated as in c were pulse-labeled with 50 μM BrdU for 2 h. Early and late S-phase cells were flow cytometrically fractionated according to their DNA content (PI-stained). BrdU-substituted DNA preparations from early and late S-phase cells were separately subjected to NGS library construction and sequencing. f Replication timing profiling by E/L Repli-seq. Mapped NGS reads from early and late S-phase cells, which were counted separately in sliding windows of 300 kb at 40-kb intervals, were used to generate a log2 (Early S reads/Late S reads) plot. Early and late replication timing domains are colored red and blue, respectively. Shown are E/L Repli-seq profiles for 2 example chromosomes (chr 9 and 12) obtained from different treatment conditions. Expanded plots for the rectangle areas are shown on the right (green horizontal lines indicate preferred initiation zones within early replicating domains of control cells). The gap seen in the chr 9 plot is an unmappable region caused by repetitive sequences. Each plot represents the average of 2 independent experimental replicates.
Replication stress induces characteristic changes in the E/L Repli-seq profile. a CPT-induced replication fork slowing revealed by DNA molecular combing. A representative image of replication fork progression visualized by consecutive pulse labeling with 100 μM IdU for 20 min (red) and CldU for 20 min (green). Labeled DNA combed on the coverslip was detected by immunofluorescence staining (see Materials and Methods). Scale bar, 10 kb. The histogram shows the measured replication fork speed (kb/min) in the control and in 1 μM CPT-treated cells (control n = 257, CPT-treated n = 195). p < 0.01 versus control. Replication labeling with IdU and CldU was performed in the presence or absence of CPT. b CPT-induced dormant origin firing revealed by DNA molecular combing. A representative image of replication origin firing visualized by the same method as in a. Scale bar, 10 kb. The histogram shows the measured ori-to-ori distances (kb) in the control and in 1 μM CPT-treated cells (control n = 128, CPT-treated n = 99). p < 0.01 versus control. c Experimental design for inducing different degrees of replication stress. d CPT-induced replication stress resulted in DNA damage and checkpoint-dependent replication inhibition. Human RPE cells were replication-labeled with 50 μM EdU for 30 min without (control) or with 1 μM CPT. γH2AX immunofluorescence staining (green) was performed simultaneously with the click chemistry-mediated detection of EdU (red). Nuclei were stained with DAPI (blue). Caffeine was also added during CPT treatment in some cases to see the effect of checkpoint inhibition. Checkpoint-dependent replication inhibition is more prominent in the longer (7 h) CPT treatment. Arrowheads indicate that γH2AX-positive cells are also positive for EdU. Scale bar, 50 μm. e Experimental overview of E/L Repli-seq analysis. Cells treated as in c were pulse-labeled with 50 μM BrdU for 2 h. Early and late S-phase cells were flow cytometrically fractionated according to their DNA content (PI-stained). BrdU-substituted DNA preparations from early and late S-phase cells were separately subjected to NGS library construction and sequencing. f Replication timing profiling by E/L Repli-seq. Mapped NGS reads from early and late S-phase cells, which were counted separately in sliding windows of 300 kb at 40-kb intervals, were used to generate a log2 (Early S reads/Late S reads) plot. Early and late replication timing domains are colored red and blue, respectively. Shown are E/L Repli-seq profiles for 2 example chromosomes (chr 9 and 12) obtained from different treatment conditions. Expanded plots for the rectangle areas are shown on the right (green horizontal lines indicate preferred initiation zones within early replicating domains of control cells). The gap seen in the chr 9 plot is an unmappable region caused by repetitive sequences. Each plot represents the average of 2 independent experimental replicates.
We also evaluated the responses to different replication stress levels by treating RPE cells with 1 μM CPT for 2 or 7 h (Fig. 1c). Cells treated with 1 μM CPT for 2 h stained strongly with an antibody against histone H2AX phosphorylation (γH2AX), and the γH2AX-positive cells showed EdU incorporation (Fig. 1c, d), a result indicative of replication-dependent DNA double-strand breaks. γH2AX accumulation was also observed in cells treated with CPT for 7 h, but in these cells EdU incorporation was significantly reduced via checkpoint-dependent origin suppression. The severe inhibitory effect on EdU incorporation was canceled in the presence of the checkpoint inhibitor caffeine (Fig. 1c, d), confirming that checkpoint-dependent origin inhibition is strongly induced by longer CPT treatment and caffeine abrogates such origin inhibition possibly through interfering ATR and CHK1 kinase signaling [Feijoo et al., 2001].
Given that the above replication stress responses are normally maintained in RPE cells, we then analyzed these cells under different treatment conditions by E/L Repli-seq (Fig. 1e, f; online suppl. Fig. 1; for all online suppl. material, see www.karger.com/doi/10.1159/000518263). The distribution patterns of the early and late replication domains along the chromosomes were globally comparable between the control and the cells treated with CPT for 2 h. However, a closer look at the individual early domains revealed damage-induced specific changes in the E/L Repli-seq profiles. In the control cells, we observed several distinct peaks in the profiles (Fig. 1f, green horizontal lines in the right panels). These peaks were highly concordant between experimental replicates (online suppl. Fig. 2) and may correspond to regions where replication initiation occurs specifically in early S phase. However, many of these peaks were not detectable after CPT treatment for 2 h, and relatively flat profiles within early domains were seen. This phenomenon is thought to be checkpoint-independent because these damage-induced changes were not alleviated in the presence of caffeine.
On the other hand, no clear domain structure was observed in the E/L Repli-seq profiles from cells treated with CPT for 7 h (Fig. 1f). This is probably caused by the much greater accumulation of DNA damage incurred by the longer treatment strongly suppressing origin firing through checkpoint activation, and the resultant inhibition of BrdU incorporation made it technically difficult to immunoprecipitate the genomic fragments replicating in early and late S phases. Simultaneous caffeine addition restored the early/late domain organization to levels comparable to the control cells (Fig. 1f), which is explained by caffeine-mediated abrogation of origin suppression and the resultant incorporation of nucleoside analogues. Thus, these results clearly demonstrate that the observed phenomenon occurs in a checkpoint-dependent manner. It should be noted, however, that simultaneous caffeine addition did not restore E/L Repli-seq peak distribution within the individual early domains. It has been proposed that, when cells encounter problems that threaten fork progression, dormant replication origin activation occurs within the already active early replication timing domains, and new origin firing within late domains is suppressed in response to the DNA damage [Rivera-Mulia and Gilbert, 2016b]. This model explains why domain-level regulation of DNA replication is important for maintaining genome integrity, but it is not consistent with our observation that checkpoint-dependent origin suppression is globally induced, regardless of whether origins are in early or late replication domains.
Replication stress-induced fork slowing activates dormant origins as a compensatory mechanism [Ge and Blow, 2010; Blow et al., 2011; Yekezare et al., 2013]. We therefore tested whether dormant origin activation is involved in E/L Repli-seq profile changes. Caffeine treatment during unperturbed S phase induces dormant origin activation [Maya-Mendoza et al., 2007], and we observed that the ori-to-ori distance decreased in caffeine-exposed cells (Fig. 2a). This clearly shows that dormant origin activation does indeed occur without any DNA damage-inducing treatment. E/L Repli-seq analysis revealed that the distribution of early/late domains is indistinguishable between the caffeine-treated cells and control (Fig. 2b). Notably, distinct E/L Repli-seq peaks within the domains lost during 2-h CPT treatment were highly conserved in caffeine-treated cells. This suggests that the replication timing changes seen within individual early domains (delineated by flattened patterns in the E/L Repli-seq profile) are not caused by dormant origin activation and may reflect previously unrecognized replication regulation in response to replication stress. Consistent with our observation, dormant origin activation seen in aging cells has been shown to have little impact on the E/L Repli-seq profile [Rivera-Mulia et al., 2018]. Currently, we do not have enough data to discuss the mechanism and biological significance of this altered replication regulation, and therefore further studies are required. For example, in order to examine the effect of DNA repair on this phenomenon, it would be necessary to perform experiments using repair-deficient mutants such as WRN patient cells [Croteau et al., 2014].
Dormant origin activation does not explain changes in the E/L Repli-seq profile. a Dormant origin activation induced solely by caffeine treatment. RPE cells were cultured in the presence of 1 mM caffeine for 2 h and the ori-to-ori distance was measured as described in Figure 1b (control n = 128, caffeine-treated n = 112). p < 0.01 versus control. b E/L Repli-seq profiles obtained from caffeine-treated cells. RPE cells were cultured with 1 mM caffeine for up to 24 h and were subjected to E/L Repli-seq.
Dormant origin activation does not explain changes in the E/L Repli-seq profile. a Dormant origin activation induced solely by caffeine treatment. RPE cells were cultured in the presence of 1 mM caffeine for 2 h and the ori-to-ori distance was measured as described in Figure 1b (control n = 128, caffeine-treated n = 112). p < 0.01 versus control. b E/L Repli-seq profiles obtained from caffeine-treated cells. RPE cells were cultured with 1 mM caffeine for up to 24 h and were subjected to E/L Repli-seq.
Collectively, our study shows that E/L Repli-seq can be used to detect genome-wide alterations in the DNA replication programs of cells under different replication stress levels. While the replication stress induced by 2-h CPT treatment resulted in the loss of preferred initiation zones within individual early domains, the more severe replication stress induced by 7-h CPT treatment caused a failure to detect clear replication timing domain organization via checkpoint-mediated replication inhibition. In more than 70 cell types from mouse and human, replication timing under normal (unperturbed) S-phase conditions has been profiled using E/L Repli-chip and E/L Repli-seq [Weddington et al., 2008; Yue et al., 2014]. However, Repli-chip and Repli-seq data from cells undergoing severe replication stress have been largely lacking. Our study is the first to show the impact of severe replication stress on Repli-seq profiling. Recently, E/L Repli-seq was used to identify a subset of genomic regions undergoing replication timing changes in cells treated with a sufficiently low enough concentration of aphidicolin to induce common fragile sites without affecting S-phase progression [Sarni et al., 2020]. By accumulating more E/L Repli-seq data from cells undergoing various forms of replication stress, it may become feasible to predict the type and degree of replication stress that cells are exposed to based on E/L Repli-seq profiles. Although more laborious than E/L Repli-seq, the use of high-resolution Repli-seq, that is sensitive enough to detect replication initiation, elongation, and termination events, may facilitate the precise assessment of replication stress [Zhao et al., 2020].
Acknowledgments
We thank Drs. S. Takahashi and I. Hiratani for their help in cell sorting.
Statement of Ethics
The immortalized cell line used in this study was purchased from ATCC. Ethical approval for the use of this cell line is not required in accordance with local/national guidelines.
Conflict of Interest Statement
The authors have no conflicts of interest directly relevant to the content of this article.
Funding Sources
This work was supported by a Grant-in-Aid for Scientific Research on Innovative Areas (grant number 20H05386 to S.T.) from the Ministry of Education, Culture, Sports, Science and Technology of Japan and a Grant from the Japan Science and Technology Agency (CREST) to S.T.
Author Contributions
S.T. designed the research; S.T. and T.H. carried out experiments, analyzed data, and wrote the paper; R.S., K.K., and K.O. supported the design and execution of the project.
Data Availability Statement
Next-generation sequencing data obtained in this study were deposited in the Gene Expression Omnibus (GEO) database under the accession number GSE167045.