Chromatin accessibility analysis reveals that TFAP2A promotes angiogenesis in acquired resistance to anlotinib in lung cancer cells.

Anlotinib, a multitarget tyrosine kinase inhibitor, is effective as a third-line treatment against non-small cell lung cancer (NSCLC). However, acquired resistance occurs during its administration. To understand the molecular mechanisms of anlotinib resistance, we characterized chromatin accessibility in both the parental and anlotinib-resistant lung cancer cell line NCI-H1975 through ATAC-seq. Compared with the parental cells, we identified 2666 genomic regions with greater accessibility in anlotinib-resistant cells, in which angiogenesis-related processes and the motifs of 21 transcription factors were enriched. Among these transcription factors, TFAP2A was upregulated. TFAP2A knockdown robustly diminished tumor-induced angiogenesis and partially rescued the anti-angiogenic activity of anlotinib. Furthermore, transcriptome analysis indicated that 2280 genes were downregulated in anlotinib-resistant cells with TFAP2A knocked down, among which the PDGFR, TGF-β, and VEGFR signaling pathways were enriched. Meanwhile, we demonstrated that TFAP2A binds to accessible sites within BMP4 and HSPG2. Collectively, this study suggests that TFAP2A accelerates anlotinib resistance by promoting tumor-induced angiogenesis.

Introduction

Lung cancer is the leading cause of cancer-related mortality worldwide [1]. Among all kinds of lung cancers, non-small cell lung cancer (NSCLC) accounts for more than 80% of lung cancer cases [2, 3]. Most patients have reached an advanced stage at the time of diagnosis and have lost the opportunity for surgery [2]. Anlotinib is a multitarget agent that inhibits tumor proliferation and angiogenesis. In ALTER randomized clinical trials (ALTER 0302 and ALTER 0303), we found that anlotinib significantly prolonged progression-free survival (PFS) and overall survival (OS) in NSCLC patients at third line [4, 5]. However, acquired resistance to anlotinib in NSCLC is inevitable. Thus, understanding the molecular mechanisms of anlotinib resistance is necessary for improving patient clinical outcomes.

Epigenetic aberrations play an important role in tumor initiation and progression [6]. The high rate of epigenetic changes in tumors results in altered gene expression patterns, which leads to the rapid evolution of tumors through the selection of drugs [7]. Epigenetic events, such as DNA methylation at the promoter region, have been found to initiate acquired resistance to antitumor drugs [8-10]. Chromatin accessibility is an important epigenetic event and recent studies have identified chromatin state changes in tumor initiation, migration, tumor metastatic progression, and drug resistance [6, 11–13].

Alterations in chromatin accessibility affect the binding of transcription factors (TFs) to their cognate genomic sequences. Many studies of tumor drug resistance have identified a role for TFs, such as MYC [14, 15], C/EBPβ [16], YAP [17], and the TEAD family of factors [18]. In addition, TFs also serve as key regulators in tumor-induced angiogenesis. In colorectal cancer, FOX1 can transcriptionally activate VEGFA to accelerate tumor-induced angiogenesis [19], and in breast tumors SOX4 directly regulates endothelin-1 expression to promote tumor-induced angiogenesis [20].

Here we performed assay for transposase-accessible chromatin using sequencing (ATAC-seq) [21] in both the parental and anlotinib-resistant (AR) lung cancer cell line NCI-H1975 to investigate chromatin accessibility changes during the process of anlotinib resistance. Compared with the parental cells, we identified 2666 regions of greater transcription accessible in AR cells, where angiogenesis-related processes and the motifs of 21 TFs were enriched. We demonstrated that TFAP2A promotes angiogenesis and accelerates anlotinib resistance.

Materials and methods

Cell culture

