In silico identification of MAPK14-related lncRNAs and assessment of their expression in breast cancer samples

Mitogen-activated protein kinase (MAP kinase) pathways participate in regulation of several cellular processes involved in breast carcinogenesis. A number of non-coding RNAs including both microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) regulate or being regulated by MAPKs. We performed an in-silico method for identification of MAPKs with high number of interactions with miRNAs and lncRNAs. Bioinformatics approaches revealed that MAPK14 ranked first among MAPKs. Subsequently, we identified miRNAs and lncRNAs that were predicted to be associated with MAPK14. Finally, we selected four lncRNAs with higher predicted scores (NORAD, HCG11, ZNRD1ASP and TTN-AS1) and assessed their expression in 80 breast cancer tissues and their adjacent non-cancerous tissues (ANCTs). Expressions of HCG11 and ZNRD1ASP were lower in tumoral tissues compared with ANCTs (P values < 0.0001). However, expression levels of MAPK14 and NORAD were not significantly different between breast cancer tissues and ANCTs. A significant association was detected between expression of HCG11 and estrogen receptor (ER) status in a way that tumors with up-regulation of this lncRNA were mostly ER negative (P value = 0.04). Expressions of ZNRD1ASP and HCG11 were associated with menopause age and breast feeding duration respectively (P values = 0.02 and 0.04 respectively). There was a trend towards association between ZNRD1ASP expression and patients’ age of cancer diagnosis. Finally, we detected a trend toward association between expression of NORAD and history of hormone replacement therapy (P value = 0.06). Expression of MAPK14 was significantly higher in grade 1 tumors compared with grade 2 tumors (P value = 0.02). Consequently, the current study provides evidences for association between lncRNA expressions and reproductive factors or tumor features.

availability of MAPK-targeting therapies 2 , identification of regulatory mechanisms of this pathway has practical significance. The interference between ceRNAs via common miRNAs characterizes a new level of gene regulation that participates in the evolution of human malignancies. Such interferences can be anticipated according to the intersection of miRNA-binding sites 7 .
In the present investigation, we aimed at identification of MAPK-related lncRNAs with putative ceRNA function. Through an in silico approach, we detected MAPK14 as the most interacting RNA with miRNAs and lncR-NAs. Consequently, we focused on this gene to identify the lncRNAs with putative interaction with it. Finally, we assessed expression of MAPK14-related lncRNAs in breast cancer samples and adjacent non-cancerous tissues (ANCTs).

