B‐cell lymphoma 2 family genes show a molecular pattern of spatiotemporal heterogeneity in gynaecologic and breast cancer

Abstract Objectives BCL2 family proteins have been widely studied over the past decade due to their essential roles in apoptosis, oncogenesis and anti‐cancer therapy. However, the similarities and differences in the spatial pattern of the BCL2 gene family within the context of chromatin have not been well characterized. We sought to fill this knowledge gap by assessing correlations between gene alteration, gene expression, chromatin accessibility, and clinical outcomes in gynaecologic and breast cancer. Materials and methods In this study, the molecular characteristics of the BCL2 gene family in gynaecologic cancer were systematically analysed by integrating multi‐omics datasets, including transcriptomics, chromatin accessibility, copy number variation, methylomics and clinical outcome. Results We evaluated spatiotemporal associations between long‐range regulation peaks and tumour heterogeneity. Differential expression of the BCL2 family was coupled with widespread chromatin accessibility changes in gynaecologic cancer, accompanied by highly heterogeneous distal non‐coding accessibility surrounding the BCL2L1 gene loci. A relationship was also identified between gene expression, gene amplification, enhancer signatures, DNA methylation and overall patient survival. Prognostic analysis implied clinical correlations with BAD, BIK and BAK1. A shared protein regulatory network was established in which the co‐mutation signature of TP53 and PIK3CA was linked to the BCL2L1 gene. Conclusions Our results provide the first systematic identification of the molecular features of the BCL2 family under the spatial pattern of chromatin in gynaecologic and breast cancer. These findings broaden the therapeutic scope of the BCL2 family to the non‐coding region by including a significantly conserved distal region overlaying an enhancer.


| INTRODUC TI ON
Each year, about 1 million new cases of gynaecologic malignant tumours are diagnosed and approximately 0.5 million women die as a result. 1,2 The formation and evolution of gynaecologic cancers are affected by many factors, including heredity, lifestyle, diet, exercise, sex and geographical environment. 1 Gynaecologic and breast cancers, which are distinct classes of tumours, share some common characteristics, such their early formation from Mullerian ducts and their developmental control by female hormones. 3 5 This especially includes the B-cell lymphoma 2 (BCL2) family proteins, which are critical mediators of the apoptotic response and can be targeted by "BH3-mimetic" small molecule compounds. 6 Since the discovery of the BCL2 family in 1984, 7 its members have expanded to include more than 15 that are involved in human cancer. 8 The BCL2 family can be categorized into three groups: pro-survival BCL2 family members; pro-apoptotic multi-BHdomain members; and pro-apoptotic "BH3-only" members. 6 Tumour cells frequently overexpress BCL2 family proteins to escape the apoptosis checkpoint. 8 Copy number variations have been detected in the anti-apoptotic members MCL1 and BCL2L1 across 26 human cancers, including gynaecologic cancers. 9 Furthermore, overexpression of MCL1 is thought to promote chemotherapy resistance and prolong cell survival in breast cancer and ovarian cancer. 10,11 However, loss of the representative pro-apoptosis members BAX, BIM and BBC3 is also a common trend in oncogenesis, facilitating tumour formation and progression through genomic deletion, silencing and mutation in several cancers. [12][13][14] Despite the clear roles for these proteins, the correlation between BCL2 family members has not been systematically characterized using integrated multi-omics data sets in gynaecologic cancers.
Previous research on the BCL2 family has been mostly concerned with overexpression or copy number variation and seldom with long-distance regulation, and the distal regulation of the BCL2 family in gynaecologic cancers remains ambiguous. 6,15 Distal DNA regulatory elements can dramatically increase the expression of oncogenes/tumour suppressor genes during tumorigenesis. 16,17 However, recognizing and characterizing their contributions can be difficult. Because of polyA selection, relatively low depth genome sequencing and specific enhancer usage in rare cell populations, some enhancers have not been identified. 16,17 Recently, the assay of transposase-accessible chromatin with sequencing (ATAC-seq) method has proven useful for identifying functional spatiotemporaland tissue-specific distal regulation mechanisms between promoters and enhancers on a genomic scale. 18 Distal regulatory regions of BCL2 have been shown to target promoters based on chromatin accessibility dynamics in breast cancer. 19,20 However, potential enhancers of the BCL2 family are lacking in systematic research of gynaecologic cancer.
Apoptosis signalling pathways form a dynamic and complex network that maintains organismal homeostasis. 21 Targeting BCL2 family proteins is an important approach in oncotherapy that strives to modulate the balance of apoptosis in order to control tumour cell death. Inhibitors of the BCL2 family (the BH3 mimetic drugs), such as ABT-263 (Navitoclax) and ABT-199 (Venetoclax), which target anti-apoptotic members, have shown promising results in clinical trials of haematological cancers, but they have not yet been approved for clinical use for solid cancers, including gynaecologic cancer. 5,22,23 Furthermore, BH3 mimetic drugs do not account for the great diversity in haematological and gynaecologic cancers, including the variable dependency on apoptotic proteins, and they vary in their safety and efficacy in the clinic. Suppression of the BCL2 family members may trigger disorganization of down/upstream signalling regardless of tumour heterogeneity, though the complex regulatory networks, especially the epigenetic signal and distal regulatory elements, are not well characterized.
To promote the development and application of drugs that target apoptosis signalling pathways, it is critical to further investigate molecular characteristics of the BCL2 family using large patient samples. Therefore, we evaluated differences in the patterns of genetic alteration by pro/anti-apoptotic genes using a multi-omics approach.
Our strategy involved analysis of the transcriptome, copy number variation, the methylome and clinical outcomes. We also evaluated associations between methylation, gene amplification and chromatin accessibility, including spatiotemporal patterns and tumour heterogeneity. Finally, we identified genetic alterations in the BCL2 family network across gynaecologic cancers. An improved understanding of the common features and variations in the regulation patterns of the BCL2 family will be beneficial for selecting molecular targets and mechanisms to accelerate programmed cell death, leading to novel targeted drugs for treating gynaecologic cancer.