NCI-H1975 cells were cultured in RPMI-1640 (Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS) (Invitrogen, Carlsbad, CA, USA) and 1% penicillin/streptomycin (Invitrogen, Carlsbad, CA, USA). All cells were cultured in a 37 °C incubator with 5% CO2 and mycoplasma contamination was checked using a TransDetect PCR Mycoplasma Detection Kit (TransGen, Beijing, China). Human microvascular endothelial cells (HMEC-1) were purchased from Shanghai Zhong Qiao Xin Zhou Biotechnology Co. Ltd (Shanghai, China). HMEC-1 cells were cultured in endothelial cell medium (ScienCell, San Diego, CA, USA) supplemented with 5% FBS, 1% endothelial cell growth supplement (ScienCell, San Diego, CA, USA), and 1% penicillin/streptomycin.

RNA-seq library preparation and data analysis

mRNA was purified using NEBNext Poly (A) mRNA Magnetic Isolation Beads (NEB, USA). Library preparation of each sample was performed using the NEBNext Ultra Directional RNA Library PrepKit (NEB, USA) and samples were sequenced at 150 bp paired-end on HiSeq X (Illumina, San Diego, CA, USA). All fastq files were passed through the quality control tool FASTQC and were mapped to the human genome (hg19) using STAR. The number of fragments aligned to each gene was obtained using featureCounts (v1.6.3, package from subread) based on the gencode v29 annotation file (gencode.v29lift37.annotation.gtf; ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/GRCh37_mapping). Normalized fragment per kilobase per million reads values were then calculated. Differentially expressed gene (DEG) analysis was performed to compare AR and parental NCI-H1975 cells using the DESeq2 R package. DEGs were defined as P-value < 0.05 and absolute log2 (fold change) ≥ 1. For DEG analysis comparing AR cells and TFAP2A knockdown AR cells, we used genes with an adjusted P-value < 0.05 and absolute log2 (fold change) ≥ 1 as significantly DEGs. Gene ontology (GO) enrichment annotation was performed using DAVID. One replicate of raw data from AR and parental cells was deposited in the EMBL database under accession numbers E-MTAB-5997 and E-MTAB-7068. Others were deposited in the GEO dataset under accession number GSE142031.

ATAC-seq library preparation and data analysis

ATAC-seq library processing was performed according to the Omni-ATAC protocol [22], which reduced the enrichment of mitochondria. The procedure generally included resuspending 50,000 viable cells and isolating nuclei; then, transposition was performed using Tn5 transposase (Vazyme, Nanjing, China), which was followed by adaptor ligation and PCR amplification (9 cycles). Next, the library was cleaned and concentrated, and quality control was performed via gel electrophoresis using an Agilent 2100 bioanalyzer. Libraries were sequenced with 150 bp paired-end on HiSeq X at 30 million raw reads each. All paired-end reads were first subjected to adaptor trimming using cutadapt (v1.18). Then, the clipped reads were aligned to the human genome (hg19) using bowtie2 (v2.3.3.1) with the parameters: -t -q -N 1 -L 25 -X 2000 no-mixed nodiscordant [23]. PCR duplicates were then removed using PicardTools (v1.119). For downstream analysis, nonuniquely mapped reads or reads mapped to the mitochondrial genome, Y chromosome, and unmapped contigs were removed [24]. To visualize the ATAC-seq signal in Integrative Genomics Viewer (IGV, v2.5.3), bam files were converted to bigwig (bw) files using deeptools (v3.3.0). To evaluate the correlation between ATAC-seq replicates, we adopted the method described by Wu et al. [23]. Briefly, the reads per kilobase of transcript, per million mapped reads (RPKM) value was calculated on a 100 bp genome window using featureCounts. The RPKM value was then summed within a 2 kb window and compared between replicates. A smooth scatter plot was then produced and the Pearson’s correlation coefficient was calculated using R (v3.6.0). Peaks were called for each sample using MACS2 (v2.1.1.20160309) in narrowPeak mode. Peaks overlapping with encoded backlist regions (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeDacMapabilityConsensusExcludable.bed.gz) were removed using bedtools (v2.26.0-148-gd1953b6). We performed differential chromatin accessibility (DA) analysis using the DiffBind package with the DESeq2 algorithm. DA was defined as a false discovery rate (FDR) < 0.05 and absolute log2 (fold change) ≥ 1. DA peaks were annotated to promoters or other genomic elements (5′-untranslated regions (5′-UTRs), 3′-UTRs, exons, introns, downstream regions, and intergenic regions) using ChIPseeker [25]. Promoter peaks were defined as −3 kb/+3 kb of nearest transcription start sites and other parameters were set to default. Distal regions were defined as upstream or downstream 2.5 kb of the defined promoter [23]. To identify TF motifs in AR-enriched peaks, a homer [26] (v4.10.4) de novo motif analyzer was used by setting parental enriched peaks as the background. AnnotatePeak.pl script in homer was used to infer the binding targets of TFAP2A.

