Received 2020-11-03

Revised 2020-12-05

Accepted 2020-12-15

 Correspondence to: Nasrollah Saleh-Gohari,Department of Medical Genetics, Afzalipour School of Medicine, Kerman University of Medical Sciences, Kerman, Iran. Telephone Number: 009834 31325829 Email Address: n_salehgohari@kmu.ac.ir

GMJ.2021;10:e2064

www.gmj.ir

Introduction

Breast cancer (BC), a frequent malignant tumor, emerges as the secondary highest mortality among all tumors in women worldwide [1]. Although previous studies have exhibited that abnormal expression of many genes are involved in the pathobiology, detection, development, invasion, and metastasis of BC, it remains a challenging topic yet [2, 3]. Therefore, it is critical to recognize the underlying mechanisms associated with initiation and progression of BC as well as uncover effective therapeutic targets for BC patients. Latterly, studies have shown that various types of non-coding RNAs (ncRNAs) play key roles in BC development and progression [4].

Much of the human genome (at least 75%) could create active and functional products that only 2% of them are protein-coding transcripts and rest of them are ncRNAs, made up of long noncoding RNAs (lncRNAs), microRNAs (miRNAs), and circular RNAs (circRNAs) [5].

CircRNAs, covalently closed single-stranded RNAs, are ncRNAs discovered recently, which are originated from precursor mRNA (pre-mRNA) via back splicing events [6]. CircRNAs are more resistant than their linear forms to digestion by exonucleases or RNase R, due to absence of 5′ caps and 3′ poly (A) tails [7].

The growing amount of data indicated that circRNAs could take part in various disease including cancers through different mechanisms, such as regulating the tumor progression by post-transcriptional mechanisms, bind to specific proteins, act as transcriptional regulators, translation, and alternative splicing [8, 9]. In recent years, circRNAs have got a lot of attentions due to acting as miRNA sponges and consequently result in upregulating the expression of miRNA-related target genes as well as the development and progression of diverse cancers [10].

In the case of BC, several reports have revealed the essential acts of circRNAs in various processes, such as proliferation, migration, invasion, apoptosis, and give a new perspective for the diagnosis and treatment of BC [11].

hsa_circ_0013958 is located at 1q14 on the plus strand with 816 nucleotides in length and is derived from the lysophosphatidic acid phosphatase 6 (ACP6) gene [12]. It is demonstrated that hsa_circ_0013958/miR-134/CCND1 axis developed cell proliferation and invasion and impeded cell apoptosis in lung adenocarcinoma (LAC) [13].

hsa_circ_0013958 has been recently identified to be upregulated in ovarian cancer (OC) and could act as an oncogene in this cancer [14]. hsa_circ_0003028 (also known as cicFUT8, hsa_circRNA_101368) was introduced as an inhibitor of the migration and invasion in bladder cancer through miR-570-3p sponging and subsequently promoted the expression of Krupple-like Factor 10 (KLF10) [15]. This circRNA originated from exon 3 of FUT8 gene located at chr14:66028054 [12].

It is suggested that hsa_circ_0003028 acted as an oncogene by modulating HMGB1/RAGE signaling pathway in hepatocellular carcinoma (HCC) [16]. Since miR-134 and its target CCND1, and miR-570-3p and its targets HMGB1/RAGE are implicated in the development of BC and on the other hand hsa_circ_0013958 and hsa_circ_000308 could sponge off these miRNAs respectively [17-20], these two circRNAs might play a role in BC.

However, to our knowledge, hsa_circ_0013958 and hsa_circ_0003028 expression and their roles in BC is unknown yet.

Hence, in the present study, we compared the expression of hsa_circ_0013958 and hsa_circ_0003028 between paired breast tumor and tumor's adjacent normal tissue samples as well as evaluation of their expression with demographic and clinicopathological features of BC patients.

In addition, since there is evidence that circRNAs can regulate expression of their target mRNAs of the different coding genes through sponging of their regulatory miRNAs in competing endogenous RNA (ceRNA) networks as circRNA-miRNA-mRNA axes, and on the other hand, deregulation of these networks have been repeatedly reported in tumorigenesis [21].

We analyzed and constructed a circRNA-miRNA-mRNA network for these two circRNAs to find potential ceRNA roles in BC.

Materials and Methods

Study population

In this study, 50 tumor tissues and paired adjacent non-cancerous tissues were collected from Shahid Faghihi hospital, Shiraz, Iran between 2017 and 2019. All BC patients who engaged in the present study did not receive radiotherapy and chemotherapy before surgery. Fresh tissue samples were promptly frozen in liquid nitrogen and preserved in -80°C for further use. The study protocol conforms to the ethical guidelines of the 1975 Declaration of Helsinki and was approved by the local ethical committee at the Kerman University of Medical Sciences (KMU) (assigned ethical code: IR.KMU.REC.1398.565). All the contributors signed a written informed consent accepting to use their clinical information and breast samples. In addition, a questionnaire was used to get the demographic and reproductive data together. Clinicopathological and demographic features of the patients are listed in Table-1 and Table-2.

RNA extraction and cDNA synthesis

