Effect of BRCA1 on the Concurrent Chemoradiotherapy Resistance of Cervical Squamous Cell Carcinoma Based on Transcriptome Sequencing Analysis

Background Cervical squamous cell carcinoma (CSCC) is the main pathological type of cervical cancer, accounting for 80%–85% of cervical cancer. Owing to concurrent chemoradiotherapy (CCRT) resistance in a subset of CSCC patients, the treatment response is often unsatisfactory. Identifying predictors and therapeutic targets related to cisplatin-based CCRT resistance in CSCC is critical. Methods We reanalyzed GSE56363, an mRNA dataset from the GEO database with 21 patients with locally advanced CSCC, to identify differentially expressed genes (DEGs) related to CCRT resistance. The hub genes were screened from the protein-protein interaction network of DEGs using cytoHubba plug-in of Cytoscape software. Transcriptome sequencing technology was used to compare differential expression between SiHa cells overexpressing BRCA1 compared with control SiHa cells. Functional annotation for DEGs and gene set enrichment analysis (GSEA) was performed to identify DEG-enriched relative signaling pathways to examine the molecular mechanisms of BRCA1 in CCRT resistance of CSCC. qPCR was used to verify the expression of key genes in SiHa/DDP cells. Results A total of 609 DEGs including 223 upregulated DEGs and 386 downregulated DEGs were identified between the complete response to CCRT (CR) and noncomplete response to CCRT (NCR) CSCC patients based on the GSE56363 dataset. Ten hub genes with the highest degrees were identified via the plug-in CytoHubba in Cytoscape: BRCA1, CDCA8, ASPM, CDC45, RAD51, HMMR, CENPF, EXO1, DTL, and ZWINT genes, and BRCA1 ranked first. Through transcriptome sequencing analysis based on GSE141558, 1344 DEGs were identified in BRCA1-overexpressing SiHa cells, including 824 upregulated DEGs and 520 downregulated DEGs. GSEA results showed that CCRT-resistance related signaling pathways, such as the JAK/STAT signaling pathway and the WNT signaling pathway, were differentially enriched in BRCA1-expressing SiHa cells. STAT1, STAT2, and CCND1 were screened as the differentially expressed target genes of BRCA1 and may correlate with resistance of CSCC. qPCR results showed that only STAT1 was significantly increased in SiHa cells with GV230-BRCA1 plasmid transfection. Conclusion BRCA1 overexpression in SiHa cells may upregulate STAT1 to activate the JAK/STAT signaling pathway, suggesting a mechanism for enhanced CCRT resistance.


Introduction
Cervical cancer (CC) ranks fourth for both incidence and mortality among cancers in women worldwide, and there were an estimated 570,000 new cases and 311,000 deaths related to CC in 2018 [1]. Cervical squamous cell carcinoma (CSCC) is the main pathological type of CC, accounting for 80%-85% of all CC cases. Approximately two-thirds of CSCC patients are diagnosed at medium or advanced stage, and the standard treatment is a concurrent chemoradiother-apy (CCRT) regimen consisting of cisplatin-based chemotherapy, external beam radiotherapy, and brachytherapy [2]. However, owing to CCRT resistance in a subset of patients, the treatment response is often unsatisfactory. Therefore, identifying efficient predictors and therapeutic targets related to CCRT is critical to improve the prognosis of CSCC patients.
With the rise of bioinformatics and high-throughput sequencing technologies, cancer development mechanisms, drug resistance, and prognosis can be better explained at the genetic level [3,4]. These findings will contribute to understanding individual differences and tumor heterogeneity in cancer patients. These technologies have been applied to investigate predictors and therapeutic targets of radiotherapy and chemotherapy resistance in many tumor types [5][6][7]. However, relatively few studies have examined CCRT resistance in CSCC [8,9].
The BRCA1 tumor suppressor responds to DNA damage by participating in DNA repair, mRNA transcription, cell cycle regulation, protein ubiquitination, and other cellular pathways [10,11]. The DNA repair ability of damaged cells improves with increased intracellular BRCA1 expression, resulting in drug resistance. Multiple studies have shown that the expression level of BRCA1 gene in cancers is negatively correlated with the sensitivity to platinum drugs [12][13][14] and sensitivity to radiotherapy [15][16][17].
In the current study, we screened the differentially expressed genes (DEGs) related to CCRT resistance of CSCC and found that BRCA1 overexpression enhanced CCRT resistance in CSCC. Subsequently, we performed transcriptome sequencing to explore the potential mechanisms by which BRCA1 may enhance resistance in CSCC. Our findings may provide theoretical support and new strategies for the diagnosis and treatment of CSCC patients with CCRT resistance.