For gene expression and chromatin accessibility correlation analysis, peaks overlapping with promoter regions and distal regions were selected. We associated the log2 (fold change) of these differential peaks with the log2 (fold change) of the same gene expression. Only those genes that were significant in both differential peak and differential expression analysis were considered. Raw data and processed data were deposited in the GEO dataset under accession number GSE142031.

TCGA survival analysis

RNA sequencing (RNA-seq) data and clinical data for NSCLC patients (lung adenocarcinoma (LUAD) and lung squamous carcinoma (LUSC)) were downloaded from the UCSC Xena (https://xenabrowser.net/datapages/). Normal controls in the expression matrix were excluded. The correlation analysis of expression and OS or disease-free survival (DFS) was performed by R package (survival). A survival plot was drawn with R (survminer).

Establishment of stable knockdown cell lines

RNA interference sequences targeting TFAP2A and a scramble control were synthesized and cloned into a linearized pLKO.1 plasmid (Addgene, USA). The interference primer sequences were as follows: shTFAP2A-F: 5′-CCGGGCTTGACCCACTTCAACCTCACTCGAGTGAGGTTGAAGTGGGTCAAGCTTTTTG-3′; shTFAP2A-R: 5′-AATTCAAAAAGCTTGACCCACTTCAACCTCACTCGAGTGAGGTTGAAGTGGGTCAAGC-3′; the scramble control was used as reported by Tan et al. [27]. Recombinant plasmids were cotransfected into 293T cells together with the packaging plasmids psPAX2 and pMD2G using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA). After 48 h, 293T supernatant was collected, concentrated with a Lenti-X concentrator (Takara, Kyoto, Japan), and used to infect AR cells. Puromycin (5 μg/mL, Thermo Scientific, Waltham, MA, USA) was used to kill nontransfected cells over the next month.

Endothelial cell migration assay

Endothelial cell migration was evaluated by wound-healing assays. HMEC-1 cells (1 × 106 cells per well) were seeded in six-well plates to form a monolayer the day before the experiment. Scratch wounds were made using a 20–200 μL pipet tip. Each well was washed with phosphate-buffered saline (PBS) three times, to remove loose cells, and 2 mL of conditioned medium (cell medium supplemented with 8 μg/mL anlotinib or PBS) was added, which was followed by a 5 h incubation at 37 °C. Then, 100 ng/mL BMP4 (R&D Systems, Minneapolis, MN, USA) was added to the shRNA negative control (shNC) medium. A phase-contrast microscope (Nikon, Tokyo, Japan) was used to capture images at 0 h and 5 h. The migration rate was calculated based on the change in wound width.

Tube formation assay

The angiogenic activity of the test medium was assessed by tube-formation assay. First, 10 μL of Growth Factor Reduced Matrigel (Becton Dickinson, UK) was added to the wells of μ-Slide Angiogenesis (Ibidi, Germany) and incubated at 37 °C for 0.5 h. Then, 2000 HMEC-1 cells suspended in 50 μL of conditioned medium were seeded in each well. After incubation at 37 °C for 10 h, images were captured using a phase-contrast microscope (Nikon, Tokyo, Japan) and tube-formation ability was evaluated by total network length, which was measured by the ImageJ Angiogenesis Analyzer plugin.