Total RNA was extracted from fresh tumors and adjacent normal tissues by using the TRIzol isolation reagent (Invitrogen, Thermo Fisher) according to the manufacturer’s procedures. The quality and integrity of RNA was evaluated by gel electrophoresis and subsequently, the quantity and concentration of extracted RNA was measured by spectrophotometer. Then, complementary DNA (cDNA) synthesis was done using Thermo ScientificTM First Strand cDNA Synthesis Kit (Fermentas, Cat.No: K1622) in agreement with the manufacturer’s instructions.

Real Time PCR

Real-time PCR was executed using RealQ Plus 2x Master Mix Green with High ROzTM (Ampliqon, Cat.No: A325402-25). β2-microglobulin (B2M) gene was applied as an internal control to normalize the data. circRNAs sequences were obtained from CircInteractome. The divergent primers included hsa_circ_0013958: GTCAGAAAGAAGGTAGAGTGGAAC (forward) and CAGGGTGGTCTCATGGTATTG (reverse); hsa_circ_0003028: GAATCTCTCCGCATGTAGAGC (forward) and CAGGTGAATAGACTTCTGTTGTTTC (reverse); and housekeeping gene (B2M) primers were AGATGAGTATGCCTGCCGTG (forward) and GCGGCATCTTCAAACCTCCA (reverse). 1 µl of cDNA, 0.75 µl of each primer, 5 µl DNase-free dH2O and 7.5 µl master mix were used to perform a 15 µl volume reaction, under thermal-cycling conditions involving 45 cycles of 95°C for 20 seconds, and then 60°C for 30 seconds. All reactions were done in duplicate. The values were measured by 2-ΔΔCT (fold change) method.

Statistical Analysis

The data were presented by mean and standard deviation for ΔΔCt, and median for fold change. The comparison of ΔΔCt and fold change between tumors and normal adjacent tissues were made by Paired sample t-test and Wilcoxon test. Mann–Whitney and Kruskal-Wallis tests were applied to assess the association between the expression of circRNAs with clinicopathologic and demographic factors. According to median, fold changes were divided into 2 groups of high and low expressions, and the comparison analyses among these groups were conducted by chi-square test and independent t test. P-value fewer than 0.05 was regarded as a statistical significance. Statistical analyses were performed in IBM SPSS 26 software.

Bioinformatics analysis Investigation of miRNAs and mRNAs related to circRNAs