Methods
In silico analyses. The total list of MAPK pathway genes were retrieved from HGNC database (https://www. genenames.org/data/genegroup/#!/group/652). The list of miRNAs identified in Homo sapiens was downloaded from Mirtarbase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) and miRNA-mRNA relationship was evaluated using this tool (based on the experimentally validated miRNA-mRNA relationship using Reporter assay and Western blotting techniques). miRNA-mRNA relationships with weak evidences were filtered. From the obtained list of miRNA-mRNA relationship with strong evidence, those associated with MAPK genes were selected. Subsequently, lncBase v2 (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=l-ncbasev2%2Findex-experimental) was used for assessment of miRNA-lncRNA associations. The identified miR-NAs from the previous step were assessed in lncBase v2 and the associated lncRNAs were retrieved. Scores>0.8 was used as the threshold. The miRNA-mRNA relationship was evaluated using Mirtarbase (http://mirtarbase. mbc.nctu.edu.tw/php/index.php) which is tool which reports these interactions based on the experimentally validated miRNA-mRNA relationship using Reporter assay and Western blotting techniques. The previous steps provided the list of lncRNA-miRNA-mRNA triplets to find the lncRNAs with potential sponging activities. Next, Expression Atlas 8 data was used to identify MAPK genes with differential expression in breast cancer tissues vs. normal tissues. Expression of the previously identified lncRNAs has been assessed in Expression Atlas as well. Finally, Co-lncRNA (http://bio-bigdata.hrbmu.edu.cn/Co-LncRNA/) tool was applied to select lncRNAs which co-express with MAPK14 in breast tissues (Fig. 1).
Patients. In the current project, we enrolled 80 female breast cancer patients aged between 36 and 60 (mean ± (SD) age: 49.59 ± 4.74). Malignant tissues and their corresponding ANCTs were obtained during surgery, promptly transferred in liquid Nitrogen to Genetic laboratory for gene expression analyses. All samples were also assessed by a pathologist to verify the diagnosis. Malignant samples included seven invasive lobular carcinomas, one papillary carcinoma, one ductal carcinoma in situ and 71 invasive ductal carcinomas. Patients were recruited from Farmanieh and Sina Hospitals during 2016-2018. All patients signed inform consent forms. The study protocol was approved by the Ethical Committee of Shahid Beheshti University of Medical Sciences. Expression analysis. RNA was extracted from all samples using the Hybrid-R 100 preps (GeneAll, Seoul, South Korea) according to the instructions. RNA samples were treated with DNAse I (Thermo SCIENTIFIC, Vilnius, Lithuania) to eliminate DNA contamination. Afterward, the RNA quantity and quality was assessed and cDNA was made from extracted RNA using Solis BioDyne kit (Estonia). Relative expressions of MAPK14 and the associated lncRNAs were quantified in all samples using RealQ Plus Master Mix Green (AMPLICON, Odense,   Table 2. The potential lncRNA-miRNA-MAPK genes interaction based on our study design. assumed for parameters with 4000 iteration and 2000 burn-outs. The 95% Highest density interval (HDI) was calculated based on the Bayesian approach. The P values were estimated from frequentist methods using quantile regression and mixed effects models. The 'quantreg' , 'ggplot2' , and 'brms' packages were used in R 3.5.2 environment. The association between tumor features and transcript levels of genes was evaluated using Chi-square test or Fisher exact test where appropriate using the Statistical Package for the Social Sciences (SPSS) v.18.0 (SPSS Inc., Chicago, IL). The significance of alteration between mean values of transcripts between discrete groups of patients was appraised using Tukey's honest significance test. The correlation between transcript levels of genes was dignified using the regression model. For all statistical tests, the level of significance was set at P < 0.05.

Results
In silico assays. There were a total 60 MAPK genes in HGNC database. These genes were assessed by Mirtarbase and lncBase to find miRNA and lncRNA associations. Table 2 shows the potential lncRNA-miR-NA-MAPK genes interaction based on our study design. As MAPK14 was found to have the greatest number of interactions with miRNAs and lncRNAs, subsequent steps were performed on this gene. We further listed miRNAs that were predicted to have associations with MAPK14 and listed the associated lncRNAs (Table 3). Co-expression analysis using GEPIA and Co-LncRNA tools revealed that NORAD, HCG11, ZNRD1ASP and TTN-AS1 lncRNAs co-express with MAPK14 in breast tissues. Consequently, we selected these four lncRNAs for expression analysis.
General data of patients. General demographic and clinical features of enrolled patients are summarized in Table 4.
Expression assays. A total of 80 breast cancer samples and 80 ANCTs were assessed. We could not detect expression of TTN-AS1 in any of malignant or non-malignant tissues, so this gene was excluded from further steps. Expression levels of MAPK14 and NORAD were not significantly different between breast cancer tissues and ANCTs. Expressions of HCG11 and ZNRD1ASP were lower in tumoral tissues compared with ANCTs (P values < 0.0001). Figure 2 and Table 5 show the results of expression analysis.
To further verify our results, we used ENCORI/Starbase v2 database to validate our findings in 1104 cancer and 113 normal samples from the TCGA project. Figure 3 shows that both HCG11 and ZNRD1ASP are down-regulated in breast cancer tissues from TCGA database.

Associations between expression levels of genes and patients' features. A significant association
was detected between expression of HCG11 and ER status in a way that tumors with up-regulation of this lncRNA were mostly ER negative (P value = 0.04). Besides, expressions of ZNRD1ASP and HCG11 were associated with menopause age and breast feeding duration respectively (P values = 0.02 and 0.04 respectively). Moreover, there was a trend towards association between ZNRD1ASP expression and patients age of cancer diagnosis in a way that expression of this lncRNA tended to be up-regulated in tumor samples from pre-menopause patients compared with their paired ANCTs (P value = 0.06). Finally, we detected a trend toward association between expression of NORAD and history of hormone replacement therapy (P value = 0.06). Table 6 and Fig. 4 summarize the results of association analysis between expression of genes and patients' data.
We also compared expression of genes among distinct categories of tumor tissues ( Table 7). Expression of MAPK14 was significantly higher in grade 1 tumors compared with grade 2 tumors (P value = 0.02). No other significant difference was detected in expression of genes among distinct categories of tumors.   www.nature.com/scientificreports www.nature.com/scientificreports/ expression data of triple negative breast cancer (TNBC) patients. They performed functional enrichment on the module that was mostly associated with Ki-67 status (Turquoise module). They also established the ceRNA network. Using this model, they have recognized correlation between two mRNAs (RAD51AP1 and TYMS) and overall survival in TNBC. Their results indicated that TNBC-specific mRNA and lncRNAs might form a complex ceRNA network which can be a putative therapeutic target for TNBC 10 . The main difference between this article and our work is inclusion of only a certain type of breast cancer in the mentioned study and assessment of the whole transcriptome.
MAPK14 codes for α subunit of p38 MAPK. This subunit is the prototypic component of the p38 MAPK proteins that has been initially recognized as a tyrosine phosphorylated protein in triggered immune cell macrophages. In addition, MAPK14 regulates production of a number of cytokines including TNF-α 11,12 . Notably, MAPK14 has an essential role in induction of cell migration and epithelial-to-mesenchymal transition (EMT) in breast cancer cells through cooperation with TGF-β 13 . The observed similar levels of MAPK14 between malignant tissues and ANCTs is in line with the previous finding that paracrine messages from tumor cells enhance the expression of nuclear EMT-transcription factors in neighboring fibroblasts leading to over-expression of EMT associated genes in tumor-adjacent tissues 14 . However, some evidences point to a tumor suppressive role of MAPK14 in breast cancer. For instance, the observed enhanced MAPK14 phosphorylation in Wip1-knockout   www.nature.com/scientificreports www.nature.com/scientificreports/ mice has been associated with lower breast tumor formation 15 . On the other hand, treatment of cancer cell lines with a certain MAPK14 inhibitor has diminished tumorigenic potential in animal models of breast cancer 16 . Notably, we detected higher levels of MAPK14 in grade 1 tumors compared with grade 2 tumors. Taken together, one could speculate different roles for MAPK14 in each step of breast tumorigenesis. Such distinct roles have also been proposed for TGF-β (a partner of MAPK14). While in early phases of breast cancer TGF-β suppresses cell cycle transition and enhances cell apoptosis, in late phases, this cytokine is associated with augmented tumor progression, greater cell motility and malignant behavior of tumor cells 17 .
We reported lower expression of HCG11 in tumoral tissues compared with ANCTs. We also detected a significant association between expression of HCG11 and ER status in a way that tumors with up-regulation of this lncRNA were mostly ER negative. Liu et al. have previously shown associations between up-regulation of HCG11 and poor breast cancer outcome. However, they did not report total expression changes between tumoral and non-tumoral tissues. Besides, they reported association between expression of this lncRNA and ER status 18 . Consistent with our results, this lncRNA has been previously shown to be down-regulated in prostate cancer cells and tissues 19 . Forced overexpression of HCG11 in prostate cancer cells has suppressed cell proliferation, invasion and migration, while enhanced cell apoptosis by regulating miR-543 expression. Besides, this lncRNA suppresses PI3K/AKT signaling pathway to inhibit progression of prostate cancer 20 . miR-543 has an inhibitory role on cell proliferation and cell cycle transition in breast cancer through modulation of ERK/MAPK 21 . Thus, the functional role of HCG11 in breast cancer might be mediated through this miRNA.
Moreover, in line with our observation, Zhang et al. have demonstrated HCG11 as an androgen-responsive lncRNA 19 . Moreover, through assessment of NONCODE data, they have detected over-expression of this lncRNA in endocrine-associated tissues such as ovary, breast and prostate, signifying its role in control of tumor evolution in these tissues 19 . Consistent with the proposed role for this lncRNA in endocrine-associated functions, we detected associations between its expression and breast feeding duration. Notably, the ceRNA network depicted by in silico assessments has shown participation of HCG11 in developmental processes, differentiation, gene expression and angiogenesis 19 . Thus, down-regulation of this lncRNA in tumoral tissues might be associated with decreased differentiation state or increased angiogenic potential.
Expression of ZNRD1ASP was lower in tumoral tissues compared with ANCTs. Besides, expression of ZNRD1ASP was associated with menopause age. Moreover, there was a trend towards association between ZNRD1ASP expression and patients' age of cancer diagnosis in a way that expression of this lncRNA tended to be up-regulated in tumor samples from pre-menopause patients compared with their paired ANCTs. This lncRNA is transcribed from the antisense strand of Zinc ribbon domain containing 1 (ZNRD1) and negatively regulates expression of the sense transcript 22 . Previous studies have shown over-expression of ZNRD1ASP in lung cancer 22 . Moreover, single nucleotide polymorphisms (SNPs) within ZNRD1ASP modulate risk of several human cancers 22,23 .
We also reported a trend toward association between expression of NORAD and history of hormone replacement therapy. This lncRNA participates in the construction of a topoisomerase complex which maintains genome stability 24 . Its over-expression in breast cancer has been associated with poor patients' survival 18 . Consistent with our data, Liu et al. did not detect any associations between its expression and ER, PR and HER2 status 18 .
Although in silico studies have shown co-expression of MAPK14 with the selected lncRNAs, we could not detect significant correlations between expression levels of lncRNAs and MAPK14 except for one case. Such lack   www.nature.com/scientificreports www.nature.com/scientificreports/ of cancer leading to an augmented dependence or association presumably similar to what has been called as oncogene-addiction. However, further experiments are needed to verify this speculation.

Conclusion
In brief, in the present study, we introduced an in silico method for identification of MAPK14-related lncRNAs with putative ceRNA role in breast cancer and assessed expression of these lncRNAs in breast cancer tissues and ANCTs. Our data supports associations between expression levels of these lncRNAs and some clinical features. Future studies are needed to elaborate the underlying mechanisms of such observations. The identified interactome comprising of MAPK14 and the 4 lncRNAs might provide new insight about the role of MAPK14 in the breast carcinogenesis and provide therapeutic targets for this cancer. As a future perspective, we can deepen the role of miRNAs in the mentioned network and assess the contribution of the selected miRNAs and their targets in the MAPK14-mediated breast carcinogenesis. Such studies would increase the insights about the regulatory mechanisms among mRNAs, lncRNAs, and miRNAs and identify promising biomarkers for breast cancer detection and treatment. Finally, this work deals with the transcriptome expression profile of MAPK14 and its associated lncRNAs. However, the effect of this interactome of MAPK14 and other interactors at the protein level were not assessed in this study which is a clear limitation of the present work. Ethics approval and consent to participate. The study protocol was approved by the Ethical Committee of Shahid Beheshti University of Medical Sciences. All patients signed informed consent forms. All steps were performed according to ethical guidelines.  Table 7. Comparison of expression levels of genes among distinct categories of tumor tissues.