Chromatin immunoprecipitation

TFAP2A binding was validated by chromatin immunoprecipitation quantitative PCR (ChIP-qPCR) using a previously established protocol [28]. A total of 1 × 107 cells were cross-linked using 1% formaldehyde for 10 min. Then, chromatin was sheared to produce 150–350 bp fragments, which were immunoprecipitated with protein A + G magnetic beads (Millipore, MA, USA) coupled with 5 μg of TFAP2A antibody (ab52222, Abcam, MA, USA); incubation occurred at 4 °C overnight and was accompanied by rotation. The ChIP DNA and input DNA were used to perform qPCR with primers (Supplementary Table S1) spanning predicted TFAP2A-binding regions in multiple genes.

Quantitative real-time PCR

Total RNA was extracted with TRIzol (Thermo Scientific, Waltham, MA, USA) and reverse transcription reactions were performed following the protocol of the PrimeScript™ RT reagent Kit with gDNA Eraser (Takara, Kyoto, Japan). Quantitative real-time PCR (RT-qPCR) was performed on ABI Step One Plus (Applied Biosystems, Waltham, MA, USA) using SYBR Green (Bio-Rad, Hercules, USA). All primer sequences used for RT-qPCR are listed in Supplementary Table S1. RT-qPCR results were analyzed by using StepOne software (v2.1, Applied Biosystems, Waltham, MA, USA).

Statistical analysis

Functional in vitro assays and qPCR were performed three times for each condition. Statistical analysis was performed with GraphPad Prism 6 (San Diego, CA, USA) using Student’s t-tests. The data are expressed as the mean ± SEM and differences were considered significant at *P < 0.05, **P < 0.01, and ***P < 0.001.

Results

Transcriptome analysis of NCI-H1975 and anlotinib-resistant NCI-H1975 cells

We have previously established an AR NCI-H1975 cell line [29, 30]. To characterize the mechanism of acquired anlotinib resistance, RNA-seq was performed with both parental and AR cells. We identified 1273 significantly DEGs (496 upregulated and 777 downregulated genes in AR cells) (Fig. 1a and Supplementary Table S2). GO analysis of 496 upregulated genes indicated that angiogenesis-related biological processes (including blood vessel development, blood vessel morphogenesis, angiogenesis, response to hypoxia, wound-healing, vasculogenesis, and positive regulation of endothelial cell migration) were statistically enriched (Fig. 1b).

Fig. 1

Angiogenesis-related genes are upregulated in AR cells.

a Volcano plot of genes that are differentially expressed between parental and AR cells. b GO biological process analysis of upregulated genes in AR cells. c Heatmap of upregulated genes related to angiogenesis in AR cells. d Quantitative real-time PCR analysis of gene expression changes in representative angiogenesis-related genes.

The change in transcription of angiogenesis-related genes is shown in Fig. 1c. We then performed RT-qPCR with parental and AR cells to verify the RNA-seq results. Compared with the parental cells, a significant upregulation was observed in AR cells for all selected genes: BMP4, HSPG2, SHH, FOFX1, ADM, EDN1, NOV, SCG2, and TGFB2 (Fig. 1d). These results suggest that angiogenesis may play an important role in anlotinib resistance.

Identifying altered chromatin accessibility regions in anlotinib-resistant cells

Recent studies have reported the role of epigenetic changes in drug resistance, including DNA methylation [31], histone modifications [16], and chromatin accessibility [13]. We next hypothesized that chromatin accessibility is involved in anlotinib resistance. To characterize changes in chromatin accessibility in parental and AR cells, we performed an ATAC-seq assay. The ATAC-seq data indicated a periodic distribution of ∼200 bp of fragment insert size (Supplementary Fig. S1a) and there was a high correlation between replicates (Supplementary Fig. S1b), suggesting qualified ATAC-seq data. Analysis of the genomic distribution of accessible chromatin regions showed that distal intergenic regions had the largest proportion, which was followed by introns and promoters (Supplementary Fig. S1c).