CircRNA-miRNA interaction was downloaded from database Circinteractome (https://circinteractome.nia.nih.gov). DIANA-miRPath (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=mirpath) and TargetScan (http://www.targetscan.org/vert_72/) databases were used to predict the interaction between microRNAs and mRNAs. But previous data didn’t specify the technique used for its predictions, so, we used mirTarBase (http://miRTarBase.mbc.nctu.edu.tw/) database that had validation methods and selected the targets related to miRNAs with strong evidence. In the last update mirTarBase, MTIs (miRNA–target interactions) validated by experimental data. To date, the database has collected >13404 validated MTIs from 11021 articles by adopting a scoring system.

Construction of circRNA-miRNA-mRNA Network

CeRNA network was constructed based on circRNA-miRNA pairs, miRNA-mRNA pairs, and mRNA-mRNA pairs. PPI (protein-protein interaction) network of mRNAs obtained from database STRING (https://string-db.org). Then, the networks were visualized using Cytoscape 3.7.2 software (https://cytoscape.org/).

Investigation of pathways Correlated mRNAs

In order to explore the molecular pathways involved in mRNAs toppgene (https://toppgene.cchmc.org/) database, it was used to predict. The Venn diagram demonstrated the common pathways associated with mRNAs related to hsa_circ_0013958 and hsa_circ_0003028 (Figure-1, Supplementary file-1). Then, the possible pathways associated with BC and their relationship with demographic factors were evaluated. Investigation of gene Gene properties Properties involved Involved in network Network

Transcription factors

Transcription factors are proteins that balance the rate of gene transcription and bind to enhancers and silencers (DNA-regulatory sequences), commonly localized in the 5'-upstream region of target genes [22]. The TF2DNA (http://fiserlab.org/tf2dna_db//search_genes) database is used to check the genes are transcription factors.

Survival analysis

Overall survival (OS) was used to predict prognosis of patients with BC. GEPIA (http://gepia.cancer-pku.cn/) database was used to assess gene survival analysis based on gene expression and uses Log-rank test.

Results

The expression of hsa_circ_0013958 and hsa_circ_0003028 in tumors and adjacent Non-cancerous tissuesTissues

After determining the expression levels of the target circRNAs in both tumor and normal samples, the Wilcoxon signed ranks test was used to compare these circRNAs relative expressions between tumors and adjacent normal tissues. The results showed that hsa_circ_0013958 expression was significantly lower in tumors (Median= 0.163) than paired adjacent normal tissues (Median= 1.057), (P-value<0. 001, Figure-2a, and 2c). In addition, we observed that the hsa_circ_0003028 is overexpressed in tumors (Median= 2.292) compared with adjacent non-cancerous tissues (Median= 0.867) (P-value=0.045, Figure-2 b, and 2d).

The relation between hsa_circ_0013958 expression and demographic and clinicopathological characteristics Characteristics

In another step of our analyses, Mann–Whitney test revealed that hsa_circ_0013958 had a significantly higher expression in patients who did not use hair dye than women used it (P-value=0.026, Figure-3a, Table-1).

In the next phase, patients were divided into 2 groups of high and low expressions according to median for hsa_circ_0013958 (Table-2). High expression level of hsa_circ_0013958 was observed in the group of participants with the age at menarche ≥14 years (P-value=0.041). We did not identify any significant association between this circRNA expression and clinicopathological features.

The relation between hsa_circ_0003028 expression and demographic and clinicopathological characteristics

The hsa_circ_0003028 expression level was meaningfully higher in the participants with the age at FFTP of ≥25 years compared with those with the age at FFTP of <25 years (P-value=0.055, Figure-3b). Furthermore, the elevated expression of hsa_circ_0003028 demonstrated a significant association with antiperspirant use (P-value=0.022, Figure-3c). The details of the analyses are shown in Table-1. Additionally, according to median, fold changes were divided into 2 groups of high and low expressions (Table-2). The high expression of hsa_circ_0003028 was meaningfully related to regular menstruation (P-value=0.018). In addition, Chi-square test demonstrated that hsa_circ_0003028 is overexpressed in the women with higher age at FFTP (P-value=0.034). We did not observe any significant associations between hsa_circ_0003028 expression and clinicopathological variables of the patients.

Potential circRNAs-mediated sponge regulatory network

CeRNA networks were generated using Cytoscape software as circRNA-miRNA, miRNA-mRNA, and PPI (protein-protein interaction). According to the results of in silico investigation and significant nodes in these networks, the potential hsa_circ_0003028/miRNA/mRNA and hsa_circ_0013958/miRNA/mRNA regulatory network were predicted. The hsa_circ_0013958/miRNA/mRNA network was formed based on 1 circRNA, 6 miRNAs, and 13 mRNAs (Figure-4). We observed that the most interactions are between the mRNAs KRAS, AKT1, SNAI2, and miRNAs including hsa-mir-637, hsa-mir-622, and hsa-mir-545 in this network. The vital role of these mRNAs and miRNAs hsa-mir-637, hsa-mir-622 in BC has been previously revealed [23-27]. The hsa_circ_0003028/miRNA/mRNA network also contained 1circRNA, 22 miRNAs, and 78 mRNAs (Figure-5). The mRNAs SMAD4, MUC1, EZH, SP1, XIAP, NOTCH1, CASP3, VEGFA, TERT, FGF2, BCL2L1, NOTCH2, IGF1R, ERBB2, GLI1, CCND1, CDH1, HIF1A, ATM, BRCA1, FOXO1, BRCA2, RELA, and DNMT1 had the most interactions in this network, whose dysregulation in BC have been previously reported. On the other hand, the hsa-mir-498, hsa-mir-421, hsa-mir-330-5p, and hsa-mir-383 showed the most interactions with hsa_circ_0013958 and related mRNAs in the predicted network. Of note, the hsa_circ_0003028/miRNA/mRNA network had 22 circRNA-miRNA pairs, 44 miRNA-mRNA pairs, and 356 PPI pairs, while hsa_circ_0013958/miRNA/mRNA network contained 7 circRNA-miRNA pairs, 13 miRNA-mRNA pairs, and 18 PPI pairs. Interestingly, it was found that these circRNAs-related ceRNA axes regulate the important biological pathways whose dysregulation occur in BC was. In this regard, our gene interaction results revealed three major molecular pathways including, PI3K-Akt signaling pathway, EGFR tyrosine kinase inhibitor resistance, and HIF1 alpha transcription network.

Transcription factor products of the mRNA member of the evaluated circRNA/miRNA/mRNA networks

We found that the mRNAs with having the transcription factor role involved in the predicted regulatory hsa_circ_0013958/miRNA/mRNA network were SNAI2 and SP7. Furthermore, the mRNAs acting transcription factor in the hsa_circ_0003028/miRNA/mRNA network included DNMT1, FOXO1, HIF1A, TWIST1, RELA, GLI1, SP1, IRF1, SMAD4, FOXO4, MECP2, BRCA1, and KLF9.

Survival Analysis

In another bioinformatic evaluation, we performed survival analysis of all genes involving in the two networks. In this way, a significant result observed for PVT1 in the hsa_circ_0003028/miRNA/mRNA network. By GEPIA web server analysis and overall survival of BRCA, it was indicated that the percent overall survival of patients with low level of PVT1 expression was significantly higher in BRCA (Logrank p= 0.04) compared to patients with elevated level of PVT1 expression (Figure-6).

Discussion

Emerging hotspot evidence demonstrates that deregulation of circRNAs play a broad impact on clinical phenotypes of the BC through their multifaceted biological roles. These ncRNAs via acting as transcriptional regulators, miRNA sponges, protein sponge, and translational regulator, behave oncogene or tumor suppressor [10, 8, 9, 10] for controlling the cell growth, proliferation, invasion, migration, and apoptosis of tumor cells in human various cancers, including BC [11, 28].

Additionally, studies have highlighted that due to their high stability, conservation, and tissue-specific expression, circRNAs are potential diagnostic and prognostic RNA markers in oncology. Herein, our results indicated that hsa_circ_0013958 was down-regulated in BC whereas expression of hsa_circ_0003028 was significantly higher in tumors compared to the paired adjacent normal tissues. A number of studies have revealed the key roles of these two circRNAs in regulating the different pathological aspects of various cancers [13, 15,- 16, 14]. In the case of hsa_circ_0013958, our data showed that it could act as a potential tumor suppressor circRNA in BC. However, this finding is inconsistent with results of the previous studies regarding its oncogenic role in other cancers. For example, Zhu et al. showed that hsa_circ_0013958, is upregulated and overactivated in LAC, and through sponging the miR-134 leads to overexpression of the CCND1 oncogene and in turn, promotes cell proliferation and invasion while inhibits cell apoptosis in LAC [13]. Moreover, in OC, hsa_circ_0013958 is overexpressed in malignant tissues as an oncogene and by intensifying the signaling of epithelial-mesenchymal transition (EMT) dictates malignant biological behavior to tumor cells. And mechanistically, its silencing suppressed the proliferation, migration, and invasion of OC cells and inversely promotion of their apoptosis [14].

Interestingly, our bioinformatics analysis indicated that the hsa_circ_0013958 can potentially act as a sponge for some miRNAs (hsa-mir-637, hsa-mir-622, hsa-mir-545) through given ceRNA networks for regulating its target mRNAs including KRAS, AKT1, SNAI2. In addition, there is clear evidence that these mRNAs and miRNAs play key roles in the physiopathology of BC [23-27]. Therefore, such observations point to the role of hsa_circ_0013958 in BC and the tumor-related function of this circRNA in BC may be the opposite of its roles which has been observed in these cancers. And this regard, there is evidence that a circRNA may play different roles in various tumors, due to its different regulatory targets, especially via acting in several ceRNA networks [21].

On the other hand, hsa_circ_0003028 (circFUT8) has been recently identified as a tumor suppressor circRNA in bladder cancer. And functionally, hsa_circ_0003028 suppressed the migration and invasion features of bladder cancer cells through inhibiting oncogenic miR-570-3p via circFUT8/miR-570-3p/KLF10 axis to remove its inhibitory effect on the tumor suppressor gene KLF10 [15].

However, the reported data has shown that this circRNA could act as an oncogene ncRNA in HCC by stimulating the migration of the HCC cells via miR-200a sponging and activating of HMGB1/RAGE signaling [16].

Again, this circRNA appears to play a dual role in inhibiting and stimulating tumor development in various cancers, and like its function in HCC cancer, it may have oncogenic activity in BC.

Bioinformatically, we predicted 24 protein-coding genes (SMAD4, MUC1, EZH, SP1, XIAP, NOTCH1, CASP3, VEGFA, TERT, FGF2, BCL2L1, NOTCH1, IGF1R, ERBB2, GLI1, CCND1, CDH1, HIF1A, ATM, BRCA1, FOXO1, BRCA2, RELA, DNMT1) and 4 miRNAs (hsa-mir-498, hsa-mir-421, hsa-mir-330-5p, hsa-mir-383) in the ceRNA networks with ceRNA activity of hsa_circ_0003028. According to our network analysis, hsa-mir-498 was a microRNA having strong interactions with other RNAs, which has been prior revealed as a potential oncomiR by targeting the tumor suppressor gene PTEN in BC [29].

The hsa-mir-421 acts as a major tumor suppressor miRNA in BC by targeting MTA1 [30]. It is noteworthy that among the above proteins involved in predicted hsa_circ_0003028/miRNA/mRNA ceRNA networks, thirteen members were key transcription factors (TFs) with known roles some of them in BC, including the DNMT1, FOXO1, HIF1A, TWIST1, RELA, GLI1, SP1, IRF1, SMAD4, FOXO4, MECP2, BRCA1, and KLF9. This suggests that the deregulation of hsa_circ_0003028 hsa_circ_0003028 may play a crucial role in the development of BC through defects in various molecular networks.

At another level, the present study showed that a higher expression of the hsa_circ_0013958 is occurred in the tumor of participants with an age at menarche ≥14 years compared with women age of menarche <14. Also, compared to women with no regular menstruation, an elevated expression of hsa_circ_0003028 was significantly observed in women who experienced regular menstruation in their life. It is well-known that reproductive factors such as age at menarche and menstruation status are related to the risk of BC through hormonal mechanisms [31, 32].

Early age at menarche is linked to the high risk of BC, due to increasing the number of ovulatory cycles and earlier exposure to ovarian hormones [33, 34]. Furthermore, there are some reports that the increased risk of BC has been also associated with regular menstrual cycles via hormonal effects [32]. Mechanistically, clinical and epidemiological evidence have strongly shown that ovarian hormones, estrogens, and progesterone, affect normal breast cell proliferation and development, and increased exposure to these hormones leads to greater random genetic errors, and subsequently occurrence of breast malignancy [35].

Interestingly, the results of previous studies have indicated that some of these factors are associated with changes in the expression of ncRNAs in breast tissue [34, 36].

Therefore, these findings may indicate that one of the possible protective mechanisms of older age at menarche against BC in women is an increase in the expression of hsa_circ_0013958, which the present study revealed that its expression decreases in breast tumors.

On the other hand, the mechanism of increased risk of this cancer by regular menstruation may be upregulation of hsa_circ_0003028, which is highly expressed in breast tumors. Another interesting finding of this study was that the hsa_circ_0003028 had a significantly higher expression in the women with age at FFTP of ≥25 years than those with the age at FFTP of <25 years. Based on numerous findings from previous studies, the younger the age of the first pregnancy, the lower the risk of BC, and the later the age of the first pregnancy, the higher the risk of BC in women. FFTP age can affect the risk of BC through hormone-dependent morphological and differentiation changes in breast tissue by decreased estrogen and prolactin [37-39].

Some studies indicated that circRNAs could be involved in reproductive development and endocrine-related pathways, and have potentially regulatory roles during oogenesis [40, 41]. Cheng et al. revealed that circRNA_103827 and circRNA_104816 were overexpressed with maternal aging, and there is a negative correlation between the expression of these circRNAs and embryo numbers. Moreover, circRNA_103827 and circRNA_104816 are involved in ovarian steroidogenesis [42].

Thus, these data suggest that circRNA hsa_circ_0003028 with potential oncogenic functions may be influenced by hormone-related BC risk factors including older FFTP, and in turn, is involved in BC initiation and development possibly by estrogen-related pathways. Nevertheless, the possible link between estrogen-related pathways and hsa_circ_0003028 is still unknown and more studies are needed.

In this study, our results also showed that the expression of hsa_circ_0013958 was significantly higher in patients who did not use hair dye than women used it. Hair dye compounds contain of some mutagenic and endocrine-disrupting chemicals (EDC) such as aromatic amines, 4-aminobiphenyl (4-ABP), and P-Phenylenediamine (PPD), which may play a role in several human cancers [43, 44].

4-ABP is carcinogenic to human, could disrupt the estrogen-related pathways, and increased the risk of BC [45].

EDCs could affect BC risk through interfering with endocrine system and subsequently, changed the mammary gland development [46], causing epigenetic changes that influenced on chromatin integrity and gene expression [47], disrupting the steroid hormones synthesis or metabolism, and acting directly on steroid hormone receptors as either agonists or antagonists [48]. EDC may also affect the expression of ncRNAs in BC; for example, it could upregulate miR-21 expression and downregulate its target genes, PDCD4 and PTEN, in BC cell lines [49], and also lncRNA HOTAIR expression was induced by endocrine-disrupting chemicals in BC [50].

This is the first study which evaluated the relation between a circRNA expression and hair dye use, but further studies are needed to demonstrate the association between endocrine-disrupting chemicals and hsa_circ_0013958 expression. And finally, we observed that the high expression of hsa_circ_0003028 demonstrated a significant association with antiperspirants use in subjects under study. Antiperspirants contain a variety of harmful chemicals, such as aluminium salts, that may be absorbed by the skin and show a genotoxic profile [51]. Aluminium could potentially change the activity of estrogen receptors in BC cells and interfere with the biological action of breast epithelial cells, thereby have a potential role in transforming normal cells to malignant cells and cause BC progression [52, 53]. Furthermore, it has been investigated that aluminium is involved in oxidative stress and DNA damaging processes through interfering with HIF1 transcription factor [54].

Our bioinformatics analysis also revealed that hsa_circ_0003028 is involved in HIF1 alpha transcription factor network. Together, we suggest that antiperspirants, which contain of aluminium, may affect the hsa_circ_0003028 expression due to regulating the HIF1 alpha transcription factor network.

Conclusion

In conclusion, we experimentally indicated that hsa_circ_0013958 was down-expressed and hsa_circ_0003028 was over-expressed in breast tumors. Moreover, we uncovered for the first time the significant association between the expression levels of hsa_circ_0013958 and hsa_circ_0003028 and some estrogen-related risk factors as well as other risk factors in BC patients. And mechanistically, our bioinformatic findings provide some clues to mechanisms of action of these two circRNAs via the circRNA-related competing endogenous RNA regulatory pathways in BC pathogenesis. Although, more researches are needed to investigate the regulatory roles of these circRNAs in BC development.

Conflict of Interest

The authors declare that they have no comping interests.

Table 2. Expression Ievels of hsa_circ_0013958 and hsa_circ_0003028 and Different Demographic and Clinicopathological Features of the BC Patients.

Features		hsa_circ_0013958 level					hsa_circ_0003028 level
		Low		High		P-value	Low		High		P-value
		N	%	N	%		N	%	N	%
Age	<50	14	63.6%	8	36.4%	0.087	9	40.9%	13	59.1%	0.254
	≥50	11	39.3%	17	60.7%		16	57.1%	12	42.9%
Age at menarche	<14	19	61.3%	12	38.7%	0.041*	18	58.1%	13	41.9%	0.145
	≥14	6	31.6%	13	68.4%		7	36.8%	12	63.2%
Regular menstruation	No	4	80.0%	1	20.0%	0.157	5	100.0%	0	0.0%	0.018*
	Yes	21	46.7%	24	53.3%		20	44.4%	25	55.6%
Age at menopause	<50	7	50.0%	7	50.0%	0.486	8	57.1%	6	42.9%	1.000
	≥50	8	38.1%	13	61.9%		12	57.1%	9	42.9%
BMI	≤25	8	40.0%	12	60.0%	0.477	11	55.0%	9	45.0%	0.826
	25-29	12	54.5%	10	45.5%		10	45.5%	12	54.5%
	≥30	5	62.5%	3	37.5%		4	50.0%	4	50.0%
Age at FFTP	<25	15	48.4%	16	51.6%	0.741	18	58.1%	13	41.9%	0.034*
	≥25	7	53.8%	6	46.2%		3	23.1%	10	76.9%
Breastfeeding duration	0	3	50.0%	3	50.0%	0.433	4	66.7%	2	33.3%	0.476
	<24	9	64.3%	5	35.7%		8	57.1%	6	42.9%
	≥24	13	43.3%	17	56.7%		13	43.3%	17	56.7%
Number of abortions	0	19	52.8%	17	47.2%	0.179	17	47.2%	19	52.8%	0.801
	1	4	33.3%	8	66.7%		7	58.3%	5	41.7%
	>1	2	100.0%	0	0.0%		1	50.0%	1	50.0%
Tumor size	<2.5	16	51.6%	15	48.4%	0.771	17	54.8%	14	45.2%	0.382
	≥2.5	9	47.4%	10	52.6%		8	42.1%	11	57.9%
Estrogen receptor	Negative	0	0.0%	1	100.0%	0.312	0	0.0%	1	100.0%	0.312
	Positive	25	51.0%	24	49.0%		25	51.0%	24	49.0%
Progesterone receptor	Negative	1	50.0%	1	50.0%	1.000	1	50.0%	1	50.0%	1.000
	Positive	24	50.0%	24	50.0%		24	50.0%	24	50.0%
HER2	Negative	15	45.5%	18	54.5%	0.370	18	54.5%	15	45.5%	0.370
	Positive	10	58.8%	7	41.2%		7	41.2%	10	58.8%

BMI: Body Mass Index; FFTP: First Full-Term Pregnancy

*P-values< 0.05 is significant.

Figure 1. The 220 pathways regarded the hsa_circ_0003028 and 110 pathways belong to the hsa_circ_0013958. The Venn diagram showing that 39 molecular pathways related to hsa_circ_0003028 and hsa_circ_0013958.

Figure 2. The expression of hsa_circ_0013958 and hsa_circ_0003028 in tumors and adjacent non-cancerous tissues. A: Lower expression of hsa_circ_0013958 in tumor tissues compared with adjacent normal tissues. B: Overexpression of hsa_circ_0003028 in tumor tissues compared with adjacent normal tissues. C,D: Line charts for comparison of hsa_circ_0013958 and hsa_circ_0003028 expression between tumor samples and paired adjacent normal tissues, according to ΔΔCt.

Figure 3. The association of hsa_circ_0013958 and hsa_circ_0003028 expression with demographic characteristics of BC patients. A: The relation between hsa_circ_0013958 expression and hair dye use. B: The relation between hsa_circ_0003028 expression and age at FFTP. C: The relation between hsa_circ_0003028 expression and antiperspirant use.

Figure 4. The potential hsa_circ_0013958-mediated sponge regulatory network based on experimentally validated and bioinformatically predicted interactions. The blue ellipse shows hsa_circ_0013958, the orange octagon nodes and green triangles present miRNAs and mRNAs, respectively.

Figure 5. The potential hsa_circ_0003028-mediated sponge regulatory network based on experimentally validated and bioinformatically predicted interactions. The blue ellipse shows hsa_circ_0003028, the orange rectangle nodes and green diamonds present miRNAs and mRNAs, respectively.

Figure 6 . Overall survival curve of PVT1 in the hsa_circ_0003028/miRNA/mRNA network. The red and blue line indicate the sample with high and low gene expression, respectively. Statistically significant was P≤0.05.

References

1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin. 2019;69(1):7-34.
2. Consortium EP. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57-74.
3. Al-Mansouri LJ, Alokail MS. Molecular basis of breast cancer. Saudi Med. 2006;27(1):9.
4. Piao H-l, Ma L. Non-coding RNAs as regulators of mammary development and breast cancer. J Mammary Gland Biol. 2012;17(1):33-42.
5. Hangauer MJ, Vaughn IW, McManus MT. Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs. PLoS Genet. 2013;9(6):e1003569.
6. Ashwal-Fluss R, Meyer M, Pamudurti NR, Ivanov A, Bartok O, Hanan M et al. circRNA biogenesis competes with pre-mRNA splicing. Mol Cell. 2014;56(1):55-66.
7. Suzuki H, Zuo Y, Wang J, Zhang MQ, Malhotra A, Mayeda A. Characterization of RNase R-digested cellular RNA source that consists of lariat and circular RNAs from pre-mRNA splicing. Nucleic Acids Res. 2006;34(8):e63-e.
8. Han B, Chao J, Yao H. Circular RNA and its mechanisms in disease: from the bench to the clinic. Pharmacol Ther. 2018;187:31-44.
9. Wu J, Qi X, Liu L, Hu X, Liu J, Yang J et al. Emerging epigenetic regulation of circular RNAs in human cancer. Mol Ther Nucleic Acids. 2019;16:589-96.
10. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK et al. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384-8.
11. Wang X, Fang L. Advances in circular RNAs and their roles in breast Cancer. J Exp Clin Cancer Res. 2018;37(1):206.
12. Glažar P, Papavasileiou P, Rajewsky N. circBase: a database for circular RNAs. Rna. 2014;20(11):1666-70.
13. Zhu X, Wang X, Wei S, Chen Y, Chen Y, Fan X et al. hsa_circ_0013958: a circular RNA and potential novel biomarker for lung adenocarcinoma. FEBS J. 2017;284(14):2170-82.
14. Pei C, Wang H, Shi C, Zhang C, Wang M. CircRNA hsa_circ_0013958 may contribute to the development of ovarian cancer by affecting epithelial mesenchymal transition and apoptotic signaling pathways. J Clin Lab Anal. 2020:e23292.
15. He Q, Yan D, Dong W, Bi J, Huang L, Yang M et al. circRNA circFUT8 upregulates Krüpple-like factor 10 to inhibit the metastasis of bladder Cancer via sponging miR-570-3p. Mol Ther Oncolytics. 2020;16:172-87.
16. Li S, Gu H, Huang Y, Peng Q, Zhou R, Yi P et al. Circular RNA 101368/miR-200a axis modulates the migration of hepatocellular carcinoma through HMGB1/RAGE signaling. Cell Cycle. 2018;17(19-20):2349-59.
17. Chen L. Exploring the Role of Circulating MIR-134 In Breast Cancer Recurrence. 2019.
18. Ahlin C, Lundgren C, Embretsén-Varro E, Jirström K, Blomqvist C, Fjällskog M-L. High expression of cyclin D1 is associated to high proliferation rate and increased risk of mortality in women with ER-positive but not in ER-negative breast cancers. Breast Cancer Res Treat. 2017;164(3):667-78.
19. Wang LL, Huang WW, Huang J, Huang RF, Li NN, Hong Y et al. Protective effect of hsa miR 570 3p targeting CD274 on triple negative breast cancer by blocking PI3K/AKT/mTOR signaling pathway. KAOHSIUNG J MED SCI. 2020.
20. Sun S, Zhang W, Cui Z, Chen Q, Xie P, Zhou C et al. High mobility group box-1 and its clinical value in breast cancer. Onco Targets Ther. 2015;8:413.
21. Abdollahzadeh R, Daraei A, Mansoori Y, Sepahvand M, Amoli MM, Tavakkoly Bazzaz J. Competing endogenous RNA (ceRNA) cross talk and language in ceRNA regulatory networks: a new look at hallmarks of breast cancer. J Cell Physiol. 2019;234(7):10080-100.
22. Lambert M, Jambon S, Depauw S, David-Cordonnier M-H. Targeting transcription factors for cancer treatment. Molecules. 2018;23(6):1479.
23. Orlandella FM, Mariniello RM, Mirabelli P, De Stefano AE, Iervolino PLC, Lasorsa VA et al. miR-622 is a novel potential biomarker of breast carcinoma and impairs motility of breast cancer cells through targeting NUAK1 kinase. Br J Cancer. 2020:1-12.
24. Riggio M, Perrone MC, Polo ML, Rodriguez MJ, May M, Abba M et al. AKT1 and AKT2 isoforms play distinct roles during breast cancer progression through the regulation of specific downstream proteins. Sci Rep. 2017;7:44244.
25. Leivonen S-K, Sahlberg KK, Mäkelä R, Due EU, Kallioniemi O, Børresen-Dale A-L et al. High-throughput screens identify microRNAs essential for HER2 positive breast cancer cell growth. Mol Oncol. 2014;8(1):93-104.
26. Zhang Y, Liu J, Wang J. KRAS gene silencing inhibits the activation of PI3K-Akt-mTOR signaling pathway to regulate breast cancer cell epithelial-mesenchymal transition, proliferation and apoptosis. Eur Rev Med Pharmacol Sci. 2020;24(6):3085-96.
27. Wang W, Hind T, Lam BWS, Herr DR. Sphingosine 1–phosphate signaling induces SNAI2 expression to promote cell invasion in breast cancer cells. FASEB J. 2019;33(6):7180-91.
28. Tang Y-Y, Zhao P, Zou T-N, Duan J-J, Zhi R, Yang S-Y et al. Circular RNA hsa_circ_0001982 promotes breast cancer cell carcinogenesis through decreasing miR-143. DNA Cell Biol. 2017;36(11):901-8.
29. Chai C, Wu H, Wang B, Eisenstat DD, Leng RP. MicroRNA-498 promotes proliferation and migration by targeting the tumor suppressor PTEN in breast cancer cells. Carcinogenesis. 2018;39(9):1185-96.
30. Pan Y, Jiao G, Wang C, Yang J, Yang W. MicroRNA-421 inhibits breast cancer metastasis by targeting metastasis associated 1. Biomed Pharmacother. 2016;83:1398-406.
31. Parsa P, Parsa B. Effects of Reproductive Factors on Risk of Breast Cancer: A. Asian Pac J Cancer Prev. 2009;10:545-50.
32. Henderson BE, Ross RK, Judd HL, Krailo MD, Pike MC. Do regular ovulatory cycles increase breast cancer risk? Cancer. 1985;56(5):1206-8.
33. Clavel-Chapelon F. Cumulative number of menstrual cycles and breast cancer risk: results from the E3N cohort study of French women. Cancer Causes Control. 2002;13(9):831-8.
34. Mansoori Y, Tabei MB, Askari A, Izadi P, Daraei A, Bastami M et al. Expression levels of breast cancer related GAS 5 and LSINCT 5 lnc RNA s in cancer free breast tissue: Molecular associations with age at menarche and obesity. Breast J. 2018;24(6):876-82.
35. Henderson BE, Feigelson HS. Hormonal carcinogenesis. Carcinogenesis. 2000;21(3):427-33.
36. Mansoori Y, Zendehbad Z, Askari A, Kouhpayeh A, Tavakkoly Bazzaz J, Nariman Saleh Fam Z et al. Breast cancer linked lncRNA u Eleanor is upregulated in breast of healthy women with lack or short duration of breastfeeding. J Cell Biochem. 2019;120(6):9869-76.
37. Tamakoshi K, Yatsuya H, Wakai K, Suzuki S, Nishio K, Lin Y et al. Impact of menstrual and reproductive factors on breast cancer risk in Japan: results of the JACC study. Cancer Sci. 2005;96(1):57-62.
38. Russo J, Moral R, Balogh GA, Mailo D, Russo IH. The protective role of pregnancy in breast cancer. Breast Cancer Res. 2005;7(3):131.
39. Abdollahzadeh R, Mansoori Y, Azarnezhad A, Daraei A, Paknahad S, Mehrabi S et al. Expression and clinicopathological significance of AOC4P, PRNCR1, and PCAT1 lncRNAs in breast cancer. Pathol Res Pract. 2020;216(10):153131.
40. Quan G, Li J. Circular RNAs: biogenesis, expression and their potential roles in reproduction. J Ovarian Res. 2018;11(1):9.
41. Zhang C, Liu J, Lai M, Li J, Zhan J, Wen Q et al. Circular RNA expression profiling of granulosa cells in women of reproductive age with polycystic ovary syndrome. Arch Gynecol Obstet. 2019;300(2):431-40.
42. Cheng J, Huang J, Yuan S, Zhou S, Yan W, Shen W et al. Circular RNA expression profiling of human granulosa cells during maternal aging reveals novel transcripts associated with assisted reproductive technology outcomes. PLoS One. 2017;12(6):e0177888.
43. Stiel L, Adkins Jackson PB, Clark P, Mitchell E, Montgomery S. A review of hair product use on breast cancer risk in African American women. Cancer Med. 2016;5(3):597-604.
44. Turesky RJ, Freeman JP, Holland RD, Nestorick DM, Miller DW, Ratnasinghe DL et al. Identification of aminobiphenyl derivatives in commercial hair dyes. Chem Res Toxicol. 2003;16(9):1162-73.
45. Hamblen EL, Cronin MT, Schultz TW. Estrogenicity and acute toxicity of selected anilines using a recombinant yeast assay. Chemosphere. 2003;52(7):1173-81.
46. Bergman Å, Heindel JJ, Jobling S, Kidd K, Zoeller TR, Organization WH. State of the science of endocrine disrupting chemicals 2012. World Health Organization; 2013.
47. Knower KC, To SQ, Leung Y-K, Ho S-M, Clyne CD. Endocrine disruption of the epigenome: a breast cancer link. Endocr Relat Cancer. 2014;21(2):T33.
48. Diamanti-Kandarakis E, Bourguignon J-P, Giudice LC, Hauser R, Prins GS, Soto AM et al. Endocrine-disrupting chemicals: an Endocrine Society scientific statement. Endocr Rev. 2009;30(4):293-342.
49. Teng Y, Manavalan TT, Hu C, Medjakovic S, Jungbauer A, Klinge CM. Endocrine disruptors fludioxonil and fenhexamid stimulate miR-21 expression in breast cancer cells. Toxicol Sci. 2013;131(1):71-83.
50. Bhan A, Hussain I, Ansari KI, Bobzean SA, Perrotti LI, Mandal SS. Bisphenol-A and diethylstilbestrol exposure induces the expression of breast cancer associated long noncoding RNA HOTAIR in vitro and in vivo. J Steroid Biochem Mol Biol. 2014;141:160-70.
51. Darbre PD. Aluminium, antiperspirants and breast cancer. J Inorg Biochem. 2005;99(9):1912-9.
52. Miller WR. Estrogen and breast cancer. Chapman & Hall; 1996.
53. Pineau A, Fauconneau B, Sappino A-P, Deloncle R, Guillard O. If exposure to aluminium in antiperspirants presents health risks, its content should be reduced. J Trace Elem Med Biol. 2014;28(2):147-50.
54. Exley C. Aluminium and Alzheimer's Disease: The science that describes the link. Elsevier; 2001.

Supplementary file1. List of common molecular pathways obtained from mRNAs linked with hsa_circ_0013958 and hsa_circ_0003028. The molecular pathways regarded to each circRNA is taken from the toppgene database.

Pathways in cancer Developmental Biology Adaptive Immune System Breast cancer Proteoglycans in cancer HTLV-I infection PI3K-Akt signaling pathway Axon guidance Pancreatic cancer HIF-1-alpha transcription factor network Focal adhesion Endocrine resistance AGE-RAGE signaling pathway in diabetic complications Thyroid hormone signaling pathway FoxO signaling pathway Melanoma EGFR tyrosine kinase inhibitor resistance Prostate cancer Apoptosis Bladder cancer

Colorectal cancer Central carbon metabolism in cancer Chronic myeloid leukemia Regulation of nuclear SMAD2/3 signaling Small cell lung cancer Toxoplasmosis Measles Hepatitis B p53 pathway feedback loops 2 Integrins in angiogenesis Endometrial cancer Renal cell carcinoma Signaling by PTK6 Prolactin signaling pathway Integrin Signaling Pathway AKT phosphorylates targets in the nucleus Telomeres, Telomerase, Cellular Aging, and Immortality Influence of Ras and Rho proteins on G1 to S Transition S phase