Inferences From Dysregulated Long Non-Coding RNA-Mediated Competing Endogenous RNAs in Various Chemotherapy Drugs and Evaluation of Drug Response in Breast Cancer

Backgroud: Differences in individual drug response, especially drug resistance, present an obstacle to the treatment of breast cancer (BRCA). Thus, the ability to predict drug response would contribute to developing novel treatment strategies. Accumulating evidence have suggested that tumor molecular proles and drug response data provide opportunities and challenges for the discovery of new molecular characteristics and mechanisms of drug response in BRCA. Methods: In the present study, an integrated pipeline was developed to explore drug response-related long non-coding RNA (lncRNA)-mediated competing endogenous RNAs (ceRNAs) motifs in BRCA. Results: Drug response-specic ceRNAs indicated that lncRNAs play an essential role in various drug treatments for BRCA. Several key drug-resistant and -sensitive dysregulated ceRNAs were identied in Adriamycin, Cytoxan, and Tamoxifen. The interactions in these ceRNAs showed strong correlations in BRCA. Most drug response-related dysregulated ceRNAs were only present in one kind of drug. A number of drug response-related ceRNAs presented diverse dysregulation patterns. We also extracted some key drug response-related lncRNAs, such as HCP5 and FAM182A. These lncRNAs were associated with certain cancer hallmarks and survival in BRCA. Conclusions: Ultimately, understanding the underlying lncRNA-mediated ceRNAs in drug responses will facilitate improved individual reactions to chemotherapy and overall outcomes of BRCA treatment.


Introduction
Among female cancer patients, breast cancer (BRCA) is the most common cancer and the second leading cause of cancer death after lung cancer [1]. The risk factors for breast cancer include certain menstrual and reproductive factors, postmenopausal obesity, use of hormone replacement therapy, family history of the disease, and carrying high-penetrance mutations [2]. Although the incidence rate of BRCA has increased slightly in recent years, the overall death rate continues to decline. This decrease in BRCA mortality could be accelerated by expanding access to high-quality prevention, early detection, and treatment services to all women. At present, surgical resection and chemotherapy are the two major modes of treatment for BRCA [3,4]. Remarkably, different drug responses to BRCA chemotherapeutic drugs lead to diverse prognoses.
Numerous studies have shown that universal resistance to chemotherapeutic drugs is the most signi cant obstacle to successful BRCA treatment [5][6][7]. Current rst-line chemotherapeutic drugs include Cytoxan, Adriamycin, and Tamoxifen. Cytoxan is an alkylating agent, which kills cells nonspeci cally via chemical reactions with DNA and RNA molecules [8]. Adriamycin is a major chemotherapeutic drug that has been widely applied to BRCA patients [9]. Tamoxifen is often used in estrogen receptor alpha-positive BRCA patients, where it competitively blocks E2-binding in order to prevent co-activator-binding pocket formation and inhibit cell proliferation [10]. Although these chemotherapeutic drugs are all effective treatments for BRCA, acquired drug resistance frequently restricts their curative effects and leads to failure in the clinic. However, most existing studies on chemotherapeutic drug resistance focus on coding genes.
Long non-coding RNAs (lncRNAs) are an important type of non-coding RNAthat lack the ability to encode proteins. An increasing number of studies have found that lncRNAs play essential roles in the occurrence and development of diseases, including many kinds of cancers [11][12][13]. Regarding tumor drug resistance, lncRNAs also are emerging as key regulators of chemotherapeutic drug responses [14]. For example, lncRNA MALAT1 modulates the chemoresistance of Docetaxel-resistant lung adenocarcinoma cells [15]. LncRNA HORAS5 mediates castration-resistant prostate cancer survival by activating the androgen receptor transcriptional program [16]. It is believed that lncRNAs participate in the process of drug resistance via competing endogenous RNA (ceRNA) motifs. CeRNAs are transcripts that crossregulate each other by competing for shared microRNAs (miRNAs) [17]. For instance, lncRNA UCA1 promotes ge tinib resistance through a ceRNA that targets FOSL2 by sponging miR-143 in non-small cell lung cancer [18]. LncRNA GSTM3TV2 upregulates LAT2 and OLR1 by competitively sponging let-7 to promote gemcitabine resistance in pancreatic cancer [19]. Thus, systematic identi cation and characterization of the roles of ceRNAs in chemotherapeutic drug responses are essential.
In the present work, an integrated computational approach was developed to identify multiple chemotherapeutic drug response-related ceRNAs. Several differentially expressed lncRNAs, genes, and miRNAs were discovered in Cytoxan-, Adriamycin-, and Tamoxifen-resistant and -sensitive breast cancer groups. Drug-resistant and -sensitive ceRNA networks were constructed for each type of drug in BRCA patients. These networks showed speci c topological features and core modules. The distributions of interactions were similar. Based on these ceRNAs, we identi ed several important lncRNAs that were related to drug response in BRCA. These drug response-related lncRNAs were associated with cancer hallmarks and thus may potentially serve as effective prognostic biomarkers in BRCA. In summary, this study found that lncRNA-mediated, drug response-speci c ceRNAs play essential roles in BRCA.