To examine the chromatin accessibility differences between these two cell lines, differentially accessible peak analysis was performed. We identified 10,614 differential peaks (FDR < 0.05 and fold change ≥ 2) (Fig. 2a and Supplementary Table S3). We then annotated these peaks to the nearest gene and performed GO analysis of 2666 peaks that exhibited greater accessibility in AR cells than in parental cells (AR > parental). The results showed that GO items of blood vessel remodeling, vasculogenesis, and angiogenesis were highly enriched, which is consistent with the results of the transcriptome analysis (Fig. 2b).

Fig. 2

Chromatin accessibility analysis of WT and AR cells.

a Differential ATAC-seq peak analysis between parental and AR cells defined 2666 peaks that were greater in AR cells (peak concentration calculated with DiffBind is shown). b GO biological process analysis of AR > parental peaks. c Correlation between distal chromatin accessible regions and gene expression. d Chromatin accessibility at distal regions and gene expression coverage maps for BMP4 and HSPG2.

To assess the relationship between chromatin accessibility and gene expression, correlation analysis on peaks at the promoter or distal regions, and gene expression on genome-wide was conducted (see “Materials and Methods”). Consistent with the previous study [24], we found a significant correlation between accessibility of the promoter or distal regions and gene expression level (Fig. 2c and Supplementary Fig. S2a). We next visualized chromatin accessibility at distal regions and expression of genes known to be involved in angiogenesis. In AR cells, we observed greater chromatin accessibility at distal regions and higher expression levels of BMP4 (Fig. 2d), a gene involved in promoting endothelial cell migration and tube formation [32]. A similar pattern was also found at the HSPG2 locus (Fig. 2d).

Transcription factors associated with anlotinib resistance

Chromatin accessibility is a prerequisite for the binding of TFs. Recent studies have revealed the role of TFs in promoting tumor migration and metastasis through ATAC-seq [11, 12]. Thus, we examined possible TF motifs located in the chromatin regions that were more accessible in AR cells than they were in parental cells. Using the de novo TF motifs discovery software HOMER [26], we identified 21 TF candidates enriched at the regions of altered accessibility in AR cells (Supplementary Fig. S3). The top ten TFs are shown in Fig. 3a. Compared with the parental cells, we found that only TFAP2A was upregulated in AR cells (Fig. 3b and Supplementary Table S2). These results suggest that TFAP2A is potentially involved in anlotinib resistance.

Fig. 3

TFAP2A is the top selectively enriched motif at enriched peaks in AR cells.

a Top 10 motifs identified by homer. b Bar graph of log2 transformed FPKM values for all significantly enriched AR cell-specific TF motifs. c Disease-free survival analysis of TFAP2A expression in the TCGA LUAD cohort. d Overall survival analysis of TFAP2A expression in the TCGA LUAD cohort.

To evaluate the clinical relevance of TFAP2A, we performed a survival analysis on NSCLC patients from TCGA cohort. Kaplan–Meier analysis indicated that the high expression level of TFAP2A was significantly correlated with poorer OS and DFS outcomes in LUAD patients (Fig. 3c, d) but not in LUSC patients (Supplementary Fig. S4a, b). Comparison of the expression level of TFAP2A among different stages in NSCLC patients showed a higher expression level of TFAP2A in late-stage LUAD patients, but the changes were not significant (Supplementary Fig. S4c).

Downregulation of TFAP2A inhibits tumor-induced angiogenesis

TFAP2A has been implicated in cancer proliferation, invasion, and epithelial-to-mesenchymal transition (EMT) [33, 34]. However, it remains unclear whether TFAP2A is involved in tumor-induced angiogenesis. To uncover the underlying roles of TFAP2A in anlotinib resistance, we performed an short hairpin RNA assay and generated a TFAP2A knockdown cell line with AR cells to examine pro-angiogenic activity.