Materials and Methods
2.1. Identification of DEGs Related to CCRT Resistance in CSCC Based on the GEO Dataset. The GSE56363 database of CSCC patients was downloaded from the NCBI-GEO database [18]. The CSCC patients were treated by CCRT, and cisplatin-based treatments were used as chemotherapeutics. Patients were divided into two groups: complete response to CCRT (CR) and noncomplete response to CCRT (NCR). The DEGs between NCR and CR patients were reanalyzed using limma package in R [19]. False discovery rate ðFDRÞ < 0:05 and |LogFoldChangeðFCÞ | >1 were considered as the cut-off criterion.

Gene Ontology (GO) Enrichment Analysis and Construction of the Protein-Protein Interaction Network.
We performed GO enrichment for the DEGs related to CCRT resistance in CSCC by clusterProfiler package [20]. The STRING database (https://string-db.org/) is an online resource that provides known direct and indirect proteinprotein interaction (PPI). We mapped the DEGs to STRING [21] to reveal the PPI information and visualized the PPI network using Cytoscape 3.7.2 [22]. The combined scores ð interaction scoresÞ > 0:4 were considered as significant. We then analyzed the PPI network to screen the hub genes using the cytoHubba plug-in of Cytoscape software [23].
We selected the top 10 DEGs with the highest degrees for further analysis. The receiver operating characteristic (ROC) curves of the hub genes were generated by the pROC package in R [24]. The violin plots of the hub gene expression were developed by the ggplot2 package in R [25]. Heat maps of the hub gene expression were developed by the pheatmap package in R [26].
2.3. Cells, Cell Culture, and Plasmid Transfection. The SiHa human CSCC cell line was purchased from Cell Bank of the Chinese Academy of Sciences, and SiHa/DDP cells were purchased from Fenghui Bio-Technology Limited. The cells were cultured in DMEM medium (Gibco, USA) supplemented with 10% fetal bovine serum (Sangon, China) and 1% antibiotic-antimitotic solution (Gibco). The GV230-BRCA1 and GV230 plasmids were purchased from Gene-Chem (Shanghai, China). SiHa cells were transfected using Lipofectamine 3000 reagent (Invitrogen, USA) following the manufacturer's protocol. Transfected cells were incubated for 48 h before further analysis.

Reverse Transcription and Quantitative Real-Time PCR.
Total RNA was extracted using TRIzol (Sangon, Shanghai, China) according to the manufacturer's protocols. The quantity and quality of RNA were measured using the NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). qRT-PCR was performed using EasyScript® First-Strand cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China) and 2x SG Fast qPCR Master Mix (Sangon Biotech, China) with the Roche LightCycler480 PCR system. The 2 -ΔΔCT method was used to calculate relative expression. The primer sequences are listed in Table 1. Experiments were performed in triplicate.

Transcriptome Sequencing of mRNAs and Differential
Expression Analysis. Sequencing libraries were produced using the Hieff NGS™ MaxUp Dual-mode mRNA Library Prep Kit for Illumina (YEASEN, Shanghai, China) based on the instructions from the manufacturer. First, total RNA was quantified using the Qubit 2.0 RNA Assay Kit (Life, USA) to determine the amount of total RNA added to the library construction. Purification and fragmentation of mRNA were performed using the mRNA Capture Beads and Frag/Prime buffer. Synthesis and purification of double-stranded DNA were conducted using 1st Strand Synthesis Reaction Buffer, 2nd Strand Synthesis Buffer, and Hieff NGS™ DNA Selection Beads before PCR. Amplification and Table 1: Primer sequences used for qRT-PCR amplification.

euKaryotic Ortholog Groups (KOG) Functional
Annotation and Gene Set Enrichment Analysis (GSEA). KOG is the NCBI annotation system based on the direct homologous relationship of genes, in which Clusters of Orthologous Groups of proteins (COG) are targeted at pro-karyotes for eukaryotes. KOG currently has 4852 classifications and can be combined with evolutionary relationships to classify homologous genes from different species into different ortholog clusters.
To explore potential mechanisms underlying the effect of BRCA1 in CCRT resistance in CSCC, GSEA was performed to detect whether a priori defined set of genes showed statistically significant differential expression between the BRCA1 overexpression group and negative control group using Broad Institute software [28]. Gene sets with a nominal P value < 0.05 and FDR < 0:25 were considered significantly enriched.
Transcriptional Regulatory Relationships Unraveled by Sentenced-based Text mining Version 2.0 database (TRRUST) (https://www.grnpedia.org/trrust/) is a database based on the literature to reflect the relationship between transcriptional

Processing of Microarray Data and Identification of DEGs
from the GEO Dataset. The GSE56363 dataset was normalized and the results are shown in Figure 1(a). A total of 21 CSCC patients at advanced stage including 9 NCR and 12 CR patients were analyzed. We obtained 609 DEGs including 223 upregulated DEGs and 386 downregulated DEGs in NCR patients compared with CR patients (Figure 1(b)).

GO Enrichment Analysis for DEGs.
To obtain the function of the DEGs related to CCRT in CSCC, enrichment analysis of DEGs was performed using clusterProfiler in R. As shown in Table 2 and Figure 2, in the biological process group, the DEGs were most enriched in Response to virus. In the cellular component group, the DEGs were most enriched in the Integral component of lumenal side of endoplasmic reticulum membrane. In the molecular function group, the DEGs were most enriched in extracellular matrix structural constituent.

Construction of the PPI Network and Screening and
Analysis of Hub Genes. The PPI of the DEGs was collected by STRING and visualized in Cytoscape software. We evaluated the degree and betweenness centrality in the PPI network through the CytoHubba plug-in and identified the top 10 hub genes with the highest degrees: BRCA1, CDCA8, ASPM, CDC45, RAD51, HMMR, CENPF, EXO1, DTL, and ZWINT ( Figure 3). Among them, BRCA1 has the highest degrees. The ROC curves of the hub genes all indicated favorable diagnostic value in CSCC patients with CCRT resistance (Figure 4).

Validation of BRCA1 Expression. The qPCR results
showed that the expression level of BRCA1 mRNA was significantly upregulated in SiHa/DDP cells compared with SiHa cells (Figure 5(a)). Furthermore, the expression level of BRCA1 was significantly increased in SiHa cell lines transfected with the GV230-BRCA1 plasmid (OV-BRCA1 group) compared with cells transfected with empty vector (NC group) ( Figure 5(b)).  After removing reads with adapters, unknown nucleotides, and low-quality reads of sequencing, we acquired a total of 15.79 Gb clean data. The Q30 percentage of all samples was more than 90% (Table S1). The transcriptome sequencing data have been uploaded to the GEO database, and the accession number is GSE141558. We normalized the read count data using the FPKM method. Using FDR < 0:05 and |LogFC | >1 as the screening criterion in the process of gene identification, the transcriptome sequencing analysis identified 1344 DEGs between the OV-BRCA1 and NC groups, including 824 upregulated DEGs and 520 downregulated DEGs (Table S2). The enrichment analysis results showed that in the biological process group, the DEGs were most enriched in response to virus; in the cellular component group, the DEGs were most enriched in the nucleosome; and in the molecular function group, the DEGs were most enriched in transferase activity transferring pentosyl groups Table 3 and Figure 6).
We also performed KOG annotation to reveal the function distribution of the DEGs to clarify the embodiment of sample differences in gene function. The results of the DEGs classified on the KOG system are shown in Figure 7.

BRCA1 May Affect CCRT Resistance via the JAK/STAT
Signaling Pathway. To more closely examine the mechanisms of BRCA1 in CCRT resistance in CSCC, we intersected the key genes in the WNT and JAK/STAT signaling pathway and target genes of BRCA1 based on the TRRUST database (Figure 9(a)). We found that BRCA1 may upregulate the expressions of STAT1, STAT2, and CCND1, three key genes in the JAK/STAT and WNT signaling pathway (Figure 9(b)). We performed qPCR verification, and the results showed that only STAT1 was significantly increased in SiHa cells with GV230-BRCA1 plasmid transfection (Figure 9(c)). These findings showed that overexpression of BRCA1 in CSCC cells may enhance resistance to CCRT through activating the JAK/STAT signaling pathway via upregulating STAT1.

Discussion
CCRT is the standard therapy for locally advanced CSCC; however, the development of intrinsic or acquired CCRT resistance dramatically decreases its clinical effectiveness. To investigate efficient biomarkers that could predict the response to CCRT and serve as potential therapeutic targets for CSCC patients, we performed comprehensive analysis of resistance in CSCC based on a GEO dataset and transcriptome sequencing. First, we reanalyzed the GSE56363 dataset and identified 609 DEGs including 223 upregulated DEGs and 386 downregulated DEGs between NCR and CR patients [18]. We performed GO enrichment for DEGs related to resistance in CSCC and identified 10 hub genes    we performed GSEA to identify CCRT resistance related signaling pathways to investigate potential molecular mechanisms of resistance in CSCC. STAT1, STAT2, and CCND1 were selected as the differentially expressed target genes of BRCA1 in the WNT and JAK/STAT signaling pathway and may correlate with the CCRT resistance of CSCC.     BioMed Research International BRCA1 (BRCA1 DNA repair associated) is a tumor suppressor gene and related to hereditary breast and ovarian cancer. BRCA1 plays a pivotal role in homologous recombination-based DNA repair and functions in both cell cycle checkpoint control and maintenance of genomic stability [30]. Genomic instability caused by DNA repair defect after BRCA1 deletion may lead to tumorigenesis [31,32]; however, it was also because of DNA repair defect in BRCA1-deficient cells, DNA damage caused by platinumbased chemotherapy or radiation may be substantially enhanced [33][34][35]. Gao et al. found that advanced or metastatic esophageal squamous cell carcinoma patients with low BRCA1 mRNA expression had increased response rate to cisplatin-based chemotherapy or chemoradiotherapy [36]. For breast cancer and ovarian cancer patients, BRCA1 mutations also increased chemosensitivity and/or radiosensitivity [37,38]. Our study showed that the expression level of BRCA1 was higher in CSCC resistant to CCRT than that in CSCC sensitive to CCRT by bioinformatics analysis and qPCR verification. ROC curves showed that BRCA1 showed a high diagnostic value in CSCC patients with CR to CCRT, which demonstrated that BRCA1 might be correlated with the sensitivity of CC cells to CCRT.
Through the analysis of our sequencing results, we found that STAT1, STAT2, and CCND1 were upregulated in SiHa cells with BRCA1 overexpression, and these findings might reveal part of the mechanisms by which BRCA1 affects sensitivity to chemoradiotherapy in CSCC. Previous studies found that BRCA1 is required for the upregulation of STAT1 and STAT2, and they play a synergistic role in the differential    13 BioMed Research International regulation of IFN-gamma target genes [39,40]. In breast cancer patients who are BRCA1 mutation carriers, CCND1 expression was significantly reduced [41]. These findings are consistent with our results, demonstrating that BRCA1 expression is positively correlated with STAT1, STAT2, and CCND1. STAT1 and STAT2 are important components of the JAK/STAT pathway, and continued activation of the JAK/STAT signal is associated with tumor progression [42,43]. CCND1 is often overexpressed in tumors, and its expression is positively correlated with activation of the JAK/STAT pathway [44]. In addition, studies showed the JAK/STAT pathway played a crucial role in cisplatin resistance of CC [45,46]. Furthermore, our qPCR validation revealed that STAT1 was upregulated in SiHa cells overexpressing BRCA1 in our study. We speculated that BRCA1 may increase the sen-sitivity of CSCC patients to cisplatin-based CCRT by upregulating STAT1 to activate the JAK/STAT pathway.

Conclusion
Our work analyzed transcriptome-related changes in CCRTresistant CSCC patients and SiHa cells based on GEO dataset analysis and transcriptome sequencing. We found that BRCA1 may enhance the sensitivity of CSCC patients to cisplatin-based CCRT by upregulating STAT1 to activate the JAK/STAT pathway. These findings may provide new therapeutic strategies to overcome intrinsic or acquired CCRT resistance in CSCC. Additional studies are needed to further elucidate the specific mechanism of BRCA1 in CCRT resistance of CSCC.