Materials And Methods
Collection of lncRNA, gene, and miRNA expression pro les in BRCA The expression pro les of whole-genome transcripts, including coding genes, miRNAs, and lncRNAs, were downloaded from The Cancer Genome Atlas (TCGA;https://gdc.xenahubs.net/download/TCGA-BRCA.htseq_fpkm.tsv.gz). The dataset included 1102 BRCA tumor and 390 control samples. All lncRNAs, genes, and miRNAs that were not expressed in all samples were excluded. The minimum value of all samples was given to any remaining expression values of 0. All expression values were transformed using log 2 (value+1). In addition, clinical information, including chemotherapeutic drugs, drug response, disease process, and survival, were also obtained. Considering the number of samples and clinical applications, only the data from BRCA patients treated with Cytoxan, Adriamycin and Tamoxifen were extracted for subsequent analysis. BRCA patients with "complete response" were considered as the drug-sensitive group. BRCA patients with "stable disease,""clinical progressive disease," or "partial response" were considered as the drug-resistant group.
Acquisition of experiment-veri ed miRNA-lncRNA and miRNA-gene interactions to construct ceRNA motifs In order to obtain accurate ceRNA motifs, experiment-veri ed miRNA-lncRNA and miRNA-gene interactions were collected. We downloaded the miRNA-lncRNA interactions from the RNA Association Interaction Database (RAID)v2.0, which contained experimental and computational prediction RNA-RNA and RNA-protein interactions from manual literature readings and other database resources. Only miRNA-lncRNA interactions with strong experimental evidence were extracted [20]. Then, gene-miRNA interactions were acquired from miRTarBase 7.0, which is a public database that includes numerous miRNA-target interactions [21]. Based on these two types of interactions, candidate ceRNA motifs were constructed.
Identi cation of drug response-associated dysregulated ceRNAs by integrating expression and candidate ceRNA interactions in BRCA Drug-resistant and -sensitive dysregulated ceRNAs for Cytoxan, Adriamycin, and Tamoxifen in BRCA were identi ed using an integrated pipeline. First, differentially expressed lncRNAs, miRNAs, and genes were identi ed between drug-resistant and -sensitive groups for diverse kinds of drugs. Second, Pearson correlation coe cients (PCCs) were calculated in the drug-resistant and -sensitive groups based on expression values of the gene-lncRNA, gene-miRNA, and miRNA-lncRNA pairs. Third, we ltered drug response-associated dysregulated ceRNAs from expression datasets based on the following rules: PCC (lncRNA-gene) > 0.2; PCC (lncRNA-miRNA) < -0.2, and PCC (gene-miRNA) < -0.2, where PCC (lncRNA-gene), PCC (lncRNA-miRNA), and PCC (gene-miRNA) represent the Pearson correlation coe cients of lncRNAgene, lncRNA-miRNA, and gene-miRNA interactions, respectively, based on expression values. Finally, if the PCC difference between drug-resistant and -sensitive groups for any one kind of interaction(gene-lncRNA, gene-miRNA, or miRNA-lncRNA) was greater than 0.3, the speci c ceRNA was considered as a drug response-related dysregulated ceRNA in BRCA.
Construction, topological feature analysis, and core module identi cation of drug response-associated dysregulated ceRNA networks Cytoscape software 3.3.0 (http://www.cytoscape.org/)was used to construct and illustrate the dysregulated ceRNA regulatory networks. Node degree was analyzed by the built-in Network Analyzer tool in Cytoscape. The degree for a node corresponds to the number of its network neighbors. Core modules were extracted from each dysregulated drug response-related ceRNA network based on the ClusterOne package in Cytoscape. ClusterOne could identify tightly connected core modules and score them. The number of nodes and scores were considered as reference conditions for screening core modules.
Functional and cancer hallmark analyses for drug response-associated dysregulated ceRNAs in BRCA The Enrich web based tool (http://amp.pharm.mssm.edu/Enrichr/) with default parameters was used to perform functional annotation of lncRNAs in drug response-associated dysregulated ceRNAs in BRCA [22]. Signi cantly enriched Gene Ontology (GO) terms (P< 0.05) were obtained. Cancer hallmark-related GO terms were acquired from a previous study [23]. The genes in these terms were downloaded from GO using AmiGo (version: 2.5.12; http://amigo.geneontology.org/amigo) [24]. All cancer hallmark-related genes were used to obtain interactions with genes in drug response-associated dysregulated ceRNAs.
Survival analysis for drug response-associated dysregulated ceRNAs in BRCA To evaluate the e ciency of drug response-associated dysregulated ceRNAs in predicting the survival of BRCA patients, an integrated risk model was constructed. Univariate Cox regression analysis was performed to demonstrate the associations between survival and expression levels of lncRNAs, genes, and miRNAs in drug response-associated dysregulated ceRNAs. The risk score for each BRCA patient was computed based on the linear combination of expression values weighted by the coe cient from the univariate Cox regression analysis: Where Cox i is the Cox regression coe cient of an lncRNA, gene, or miRNA node in a speci c drug response-associated dysregulated ceRNA, N is the number of nodes in a ceRNA motif, and Exp(i) is the expression value of node i. All BRCA patients were divided into high-and low-risk groups based on the median risk score. The patients with high-and low-risk scores were expected to have poor and good survival outcomes, respectively. Kaplan-Meier survival analysis was performed for the two groups. The log-rank test was used to assess statistical signi cance (P< 0.05). R 3.6.2 statistical software(The R Foundation, Vienna, Austria) was utilized for all analyses.

Results
Differentially expressed key genes between drug-resistant and -sensitive BRCA patients treated with diverse kinds of chemotherapy drugs In order to characterize the molecular features in BRCA patients with diverse drug responses, we identi ed drug-resistant and -sensitive differentially expressed lncRNAs, genes, and miRNAs. In Cytoxan, Adriamycin, and Tamoxifen, the quantity of differentially expressed lncRNAs, genes, and miRNAs was variable ( Figure 1A). We discovered that most of these differentially expressed molecules were downregulated. Some of these lncRNAs, genes, and miRNAs were differentially expressed in multiple kinds of chemotherapy drugs ( Figure 1B). We identi ed 2446, 579, and 110 common drug responserelated lncRNAs, genes, and miRNAs, respectively, in three kinds of drugs. However, several of these common drug response-related lncRNAs, genes, and miRNAs showed different expression patterns ( Figure 1C). For example, 9.73%, 8.63% and 12.73% of molecules showed opposite expression patterns among Cytoxan, Adriamycin, and Tamoxifen, respectively. We also focused on a number of key lncRNAs, which also displayed different dysregulated expression patterns ( Figure 1D). LncRNA CTB-164N12.1 was upregulated in Cytoxan and Adriamycin but downregulated in Tamoxifen. Several other lncRNAs, including LINC01028, CTD-2128A3.2, and RP11-167H9.4,revealedsimilar dysregulated expression patterns. These results indicate that certain lncRNAs, genes, and miRNAs play essential roles in the drug responses of BRCA patients. Furthermore, several lncRNAs showed more complex regulatory patterns in BRCA drug responses.
Speci c features of drug-resistant and -sensitive dysregulated ceRNA networks An integrated computational pipeline was designed to identify drug-resistant and -sensitive dysregulated ceRNAs. We found 1574, 1630, and 2990 drug-resistant dysregulated ceRNAs in Cytoxan, Adriamycin and Tamoxifen, respectively. Furthermore, 48, 44, and 213 drug-sensitive dysregulated ceRNAs were identi ed in Cytoxan, Adriamycin and Tamoxifen, respectively (Figure 2A). In order to better describe drug responserelated dysregulated ceRNAs, drug-resistant and -sensitive dysregulated ceRNA networks were constructed for the aforementioned three drugs. The drug-resistant dysregulated ceRNA network of Tamoxifen had most edges and nodes ( Figure 2B). Three drug-resistant dysregulated ceRNA networks showed scale-free topological features ( Figure 2C and D). Only a small number of nodes in these networks had great degrees. This nding indicates that drug response-related dysregulated ceRNA networks are informative biological networks. We also extracted some core modules from the ceRNA networks of Cytoxan, Adriamycin, and Tamoxifen ( Figure 2E). The core module in the Adriamycin network had 9 nodes (1 lncRNA, 1 miRNA, and 7 genes) and 15 edges. IGF2BP2-AS1 is an lncRNA located in the antisense chain of coding gene IGF2BP2. A previous study reported that IGF2BP2-AS1 is associated with survival in lung squamous cell carcinoma [25]. Accumulating evidence has also suggested a link between dysregulation of IGF2BP2 and cancer [26]. IGF2BP2 has been hypothesized to serve as an N6methyladenosine reader for the promotion of cancer stemness-like properties and pancreatic cancer pathogenesis [27].

Evaluation ofpatterns in drug-resistant and -sensitive dysregulated ceRNAs
We used expression correlations between any two molecules in each dysregulated ceRNA to evaluate the patterns of drug-resistant and -sensitive dysregulated ceRNA networks. We found that the co-expressed correlations between miRNA-lncRNA, miRNA-gene, and gene-lncRNA interactions in both drug-resistant and -sensitive dysregulated ceRNA networks showed similar distributions ( Figure 3A). PCC differences between drug-resistant and -sensitive groups for any one kind of interaction (gene-lncRNA, gene-miRNA,or miRNA-lncRNA) were also analyzed. We discovered that most of these differences were concentrated between 0.2 and 0.6 in the three kinds of drugs ( Figure 3B). In addition, we analyzed the number of dysregulated interactions. Among the three kinds of drug-resistant dysregulated ceRNAs, two dysregulated interactions were most commonly observed. For example, there were 211, 746 and 617 ceRNAs with one, two, and three dysregulated interactions, respectively, in Adriamycin-resistant ceRNAs. In drug-sensitive dysregulated ceRNAs, most ceRNAs only displayed one dysregulated interaction ( Figure  3C, D). These results indicate that our method was accurate and reliable for identifying drug-resistant and -sensitive dysregulated ceRNAs.

Common and speci c drug-resistant and -sensitive dysregulated ceRNAs in BRCA
We also evaluated common and speci c features of drug-resistant and -sensitive dysregulated ceRNAs. More than 90% of drug response-related ceRNAs were only dysregulated in one kind of drug ( Figure 4A). However, a small number of ceRNAs were dysregulated in two or all three kinds of drugs. For example, ceRNA UCP1/miR-130b-3p/RMST were dysregulated in all three drugs ( Figure 4B), and ceRNA FOXO1/miR-15a-5p/CLRN1-AS1 were dysregulated in two kinds of drugs. Although some drug responserelated ceRNAs were dysregulated in multiple kinds of drugs, their dysregulation patterns were diverse. For example, ceRNA UCA1/miR-16-5p/BACE1 were dysregulated in the three kinds of drugs and all displayed three dysregulated interactions ( Figure 4C). More ceRNAs showed different dysregulation patterns than similar ones. Dysregulated interactions were diverse for ceRNA OTX2-AS1/miR-16-5p/TPPP3 in all three kinds of drugs. In addition, we analyzed the frequency of lncRNAs, genes, and miRNAs to identify key molecules ( Figure 4D). For example, key lncRNAs like HCP5, FAM182A, and NPHP3-AS1 were discovered. Newly identi edkey genes included PTEN, HMGA2, and CCND1. We also found that26.28%, 61.05%, and 27.51% of lncRNAs, miRNAs, and genes occurred in three, two, and one kind(s) of drugs, respectively ( Figure 4E). Notably, lncRNA HCP5, which had the highest frequency, was only dysregulated in Tamoxifen ( Figure 4F). A number of other lncRNAs with high frequency were dysregulated in multiple drugs.

Cancer hallmarks in drug-resistant and -sensitive-dysregulated ceRNAs
In order to depict the roles of drug-resistant and -sensitive dysregulated ceRNAs, we performed a functional analysis of lncRNAs. These lncRNAs were enriched in several essential GO terms, including regulation of biosynthetic processes, response to gonadotropin, polyphosphate metabolic processes, and regulation of inositol phosphate biosynthetic processes ( Figure 5A). Some of these GO enrichment terms were related to the development and progression of BRCA. For example, gonadotropin is a glycoprotein hormone that regulates the gonadal development of vertebrates as well as promotes the production and secretion of sex hormones. A large number of studies have reported that BRCA is well known for being strongly in uenced by female steroids [28][29][30]. Inositol phosphate is also thought to play an essential role in BRCA [31]. In addition, we explored the relationships between drug-resistant and -sensitive dysregulated ceRNAs and cancer hallmarks. We found that some cancer hallmarks, such as tissue invasion and metastasis, self-su ciency, insensitivity to antigrowth signals, and evasion of apoptosis, were associated with multiple genes in drug-resistant and -sensitive dysregulated ceRNAs ( Figure 5B). Together, these results indicate that drug-resistant and -sensitive dysregulated ceRNAs have speci c functions in BRCA.
E cacy of drug-resistant and -sensitive dysregulated ceRNAs as prognostic biomarkers in BRCA Finally, we utilized a comprehensive computational method to evaluate the prognosis of drug-resistant and -sensitive dysregulated ceRNAs. Each drug type displayed a different percentage of survival-related drug-resistant and -sensitive dysregulated ceRNAs ( Figure 6A). For example, we observed that 86.36% of drug sensitive-related dysregulated ceRNAs in Cytoxan were related to survival ( Figure 6B). The HR values for genes, lncRNAs, and miRNAs in a given dysregulated ceRNA were diverse ( Figure 6C). These values represent relative protection and risk for molecule survival. We also evaluated the survival of some key drug-resistant and -sensitive dysregulated ceRNAs ( Figure 6D). Dysregulated ceRNA PURA/LINC02120/miR-15a-5p was signi cantly associated with survival (P<0.001). Dysregulated ceRNA MTDH/LINC00092/miR-542-3p was also signi cantly associated with survival (P< 0.001). In addition, BRCA patients with high-risk scores showed poorer prognoses. All of these results suggest that drugresistant and -sensitive dysregulated ceRNAs could serve as effective prognostic biomarkers in BRCA.

Discussion
Chemotherapy is currently the primary and most effective treatment for various cancers, including BRCA, but its e cacy is limited by individual differences in drug response. Drug resistance is one of the notable reasons for BRCA treatment failure. Thus, an understanding of individual drug responses for BRCA patients is urgently needed to improve BRCA treatment. Recent studies have indicated that a number of lncRNAs play essential roles in drug response. In the present work, we identi ed drug response-related lncRNAs based on ceRNA motifs in BRCA. An integrated and computational pipeline was designed to identify drug response-related ceRNAs for multiple drug types according to lncRNA, miRNA, and gene expression levels in BRCA. Thus, we extracted a number of drug-resistant and -sensitive ceRNAs for Adriamycin, Cytoxan, and Tamoxifen in BRCA. Drug response-related ceRNAs displayed diverse dysregulation patterns in multiple drugs. Several key lncRNAs, miRNAs, and genes, such as FAM182A, miR-155-5p, and PTEN, were identi ed based on these dysregulated drug response-related ceRNAs. These ceRNAs were also associated with certain notable cancer hallmarks and thus have the potential to serve as prognostic biomarkers in BRCA.
CeRNAsare an important regulatory mechanism for lncRNAs in cancer [32]. In 2011, Salmena et al. proposed the ceRNA hypothesis on bona de communications between coding and noncoding RNAs. Recent evidence has highlighted the crucial regulatory roles of ceRNA networks in many kinds of cancers [17,33,34]. The functions of ceRNAs in BRCA were also investigated. For example, FOXO1 3'UTR may function as a miRNA inhibitor for modulating metastasis [35], and lncRNA PVT1 is thought to form ceRNA with miR-200 in BRCA [36]. CeRNAs also play speci c roles in BRCA drug resistance. LncRNA NONHSAT101069 was upregulated in BRCA tissues and promoted epirubicin resistance via regulation of the NONHSAT101069/miR-129-5p/Twist1 axis, highlighting its potential as an oncogene and therapeutic biomarker for BRCA [37]. The lncRNA regulator of reprogramming promotes the progression of BRCA and decreases sensitivity to rapamycin through miR-194-3p targeting of MECP2 [38]. In our study, we focused on three kinds of drugs, including Adriamycin, Cytoxan, and Tamoxifen,for BRCA. In these three drug types, several drug-resistant and -sensitive dysregulated ceRNAs were identi ed and analyzed. Thus, our work provides a number of candidate drug response-related ceRNAs in BRCA.
We further analyzed drug-speci c dysregulations for theceRNAs in BRCA. Most of them were only dysregulated in one kind of drug. This result indicates that lncRNAs may exert their diverse roles through ceRNA to in uence multiple types of drug responses in BRCA. Although a small number of drug responserelated ceRNAs were found to be dysregulated in multiple kinds of drugs, they still showed diverse dysregulated patterns. For example, the interactions of drug response-related ceRNAs miR-16-5p/OTX2-AS1/BACE1 were different. In Adriamycin, the interaction between miR-16-5p and OTX2-AS1 was impacted. In Cytoxan, two interactions (miR-16-5p and OTX2-AS1 as well as OTX2-AS1 and TPPP3) were dysregulated. In Tamoxifen, two interactions were also disrupted (miR-16-5p and OTX2-AS1, miR-16-5p and TPPP3). These different dysregulated interactions contribute to drug-speci c responses in BRCA ceRNAs.
In addition, we identi ed several key drug response-related lncRNAs, miRNAs, and genes. For example, a number of key lncRNAs like HCP5, LINC01085, FAM182A, LINC01559, NPHP3-AS1, and FGF13-AS1 were identi ed. Some of these lncRNAs have been reported to be associated with drug response. For example, lncRNA HCP5 promoted fatty acid oxidation through the miR-3619-5p/AMPK/PGC1/CEBPB axis to promote stemness and chemoresistance in gastric cancer, indicating that targeting HCP5 was a novel approach to enhancing chemotherapy e cacy for gastric cancer. In our work, we found that HCP5 was present in 339 drug response-related dysregulated ceRNAs, all of which were in Tamoxifen. Thus, this result indicates that lncRNA HCP5 maybe a key drug response-related lncRNA for Tamoxifen resistance in BRCA.

Conclusion
In the present study, a number of key drug response-related ceRNAs were identi ed and characterized in BRCA. Our work provides novel insights for leveraging publicly available molecular data to evaluate clinical drug responses and contribute to personalized cancer medicine. We propose that lncRNAmediated ceRNAs have the potential to serve as biomarkers for drug resistance, thereby providing an avenue for improved individual treatment of BRCA.

Declarations
Authors' contributions ZQY conceived and designed the experiments. ZZW, LDB, ZH and QQ analyzed data. ZZW wrote this manuscript. All authors read and approved the nal manuscript.    Survival curves for two prognostic drug response-related dysregulated ceRNAs. Blue and yellow represent low-and high-risk groups, respectively.