The pro-angiogenic effects were first investigated by a wound-healing assay. Compared with these cells cultured in shNC medium, HMEC-1 cells cultured with medium from TFAP2A knockdown AR cells (shTFAP2A) had a higher migration ability (Fig. 4a, b). After the introduction of anlotinib, the migration of HMEC-1 cells was remarkably decreased, whereas no significant decrease was observed between cells treated with shNC alone and those treated with shNC that were supplemented with anlotinib (Fig. 4a, b).

Fig. 4

TFAP2A promotes tumor-induced angiogenesis and anlotinib resistance.

a Representative images of wound-healing assays performed on HMEC-1 cells treated with conditioned medium (shNC medium + PBS, shTFAP2A medium + PBS, shNC medium + anlotinib, or shTFAP2A medium + anlotinib). b Quantification of wound closure in HMEC-1 cells treated with conditioned media. c Representative images of tube-formation assays performed with HMEC-1 cells in the presence of the indicated medium conditions. d Quantification of network length relative to HMEC-1 cells treated with the indicated medium. *P < 0.05, **P < 0.01.

To further evaluate the pro-angiogenic effects of TFAP2A, tube-formation assays were performed. Similar to the results of the wound-healing assay, treatment of HMEC-1 cells with medium from TFAP2A knockdown AR cells resulted in a decrease in the relative network length compared with cells cultured in shNC cell medium (Fig. 4c, d). Introduction of anlotinib to medium from TFAP2A knockdown cells decreased the tube-formation ability of HMEC-1 cells, whereas no significant decrease was observed between cells treated with shNC alone and those treated with shNC that were supplemented with anlotinib (Fig. 4c, d).

Validation of TFAP2A targets

To determine whether TFAP2A modulates the expression of genes involved in tumor-induced angiogenesis, we performed RNA-seq analysis of TFAP2A knockdown AR cells. DEG analysis showed that 2280 protein-coding genes were downregulated upon TFAP2A knockdown (Supplementary Table S4 and Supplementary Fig. S5a). These downregulated genes were enriched in pro-angiogenesis GO terms, such as platelet-derived growth factor receptor signaling pathway (PDGFR), transforming growth factor-β receptor signaling pathway (TGF-β), vascular endothelial growth factor receptor signaling pathway (VEGFR), and sprouting angiogenesis (Fig. 5a). We then computationally identified 1073 binding target genes of TFAP2A in the regions with increased accessibility in AR cells (Supplementary Table S5).

Fig. 5

Identification and validation of TFAP2A-binding targets.

a GO analysis of downregulated genes with TFAP2A knockdown in AR cells. b Venn plot of downregulated genes with TFAP2A knockdown and putative TFAP2A-binding targets. c IGV plot for ATAC-seq at the BMP4 and HSPG2 loci. Putative TFAP2A-binding targets are shown with motif logos. d Chromatin immunoprecipitation PCR validation of TFAP2A-binding targets. e Quantitative real-time PCR validated the downregulation of BMP4 and HSPG2 when comparing the TFAP2A knockdown and shNC groups.

Intersecting with RNA-seq data, we found an overlap of 171 genes that were putatively bound by TFAP2A and downregulated in TFAP2A knockdown AR cells (Fig. 5b), including angiogenesis-related genes, such as BMP4, HSPG2, and SHH. To validate whether these potential targets are physically bound by TFAP2A, we performed a ChIP-qPCR assay and observed the significant binding enrichment of TFAP2A at the selected accessible regions of BMP4 (distal) and HSPG2 (intron) (Fig. 5d). Accordingly, we found that the expression of the TFAP2A-binding targets BMP4 and HSPG2 was downregulated (Fig. 5e). BMP4 has been reported to induce tube formation and the migratory efficiency of microvascular endothelial cells [32]. Here we also found that the reduced pro-angiogenic ability of TFAP2A knockdown AR cells was partially rescued by supplementation with exogenous BMP4 (Supplementary Fig. S6).

Discussion