| Determination and comparison of genetic alterations
In this study, gynaecologic cancer mainly included breast invasive carcinoma (BRCA), high-grade ovarian serous cystadenocarcinoma (OV), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), uterine corpus endometrial carcinoma (UCEC) and uterine carcinosarcoma (UCS). Our analysis using cBioPortal (The cBio Cancer Genomics Portal) was performed in four stages: select data sets, select genomic profiles, define sample sets and type the interesting genes. 8 Overall survival and BCL2 family gene alterations in pan-cancers, including somatic mutations, amplification, deep deletion and DNA methylation, were measured through the cBioPortal database. We chose the "individual cancer study" query option within cBioPortal to explore genomic alterations for the following genes:

| Exploring gene expression levels between normal tissues and tumour samples
Gene Expression Profiling Interactive Analysis (GEPIA) is a free and interactive online database for cancer mRNA data, which consists of 9736 tumours and 8587 normal samples from The Cancer Genome Atlas (TCGA) and Genotype-tissue Expression dataset projects (GTEx). 24 The expression levels of BCL2 family genes in gynaecologic tissues and normal samples were measured using the Boxplot module in GEPIA. 25

| Evaluating chromatin accessibility of the BCL2 gene family
The assay of transposase-accessible chromatin with sequencing (ATAC-seq) can sensitively map chromatin accessibility of active genes and also can identify functional spatiotemporal and tissuespecific regulatory mechanisms between promoters and distal sites within a 500 kb genomic region surrounding the transcriptional start site (TSS). 26 To explore and affirm the chromatin state of BCL2 family gene loci in gynaecologic cancer, we compared differential chromatin-accessibility (ATAC-seq) data sets from the UCSC Xena browser with published histone modification ChIP-seq datasets from gynaecologic tumour samples, including EP300, TEAD4, H3K4me1 and H3K27ac (Table S1). To highlight both differences and similarities between BCL2 family gene characteristics in gynaecologic cancers, we selected male-specific tumour testicular germ cell tumours as a non-gynaecologic cancer. The published histone modification data sets were collected in public databases including Gene Expression Omnibus (GEO), Encyclopedia of DNA Elements (ENCODE), and the Roadmap Epigenetics Project. The analysis of ChIP-seq data sets was divided into three parts: mapping the reads, calling the peaks and annotating the peaks. The raw data of SRA files that were downloaded from GEO (ftp://ftp-trace.ncbi.nih.gov/sra/srain stant /reads / ByRun /sra/SRR/) used fastq-dump (https://ncbi.github.io/sra-tools /fastq -dump.html) to convert SRA files to FASTQ files. First, we mapped these FASTQ files to the human genome (hg38) using BWA with the BWA-MEM algorithm. 27 Then, we used MACS2 to call the statistically significant peaks with the macs2 callpeak command. 28,29 Finally, we annotated the open-/closed peaks, including average conservation regions using HOMER with annotatePeaks.pl module. 30 The distal chromatin accessibility of the BCL2 family in gynaecologic cancer was analysed and visualized with the UCSC Xena Browser and the Integrative Genomics Viewer. 31

| Network analysis of the BCL2 family
To explore the BCL2 family signalling networks in gynaecologic cancer, we queried the BCL2 gene in cBioPortal. 32

| Statistical analysis
Statistical analyses of mRNA expression levels were performed in GEPIA. Differences between two groups were compared using t test, and differences between multiple groups were compared using one-way ANOVA. Family genes with P values <0.01 are considered significantly differentially expressed genes. Pearson/Spearman tests and linear regression equations were employed to investigate the correlation between DNA methylation and mRNA expression. If the Pearson/Spearman coefficient is <0, it indicates that DNA methylation and mRNA expression are negatively correlated. Linear regression equations can be utilized to show a gradual downward tendency.
Kaplan-Meier curves were prepared to evaluate overall survival, and the log-rank test was used to calculate P values. P < .05 was regarded as statistically significant in Kaplan-Meier survival curve analyses.

| Gene alteration of the BCL2 family in pancancer
BCL2 family proteins have been referred to by different names in different studies. To facilitate our study of the BCL2 family, we adopted a unified terminology with gene names from NCBI (Table S2) Therefore, we surveyed their genomic alterations in TCGA pan-cancer datasets. The gene alteration frequencies of the BCL2 family, including mutations, deletions and amplifications, are shown ( Figure 1A; Tables S3). The data show that gene alteration frequencies in the BCL2 family were more widely found in gynaecologic cancer than many other cancers. In general, BCL2 family amplification was more frequent than deep deletion and missense mutations across pan-cancer. This was particularly notable for gynaecologic cancers, including F I G U R E 1 Gene alteration of the BCL2 family in pan-cancer. A, The alteration frequencies of the BCL2 family across the TCGA PanCancer atlas. The horizontal axis represents the types of cancer, and the vertical axis represents the alteration frequencies of the BCL2 family. Green coding indicates non-synonymous mutations, purple coding indicates gene fusions, red coding indicates gene amplification, blue coding indicates deep deletions and grey coding indicates multiple alterations. Gynaecologic cancers are marked in red font. B, Genetic alteration of the BCL2 family in gynaecologic cancer. The functional classification of the BCL family member is indicated: a Activators; b Effectors; c Guardians; d Initiators; e Sensitizers UCS (35.09%), OV (18.46%) and UCEC (11.91%). More than 20% of cases in these data sets had BCL2 family amplification. Given the roles of BCL2 family proteins in regulating apoptosis, these results suggest that they are integrally involved in gynaecologic cancer.
To further explore the genetic alteration of the BCL2 family in gynaecologic cancer, we generated an Oncoprint of their alteration frequencies using cBioPortal. UCS, CESC, OV, BRCA and UCEC showed 8.92%,7.71%, 7.67%, 5.93% and 3.55% average alteration rates, respectively ( Figure 1B). Interestingly, we observed low mRNA expression in some patient samples of OV and BRCA, especially for the pro-apoptotic member. On the contrary, the alterations in the anti-apoptotic members, including BCL2L1 and MCL1, predominantly included mRNA upregulation and gene amplification. Consistently, driver mutations and downregulation of mRNA seldom were observed in gynaecologic patient samples for BCL2L1 and MCL1. These results suggest that for the gynaecologic cancers, anti-apoptosis BCL2 family members tended to be upregulated rather than downregulated, with gene amplification at higher frequency for the anti-apoptosis compared with the pro-apoptosis members.

| The mRNA expression of the BCL2 family is coupled with widespread chromatin accessibility changes in gynaecologic cancer
To assess the overall consequence of BCL2 family upregulation on chromatin accessibility, we first evaluated the expression changes in gynaecologic cancer using TCGA and GTEx databases. The different patterns of mRNA expression between the pro/anti-apoptotic members in gynaecologic cancer are shown in Figure 2 and showed upregulation in gynaecologic cancer. Therefore, there was a unique pattern of regulation of BCL2 family gene expression for the pro/anti-apoptotic members, with the pro-apoptotic members more highly expressed in gynaecologic cancer compared to normal tissues.
The presence of chromatin accessibility changes accompanied by H3K4me1 or H3K27ac histone modification outside promoter regions is indicative of an active enhancer. 33-37 Therefore, we explored whether the links between distal regulatory elements and promoters were specifically correlated with histone modifications.
To this end, we used ATAC-seq, which accurately identifies tissue-or stage-specific distal links regardless of the number of samples. We

| The relationship between promoter DNA methylation and BCL2 gene expression in gynaecologic cancer
To further investigate the regulatory mechanism of BCL2 family expression in gynaecologic cancer, we examined the correlation be-

| Prognostic analysis of selected genes with the enhancer signature in gynaecologic cancer
To explore the effect of BCL2 family alteration on patient survival in gynaecologic cancer, we compared the overall survival for cases with and without gene alteration ( Figure 5). Kaplan-Meier analysis indicated that gene alteration in BOK was significantly associated with favourable prognosis in OV (P < .1), which is consistent with the results of the genetic alteration analysis (Figure 1). Similarly, favourable prognosis was found for BAK1 in CESC (P < .05) and HRK in BRCA (P < .05). Nevertheless, gene alterations in BAD and BIK were associated with poor survival in UCEC (BAD, P < .05; BIK, P < .01). These results are supported by the high expression of BAK1 in CESC and BAD and BIK in UCEC ( Figure 2). However, no correlation between mRNA expression and survival was observed for BCL2L1.

| The regulatory network of the BCL2 family in gynaecologic cancer
The dynamic regulatory network of the BCL2 family protein suggests the importance of mitochondrial outer membrane permeabilization (MOMP), but the mechanisms behind the interactions remain controversial. 15,21,38 To identify the pattern of BCL2 family protein inhibition of MOMP, we comprehensively inferred the pan-cancer BCL2 family interaction signal network in gynaecologic cancer ( Figure 6A-E). The resulting networks from cBioportal analysis were highly structured across gynaecologic cancer

| D ISCUSS I ON
Apoptosis is a continuous programmed cell death event that is responsible for periodic control of damaged cells during normal development, maintaining a balance of organismal homeostasis and preventing pathological autoimmunity and tumorigenesis. 6 Unfortunately, apoptosis inhibition or escape in tumour cells leads to abnormal survival and accumulation of dysfunctional cells. 8 Apoptosis is controlled by the BCL2 protein family, which includes pro-apoptotic and pro-survival members. Some members of the family, such as BCL2, are overexpressed in tumour cells, thus breaking the balance between life and death and resulting in unlimited proliferation 8,19 ( Figure 6I). Furthermore, the regulatory relationships among these proteins and DNA amplification, mutation, methylation, overall survival and chromatin accessibility have not been comprehensively examined in gynaecologic cancers.
Here, we evaluated the molecular characteristics of BCL2 family genes under the spatial pattern of chromatin in gynaecologic cancers. As noted, most prior studies have identified initiators, effectors and guardians of the BCL2 family within a unified dynamic model in a variety of cancers, leading to anti-apoptosis and drug resistance in tumour cells 21 ( Figure 6I). Recent systematic analyses have uncovered the dynamics of BCL2 family regulation, particularly for anti-apoptotic members, as the most frequent somatic copy-number alteration in human cancer. 9,15,21 We found that BCL2L1 and MCL1 show higher expression levels as compared to other BCL2 members, which is consistent with the genetic alteration results. Our observation that anti-apoptotic members are not widely expressed in gynaecologic cancers agrees with previous studies, and low expression level and deep deletion was also observed in some patient samples for anti-apoptotic members. 39,40 In terms of distal non-coding chromatin accessibility, both BCL2 and BCL2L1 have long-range promoter interactions that affect gene expression. Conversely, BCL2, BCL2L2, and MCL1 show downregulation across 2-3 different cancer types, despite distal BCL2 regulation. 41 These phenomena may be due to the enforced expression or activation of miRNAs that normally suppress BCL2 family expression, such as miR-15a or miR-16.1 targeting of BCL2; miR-29, miR-125 and miR-193 targeting of MCL1; or let-7 targeting of BCL2L2. [41][42][43][44] Of the remaining overexpressed members, BAX, NOXA, BIK, BID and BAK are pro-apoptotic. Enhancer signatures can be detected, and the mRNA expression levels of these pro-apoptotic genes may be affected by these distal signatures.
However, active effectors (BAX and BAK) or initiator (PMAIP1, BIK and BID) cannot lead to apoptosis in tumours, and these pro-apoptotic proteins may be sequestered by guardians (BCL2 and BCL2L1) that have more potent enhancer activity to propel the anti-apoptotic mechanisms.
In our study of overall survival, BOK in OV, BAD and BIK in UCEC, BAK1 in CESC, and HRK in BRCA show noteworthy associations with overall survival. BAD and BIK regulate tumour growth in many cancers. [45][46][47][48] Furthermore, worse overall survival may be triggered by activated guardians (the antiapoptotic proteins eg BCL2L1), which can bind and neutralize BAD and BIK to mediate the inhibition of apoptosis. 15,21,38 BAD phosphorylation has recently been reported to promote tumour cell survival, and post-translational modification might therefore contribute to impaired pro-apoptotic proteins. 49 Overexpression of BAK1 has previously been associated with a favourable prognosis in breast cancer, 50 while BOK, BAK1 and HRK as effectors are inserted into the mitochondrial outer membrane, resulting in MOMP without sequestration by the guardians. 15,21,38 F I G U R E 5 The overall survival in gynaecologic cancer with BCL2 family alteration. Alterations in BOK in OV, BAD and BIK in UCEC, BAK1 in CESC, and HRK in BRCA were significantly associated with overall survival Overall, our results imply that BAD, BIK and BAK1 may be prognostic genes for clinical effects in gynaecologic cancer.
It is noteworthy that deep deletion and mutation were not observed widely in gynaecologic cancers, especially for pro-apoptotic genes, which differs from previous studies. 6 Beyond a limited number of genetic differences, the network from cBioportal analysis was highly structurally similar for different gynaecologic cancers. Our discovery of co-mutation signatures for TP53 and PIK3CA with BCL2L1 is novel and was not revealed in a previous study of gynaecologic cancer. 55,56 Co-mutation of TP53 and PIK3CA is a primary mediator of anti-apoptotic protein inhibition of MOMP that lead to a more aggressive phenotype with a worse prognosis in breast cancer. 56 These proteins are connected by co-mutation signatures; MCL1 and BCL2L1 transcriptional activity are consistent with chromatin accessibility of the region surrounding the gene loci with more than one predicted distal enhancers.
Therefore, we hypothesize that the co-mutation signature of TP53 and PIK3CA acts to stimulate the formation of chromatin loops of BCL2L1 gene loci to strengthen its transcriptional activity. BCL2L1, as the only upregulated anti-apoptotic member, prevents the activation and oligomerization of pro-apoptotic members that act as guardians against the effectors or activators on the mitochondrial membrane, which maintains the balance of mitochondrial membrane potential and thus prevents the release of cytochrome c to suppress caspases. 39,57-59 Using BCL2L1, TP53/TP63 and PIK3CA as possible potential prognosis markers may therefore be an effective approach for diagnosis and treatment of gynaecologic cancer.
The resulting network models showed specific differences in YWHAZ amplification in BRCA, PIK3CB amplification in CESC, PIK3R1 mutation in UCEC, and downregulation of XRCC6, NMT1 and CASP3 in OV, which suggests that the BCL2 family protein network can be used to identify different types of gynaecologic cancer. YWHAZ has been shown to stimulate lung cancer cell proliferation and metastasis and promote the invasion of breast cancer cells, 60 suggesting that it might serve as a therapeutic target of breast cancer. 61 Mutation of PIK3CB, as the catalytic subunit in the PI3K signalling pathway, drives tumour cell growth and migration. 62 PIK3CB has been reported as a selective survival factor in glioblastoma. 63 Furthermore, co-mutation of PIK3R1 and PIK3CA is associated with oncogenesis and hyperactivity of the PI3K signal pathway in breast cancer, supporting an oncogenic role of the co-mutation pair. 64 Loss of PIK3R1 is an effective therapeutic mechanism for PIK3CA-positive breast cancers. 65 On the other hand, activation of CASP3 is involved in the initiation of cell apoptosis, 66 inhibition of NMT1 regulates breast cancer oncogenesis by the JNK pathway, 67 and inactive XRCC6 fails to protect genomic integrity. 68 Therefore, our findings further validate previous studies demonstrating that downregulation of XRCC6, NMT1 and CASP3 is significantly associated with tumorigenesis.
Long-range enhancer-promoter gene expression is facilitated and constrained by the 3D architecture of mammalian genomes, which plays a key role in disease. 69 We demonstrated that the sig- In conclusion, as the first systematic analysis of molecular feature of the BCL2 family under the spatial pattern of chromatin in gynaecologic cancer, our study broadens the therapeutic scope of the BCL2 family to the distal non-coding region. We demonstrated that differential expression of BCL2 family members occurs at different frequencies. Furthermore, we identified the relationship between overall survival, enhancer signature, gene amplification and DNA methylation. Our results also establish a shared protein regulatory network in which the co-mutation signatures of TP53 and PIK3CA interact with BCL2L1, which provides a new strategy for biomarker identification in oncotherapy.

CO N FLI C T O F I NTE R E S T
All authors declare that they have on conflict of interest.

AUTH O R CO NTR I B UTI O N S
Jiajian Wang and Hongli Du planned the study. Jiajian Wang and Sidi Li performed the gene expression and chromatin accessibility F I G U R E 6 The regulatory network of the BCL2 family in gynaecologic and breast cancer. A-E, The regulatory network of the BCL2 family members with 30% alteration frequency. The networks from cBioportal analysis show alteration frequency, mRNA expression, copy number and mutation coding by various colours in different circular areas across gynaecologic cancers. F-H, The co-mutation/amplification signatures for TP53, PIK3CA and BCL2L1, for which gene mutation or amplification/upregulation was found more frequently. I, The BCL2 family-mediated mitochondrial apoptosis pathway in gynaecologic and breast cancer. In the death signal pathway, upstream signalling unleashes initiator proteins, including sensitizers (BIK, PMAIP1, etc) and activators (BID, etc). Subsequently, these proteins are bound and immediately sequestered by guardian proteins (including BCL2L1). Because the transcriptional expression of the BCL2L1 gene tends to have more amplification events and chromatin accessibility, BCL2L1 protein never reaches a state of saturation. Guardian proteins prevent effector proteins from oligomerizing and causing mitochondrial outer membrane permeabilization (MOMP), which prevents cytochrome c release from the mitochondria. Without cytochrome c, APAF1 cannot dimerize with procaspase-9 to form the apoptosome. Thus, guardian proteins mediate the immortalization of cells of BCL2 family, generated the results of DNA methylation and overall survival, constructed the regulatory network, designed the ana-