Integrated analysis of ATAC-seq and RNA-seq reveals ADSCP2 regulates oxidative phosphorylation pathway in hypertrophic scar fibroblasts

Ethical approval and consent to participate

Human tissue samples were procured from patients who underwent surgical procedures at the Women’s Hospital of Nanjing Medical University (Nanjing Women and Children’s Healthcare Hospital). All participants furnished written informed consent. In accordance with the Declaration of Helsinki, the study was sanctioned by the Ethics Committee of the Women’s Hospital of Nanjing Medical University (Nanjing Women and Children’s Healthcare Hospital) (Approval number: 2022KY-162).

Primary hypertrophic scar fibroblasts cell culture

Utilizing a previously established protocol (Li et al., 2023), we successfully isolated fibroblasts derived from hypertrophic scars. These hypertrophic scar-derived fibroblasts were cultured from patients (n = 3) who had not undergone any prior treatment for hypertrophic scars before their surgical excision. The cell collection process involved the use of 5% Dispase II (Sigma, Burbank, CA, USA), followed by 0.2% collagenase I (Sigma, Burbank, CA, USA). The cells were then cultured in Fibroblast Medium (Catalog #2301; ScienCell, Carlsbad, CA, USA), which comprised of basal medium, a minimal concentration of fetal bovine serum (2%, FBS, Catalog #0010; ScienCell, Carlsbad, CA, USA), and fibroblast growth supplement (FGS, Catalog #2352; ScienCell, Carlsbad, CA, USA) to aid in the preservation of the fibroblastic phenotype. The cell cultivation took place at a temperature of 37 °C in an incubator with 5% CO₂ concentration. In the course of this research, fibroblasts from hypertrophic scars, specifically between passages three and seven, were utilized.

Cell treatment with ADSCP2

As previously delineated (Li et al., 2023), all peptides, with a purity exceeding 98%, were synthesized by Shanghai Science Peptide Biological Technology Co., Ltd and subsequently dissolved in distilled water (15230162; Thermo Fisher Scientific, Waltham, MA, USA). The peptide sequence for ADSCP2 was DENREKVNDQAKL, while the scrambled peptide (Sc) for ADSCP2 was NEVQADRKKELND.

ATAC-seq and sequencing data analysis

The hypertrophic scar fibroblasts were exposed to either 25 μM Sc (serving as the control) or 25 μM ADSCP2 in triplicate, after which they were consolidated into one control and one ADSCP2 sample for ATAC-seq analysis utilizing the BGI platform. ATAC-seq libraries were produced in accordance with a previously established protocol (Buenrostro et al., 2013). The process involved the isolation and purification of nuclei, which were subsequently incubated with the Tn5 transposase and supplemented with adaptors. The ATAC-seq libraries were then amplified through PCR, incorporating two distinct barcodes. Following the PCR reaction, the libraries underwent purification using Agencourt Ampure XP beads and were subsequently sequenced on a MGISEQ-2000 sequencer (BGI, Shenzhen, China) utilizing 50 bp paired-end sequencing. The raw sequencing reads were then subjected to trimming to eliminate low-quality reads and adapters, a process facilitated by fastp (version 0.22.0). The resulting clean data were aligned to the hg19 reference genome with the aid of Bowtie2 (version 2.4.2).

For the purpose of peak calling, unique reads from fragments that were less than 150 base pairs were utilized via the MACS (model-based analysis for ChIP-Seq, version: 2.1.0) algorithm. Differential peaks between samples were identified using the MAnom algorithm. The significance of each peak was determined through a Bayesian model, with significant regions being selected if the absolute value of M was greater than or equal to 1 and the p-value was less than or equal to 10⁻⁵.

Data visualization was facilitated through the use of the integrative genomics viewer (IGV 2.16.2).

RNA-seq and sequencing data analysis

In the current investigation, we used previously identified genes with the RNA-seq from hypertrophic scar fibroblasts subjected to treatment with either 25 μM Sc (control) or ADSCP2 in triplicate, followed by collection for the preparation of a transcriptome library and sequencing (Illumina NovaSeq6000, in PE150 mode), a process carried out by the ShanghaiBio Corporation (Li et al., 2023). The statistical power of this experimental design, calculated in RNASeqPower is 0.84. Here, a gene is classified as differentially expressed if it exhibits a fold change greater than 1.2 or less than 0.83, and a p-value less than 0.05 when comparing two groups. The interconnections between differential genes were scrutinized utilizing the STRING analysis tool.

Bioinformatics enrichment analysis

Differentially expressed genes (DEGs) or genes of interest were enriched using the gene ontology (GO, http://www.geneontology.org/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Based on Fisher’s exact test, the top 30 or 20 enriched GO terms or pathways were listed.

Reverse transcription and quantitative polymerase chain reaction

Hypertrophic scar fibroblasts were subjected to treatment with either 25 μM Sc (control) or ADSCP2 in triplicate. Subsequently, total RNA was extracted, reverse transcribed into complementary DNA (cDNA), and amplified utilizing the SYBR Green qPCR Master Mix (Vazyme, Nanjing, China) on an ABI Vii7 PCR Detection System (ABI, CA, USA). The expression levels of COX6B1 and NDUFA1 were normalized to 18s rRNA using the 2^(−ΔΔCt) method. Primer sequences for NDUFA1 were as follows: forward, 5′-GCG TAC ATC CAC AGG TTC ACT-3′; reverse, 5′-GCG CCT ATC TCT TTC CAT CAG A-3′. Primer sequences for COX6B1 were as follows: forward, 5′-CTA CAA GAC CGC CCC TTT TGA-3′; reverse, 5′-GCA GAG GGA CTG GTA CAC AC-3′.

Determination of ATP and lactic acid levels

Hypertrophic scar fibroblasts were exposed to either 25 μM Sc (serving as the control) or ADSCP2 in triplicate for a duration of 24 h. ATP concentrations were quantified utilizing the ATP Assay Kit (Catalog No. S0026; Beyotime, Shanghai, China) in accordance with the manufacturer’s protocol. Lactic acid concentrations were measured using the lactic acid LD detection kit (Catalog No. KGA7402-48; KeyGEN, Rockville, MD, USA). Protein concentrations were determined via the enhanced BCA protein assay kit (Catalog No. P0010; Beyotime, Shanghai, China), and ATP and lactic acid levels were normalized to protein content.

Chromatin accessibility changes following ADSCP2 treatment

The analysis of ATAC-seq data was conducted to delineate the alterations in chromatin accessibility subsequent to ADSCP2 treatment. The control and ADSCP2 treated hypertrophic scar fibroblasts for ATAC-seq were assembled into two distinct pools. Open chromatin signals were procured through alignment, duplication removal, and peak calling (Table 1).

Using MACS (Zhang et al., 2008), we discerned 19,079 accessible chromatin regions (ACRs) in the control sample and 18,061 ACRs in the ADSCP2 sample (Fig. 1A). These ACRs were predominantly concentrated in the intergenic and intron regions in both control and ADSCP2 treated hypertrophic scar fibroblasts (Fig. 1B). In the control group, 50.1% and 44% of ACRs positioned at the intergenic and intron regions respectively (Fig. 1B). In the ADSCP2 group, 49.5% and 44.2% of ACRs positioned at the intergenic and intron regions respectively (Fig. 1B).

Utilizing the MAnorm method (Shao et al., 2012), a total of 7,805 distinct ACRs were identified (Fig. 1A). The differential peaks observed between the control and ADSCP2 were primarily situated in the intergenic and intron regions (Fig. 1C). A comprehensive view of the open regions across all chromosomes was presented for the control and ADSCP2 group (Fig. 1D). A total of 3,176 genes were associated with 7,805 differential peaks, and genes exhibiting more than five differential peaks were cataloged (Fig. 1E). GO analysis indicated that the 3,176 genes associated with the 7,805 differential peaks were predominantly related to cellular processes, biological regulation, cell structure, organelle function, binding, and catalytic activity (Fig. 1F). Furthermore, the KEGG pathway analysis revealed that these differential peak genes were primarily linked with focal adhesion, calcium signaling pathway, and regulation of actin cytoskeleton (Fig. 1E).

Analysis of genes at ACRs following ADSCP2 treatment

Following the treatment with ADSCP2, the Venn analysis revealed an increase in accessibility for 640 peaks, while 1,658 peaks demonstrated reduced accessibility (Fig. 2A). Genes associated with close-ACRs were predominantly linked to the response to xenobiotic stimuli, embryonic organ development, and class II regionalization (Fig. 2B). The KEGG pathway analysis indicated that the close-ACRs genes were primarily enriched in oxidative phosphorylation, ribosome, and neuroactive ligand-receptor interaction (Fig. 2B). Conversely, genes associated with open-ACRs were predominantly linked to the development of the reproductive system, reproductive structure, and ribosome biogenesis (Fig. 2B). Analyses of KEGG pathways revealed that cytokine-cytokine receptor interaction, ribosome, and JAK-STAT signaling pathway were predominantly enriched for the open-ACRs genes, as depicted in Fig. 2B.

Integrated analysis of ATAC-seq and RNA-seq

In a previous study, we examined the RNA-seq data of ADSCP2, utilizing a parameter of fold change ≥2, P < 0.05, which led to the identification of 160 upregulated genes and 170 downregulated genes (Li et al., 2023). In the current investigation, we used previously identified genes, modified the parameter to fold change ≥1.2, P < 0.05, with the aim of detecting the minimal changes associated with the functions of ADSCP2. This adjustment resulted in the identification of 345 upregulated genes and 399 downregulated genes, as shown in Fig. 3A.

The study presents the top 30 upregulated and downregulated genes (Fig. 3B). An analysis of the top seven enriched GO terms, encompassing biological process, cellular component, and molecular function, revealed a predominant enrichment in cellular response to arsenic-containing substance, cytoplasmic translation, cytosolic ribosome, proteasome binding, and structural constituent of ribosome (Fig. 3C). Furthermore, an examination of the top seven enriched KEGG pathways indicated a significant enrichment in ribosome and oxidative phosphorylation (Fig. 3C). A protein-protein interaction analysis, conducted using the search tool for the retrieval of interacting genes/proteins (STRING), elucidated the intricate interaction among these genes.

An integrated analysis of ATAC-seq and RNA-seq, utilizing venn analysis, revealed that in the ADSCP2 treatment group, six genes were upregulated and exhibited increased accessibility, while 27 genes were downregulated and demonstrated decreased accessibility (Fig. 4A). The differential peaks coincided with 19 downregulated genes and 26 upregulated genes (Fig. 4A). The KEGG pathway analysis indicated that the downregulated genes were primarily associated with oxidative phosphorylation (Fig. 4B), whereas the upregulated genes were predominantly related to the spliceosome (Fig. 4C). The joint assessment of ATAC-seq and RNA-seq data identified two oxidative phosphorylation-related genes, COX6B1 (cytochrome c oxidase subunit 6B1) and NDUFA1 (NADH dehydrogenase (ubiquinone) alpha subcomplex-1). The integrative genomics viewer and RNA-seq data of COX6B1 and NDUFA1 were further analyzed (Fig. 4D).

As depicted in Fig. 4D, the control treatment group exhibited three peaks for COX6B1 and two peaks for NDUFA1, while the ADSCP2 treatment group showed one peak for COX6B1 and none for NDUFA1. This suggests reduced accessibility of COX6B1 and NDUFA1 following ADSCP2 treatment. A subsequent analysis of the RNA-seq data revealed a downregulation of both COX6B1 and NDUFA1 in the ADSCP2 group. According to the UCSC genome browser analysis, the promoter regions of COX6B1 and NDUFA1 are characterized by a high prevalence of histone modifications (specifically, histone H3 lysine 27 acetylation, or H3K27Ac) and multiple transcription factor binding sites (Figs. 5, 6).

Peak 1 is situated at the promoter region, while peaks 2 and 3 are found at the intron regions of COX6B1 (Fig. 5B). Similarly, for NDUFA1, peak 1 is located at the promoter region, and peak 2 is found at the intron region (Fig. 6B). The JASPAR transcription factor prediction analysis suggests that transcription factors such as ELK4 (ETS transcription factor), FLI1 (friend leukemia virus integration 1), ZNFs (zinc finger protein), and KLFs (Kruppel-like factors) may potentially bind with the COX6B1 promoter region, thereby influencing its transcription (Fig. 5B). Additionally, SP1 (specificity protein 1), KLFs, ZNFs, and SREBFs (sterol regulatory-element binding proteins) are likely to bind with NDUFA1 promoter regions (Fig. 6B). These findings show that ADSCP2 reduces the transcriptions of COX6B1 and NDUFA1 genes useful for the OXPHOS process. The mechanisms may be that ADSCP2 functions by regulating histone modifications or transcription factor binding sites on those promoter regions.

ADSCP2 inhibits the mRNA expression levels of COX6B1 and NDUFA1, as well as reduces cellular ATP and lactic acid levels

To investigate the impact of ADSCP2 on the transcription of COX6B1 and NDUFA1, we conducted quantitative PCR analysis. The results indicated a reduction in mRNA levels of both COX6B1 and NDUFA1 in the ADSCP2 treatment group (Fig. 7A). Additionally, we assessed ATP and lactic acid levels in cell pellets and culture medium following ADSCP2 treatment (Fig. 7B). The findings revealed a decrease in intracellular ATP levels, while ATP levels in the culture medium were elevated, albeit marginally, post-ADSCP2 treatment (Fig. 7C). Furthermore, intracellular lactic acid levels were significantly reduced, whereas lactic acid levels in the culture medium remained unchanged after ADSCP2 treatment (Fig. 7D).