Anlotinib is a multitarget agent used for blocking tumor angiogenesis and proliferative signaling. The targets of anlotinib include VEGFR (2/3), PDGFR (α/β), and fibroblast growth factor receptor (1–4) [35, 36]. Our previous studies have reported that anlotinib as a third-line treatment can prolong PFS and OS in NSCLC patients [4, 5, 37, 38]. However, acquired drug resistance is inevitable in advanced NSCLC patients.

Angiogenesis, the formation of new blood vessels from the preexisting vascular network, is a hallmark of cancer [39]. During cancer progression, tumor cells secrete pro-angiogenic factors (including VEGF, PDGFB, TGF-β, matrix metalloproteinases, and chemokines), which directly or indirectly induce the sprouting, migration, and proliferation of endothelial cells [40]. Our previous study reported that the overexpression of the ELR + CXC chemokine CXCL2 in AR cells is involved in anlotinib resistance [29]. In this study, we observed the upregulation of some pro-angiogenic factors in AR cells, including BMP4, HSPG2, SHH, FOFX1, ADM, EDN1, NOV, SCG2, and TGFB2 (Fig. 1c, d). These results suggest that some upregulated pro-angiogenic factors are secreted and contribute to the pro-angiogenic ability of AR cells.

The alteration in chromatin accessibility has been reported to be involved in acquired resistance against tumor drugs [13]. In this study, we performed chromatin accessibility analysis in parental and AR cells via ATAC-Seq. Differential accessibility analysis revealed that greater accessibility regions in AR cells were enriched in angiogenesis-related processes, such as TGF-β signaling, blood vessel development, PDGF signaling, and the cellular response to hypoxia (Fig. 2b). TFs play an important role in tumoral drug resistance [14-18]. By analyzing the accessibility of chromatin, we can not only identify the regulatory regions but also infer TF motifs. We performed de novo motif analysis of the greater accessible regions in AR NCI-H1975 and identified motifs of 21 TF candidates. We found that only TFAP2A is upregulated in AR cells.

Previous studies reported that TFAP2A promotes EMT procession by regulating TGF-β signaling in cancer cells [33] and TFAP2A regulates tumor growth via hypoxia inducible factor-1a (HIF-1a) signaling in nasopharyngeal carcinoma and NSCLC [34, 41], whereas the HIF-1a pathway has been confirmed to promote angiogenesis [42]. To validate the role of TFAP2A in anlotinib resistance, we performed TFAP2A knockdown in AR cells and found that tumor-induced angiogenesis was decreased in AR cells by TFAP2A knockdown. Meanwhile, the anti-angiogenic ability of anlotinib was partially rescued in TFAP2A knockdown cells. We identified 2280 downregulated genes in TFAP2A knockdown AR cells and found that they were enriched in the PDGF, EGFR, TGF-β, and VEGF signaling pathways (Fig. 5a), suggesting the importance of TFAP2A in angiogenesis. To understand how TFAP2A regulates its direct transcriptional targets related to angiogenesis, we computationally identified 1073 putative binding genes and found that the expression of 171 genes might be regulated by TFAP2A (Fig. 5b), including the angiogenesis-related genes BMP4, HSPG2, and SHH. We finally demonstrated that TFAP2A binds and regulates BMP4 and HSPG2. Previous studies have reported that BMP4 induces tube formation and the migratory efficiency of microvascular endothelial cells [32], and BMP4 induces an increase in VEGFR2 expression, which mediates endothelial cell activation [43]. Our results also showed that the decreased pro-angiogenic ability of TFAP2A knockdown with AR cells was partially rescued by supplementation with BMP4, which suggests that TFAP2A promotes angiogenesis by directly upregulating BMP4; in turn, this upregulation causes the migratory efficiency of microvascular endothelial cells and further promotes angiogenesis and anlotinib resistance.

Briefly, in our chromatin accessibility analysis, we found that the TF TFAP2A plays an important role in anlotinib resistance by directly regulating the expression of BMP4.

Supplementary materials

Supplementary Table S2

Supplementary Table S3

Supplementary Table S4

Supplementary Table S5