Integrated High-Throughput Bioinformatics (Microarray, RNA-Seq, and RNA Interaction) and qRT-PCR Investigation of BMPR1B Axis as a Potential Diagnostic Biomarker of Isfahan Breast Cancer.

During the last 50 years, the bulk of cancer research has concentrated on discovering how tumor cells vary from normal cells in gene expression.[ 1 ] The occurrence of abnormality in the gene expression levels could provide various human diseases, including Alzheimer's,[ 2 ] multiple sclerosis,[ 3 ] hepatitis,[ 4 ] diabetes,[ 5 ] and several cancer types, including retinoblastoma,[ 6 ] colorectal cancer,[ 7 ] lung cancer,[ 8 ] head and neck cancer,[ 9 ] hepatocellular carcinoma,[ 10 ] and breast cancer (BC).[ 11 12 13 14 ] An in-depth examination of gene expression patterns in a number of pathogenic conditions, such as BC, can aid in therapy and disease prediction. The leading cause of cancer-related deaths among women globally is BC.[ 15 16 17 18 ] According to Cancer Statistics 2018, BC is the most frequent female malignancy and the main cause of cancer mortality in women, with more than 2.1 million females diagnosed yearly and more than 62,000 fatalities. More than 60% of BC fatalities occur in poorer nations.[ 19 ] Despite improvements in important outcomes for BC patients because of early identification and recent developments in anti-cancer medication, the disease's recurrence rate remains significant.[ 20 21 22 ] Overall, measuring the expression of genes linked to BC, identifying diagnostic and prognostic biomarkers, and understanding gene expression patterns in different clinical and pathological situations linked to BC could provide useful information about the disease and aid in its prevention. Scientists have developed a number of approaches for measuring gene expression. Two solid tools for evaluating gene expression are real-time PCR and microarray. Gene expression profiling and genome-wide gene expression analysis using DNA microarray might offer information on the amount of expression and relative expression of genes and RNAs in various groups, such as “tumor” and “normal,” or “treated” and “untreated.”[ 23 ] The bone morphogenic protein (BMP) receptor family of transmembrane serine/threonine kinases includes bone morphogenetic protein receptor type 1B ( BMPR1B ), which is a member of the transforming growth factor (TGF) superfamily, whose members are dynamically expressed in the endometrium during menstruation, pregnancy, and endometriosis.[ 24 25 ] According to recent research, BMP signaling appears to be a critical regulator during embryonic, cardiac, and brain development.[ 26 ] BMP members and their receptors are essential for the development and proper reproductive function of the ovaries as well as tumor-suppressors against ovarian cancer.[ 27 ] Therefore, in this research, we aimed to find the expression level of BMPR1B in the Isfahan BC samples as a potential prognostic biomarker. Furthermore, demonstrating the biological RNA and protein interactions with the downstream genes, the correlation of the BMPR1B expression with the survival rate of BC patients, and relevant signaling pathways are the main goals of this study.

Our microarray data analysis on three GSE31448 dataset revealed that BMPR1B could have a significant downregulation in BC patients compared to normal individuals and could be considered as a potential biomarker for BC [ Figure 1 ]. P value and adjusted P value were considered as the statistical significance parameters. For both parameters, 0.05 was considered as the significance level. The related plots of microarray data analysis. (a) Volcano plot is indicating the distribution of each gene analyzed in the GSE31448 microarray dataset, based on logFC and –log 10 (adj. P . value). The red color indicates upregulated genes and the green color indicates downregulated genes. BMPR1B is indicating in the plot by a black point as a significantly downregulated gene in the dataset (b) heatmap of top 200 up differentially expressed genes in the dataset After choosing the BMPR1B as the target gene, the KEGG database was used for the analysis of molecular pathways. Pathway enrichment analysis by the KEGG database revealed that BMPR1B is involved in several biological pathways, such as the TGF-beta signaling pathway. Furthermore, the alternative splicing pathway and pathway in cancer are the most relevant pathways to BMPR1B . According to KEGG result, BMPR1B gene family is one of the gene families involved in the pathways leading to cell tumor. Based on the Reactome pathway analysis, BMPR1B presents in signaling by BMP pathway [ Figure 2 ]. Reactome pathway analysis. Based on pathway enrichment analysis, BMPR1B has a significant role in signaling by bone morphogenic protein pathway. The signaling pathway picture provided by Reactome online database Analysis of miRNA–mRNA interaction represents that the most regulatory effect of hsa-miR-181a-5p is on the target genes that have a crucial role in tumorigenesis pathways, such as escape from apoptosis and aberrant reproduction. In addition, BMPR1B is one of these target genes that have an important role in these pathways. Survival analysis based on the GEPIA2 result revealed that the low expression of BMPR1B was correlated to the less survival rate in BC samples. Furthermore, the mRNA–miRNA interaction analysis showed that hsa-miR-181a-5p could have various biological interactions with some other genes, such as RAD9-HUS1-RAD1 interacting nuclear orphan 1 ( RHNO1 ), sterol O-acyltransferase 1 ( SOAT1 ), cell division cycle 5 like ( CDC5L ), UDP-glucuronate decarboxylase 1 ( UXS1 ), and TNFAIP3 interacting protein 1 ( TNIP1 ). This hub gene could be remarkable biomarkers for BC disease [ Figure 3 ]. The mRNA–miRNA interaction analysis. BMPR1B is indicated by a red node in the network The protein–protein interaction analysis was performed using STRING online software. This analysis revealed that BMPR1B could have correlations with some other mRNAs, such as bone morphogenetic protein 2 ( BMP2 ), bone morphogenetic protein 6, SMAD Family Member 4 ( SMAD4 ), SMAD Family Member 5 ( SMAD5 ), and SMAD Family Member 6 (SMAD6). These hub genes could have some important effect on the expression of BMPR1B and hsa-miR-181a-5p in various human diseases, including BC [ Figure 4 ]. The protein interaction analysis result of STRING online software Survival analysis of BMPR1B expression revealed that the high expression of BMPR1B has a not-significant relation with low survival rate of patients [ Figure 5 ]. Survival analysis of BMPR1B , based on GEPIA2 online database. The x axis of this plot represents the time (month) and the y axis represents the precent survival of the patients. Red color indicates the patients with high expression level and the blue color indicates the patients with low expression level. According to this analysis, there is a non-significant correlation between the low expression of BMPR1B and low survival rate of the patients BMPR1B is downregulated in BC The expression of BMPR1B was evaluated by GEPIA2 database and real-time PCR experiment. GEPIA2 [ Figure 6 ] and our experiment [ Figure 7a ] revealed that BMPR1B had a significantly low expression in BC patients (logFC: –2.335, P value = 0.0062). Furthermore, ROC (Receiver operating characteristic) analysis was performed on the expression data of BMPR1B . The result of ROC analysis revealed that BMPR1B could be a suitable diagnostic biomarker for BC patients (AUC (Area Under Curve) = 0.7, P value = 0.003, Figure 7b ). Expression analysis of BMPR1B by GEPIA2 revealed that BMPR1B has a significant low expression in the breast cancer samples, compared to control qRT-PCR analysis of BMPR1B expression level. (a) Relative expression analysis of BMPR1B revealed that BMPR1B has a significant decrease in the breast cancer (BC) tissue's expression (b) ROC (Receiver operating characteristic) analysis revealed that BMPR1B could be a suitable prognostic biomarker for BC Clinicopathological analysis Analysis of clinicopathological features of BC patients was performed on expression data of BMPR1B . Different stages, age, tumor size, lymph node metastasis, menopausal, ER (Estrogen receptor), PR (Progesterone receptor), and Her2 receptor status were analyzed, based on the expression level of gene and miRNA. Our analysis revealed that the BMPR1B expression did not correlate with different pathological statuses in BC.

To design a reasonable experiment and have a logical hypothesis, we performed an integrated high-throughput bioinformatics analysis to obtain BC's hub potential prognostic and diagnostic biomarker. Furthermore, a network analysis was performed to find the correlated transcript factors and the possible effect of coding and non-coding RNAs on each other. Finally, we evaluated the expression level of BMPR1B as the potential prognostic biomarkers that could have a relative regulatory effect on each one by binding to the has-miR-181a-5p in an RNA interaction network. The correlation of the BMPR1B expression with the BC patients’ survival rate was also evaluated. In summary, we demonstrated that BMPR1B could be a significant prognostic biomarker in the BC samples of the Isfahan population, with a decrease in the expression level. ROC analysis revealed that BMPR1B could be considered as a BC diagnostic biomarker in the population. Furthermore, the expression level of BMPR1B had a non-significant positive correlation with the survival rate of the patients. Based on our bioinformatical analysis, we understand that some other mRNAs could be affected by has-miR-181a-5p, such as ERBIN , SOAT1 , RHNO1 , CDC5L , THBS3 , and TNIP1 . Also, during the mRNA–mRNA interaction analysis, we found that BMPR1B could have a remarkable interaction with some other genes, such as SMAD4 , SMAD6 , SMAD7 , SMAD9 , BMP2 , and ACVR2A . Furthermore, we find that SMAD2 , SMAD3 , SMAD6 , and SMAD7 are involved in the TGF-beta signaling pathway; the same pathway that BMPR1B was involved in it. As a suggestion, future experiments could be performed on the expression level of these genes found. In addition, for a better understanding of the interactions and the effects of these biological factors on the expression level of each other, the correlation of these genes with each other and especially with BMPR1B could be examined. Previous studies have revealed that BMPR1B could be affected by related proteins such as SMAD4 .[ 24 ] Based on this study, luciferase and chip assay represented that SMAD4 could enhance PII promoter activity of the ovine BMPR1B by interacting with the SBE1 motif directly. BMPs, members of the TGF-β superfamily, bind to their receptors, including BMPR1B , resulting in phosphorylation of SMAD1 /5/8, which activates the SMAD -dependent pathway.[ 30 ] Our bioinformatics analysis of the signaling pathway of BMPR1B proved the previous studies about this gene and related pathways. About the role of different polymorphisms of BMPR1B and their role in the development of BC, Zheng et al . 2022[ 31 ] revealed that BC susceptibility was shown to be elevated by the T allele of rs1434536 and rs1970801. Based on the study of Bokobza et al . in 2009,[ 32 ] reducing the BMPR1B expression level has a significant correlation with poor prognosis and proliferation of BC cells. Our result supports mentioned result about the expression level of BMPR1B . About hsa-miR-181a-5p, previous studies revealed that this miRNA could have a remarkable effect on cell proliferation in gastric cancer,[ 29 ] reducing oxidation resistance in osteoarthritis,[ 30 ] papillary thyroid cancer,[ 31 ] and suppressing the proliferation and GC (Gastric cancer) cells migration.[ 32 ] It is highly recommended to perform the luciferase assay experiment to investigate the RNA interaction of hsa-miR-181a-5p and BMPR1B . Also, finding the novel relevant single nucleotide polymorphisms to the risk of BC in the BMPR1B and performing the HRM (High resolution melt) qRT-PCR experiment helps better understand the main mechanisms of BMPR1B as a diagnostic biomarker in BC patients. Pathway-specific studies also may help find the upstream and downstream factors regulated by BMPR1B . The authors certify that they have obtained all appropriate patient consent forms. In the form the patient(s) has/have given his/her/their consent for his/her/their images and other clinical information to be reported in the journal. The patients understand that their names and initials will not be published and due efforts will be made to conceal their identity, but anonymity cannot be guaranteed. Nil. There are no conflicts of interest.

All research methods in this study involving human samples were authorized by the Ethics Committee of Al-Zahra Hospital, Isfahan University of Medical Sciences, and all patients signed written consent forms. The survey protocol was authorized via the ethical committee of Islamic Azad Rasht Branch University (the bioethical code number: IR.IAU.RASHT.REC.1398.055). In case–control research, the BC and surrounding normal breast tissue samples were examined from 50 individuals with BC. None of the patients had radiation or chemotherapy experiences in the past. Tissue samples were washed in distilled water before being immersed in RNA later solution (Invitrogen, USA) and frozen in liquid nitrogen for pathologist review. Table 1 lists the clinicopathological features of BC patients. Clinicopathological features of samples ER, PR: Estrogen receptor, Progesterone receptor BC tissue samples and their normal breast tissue equivalents from the same people were obtained, and their total RNA content was extracted using an RNA extraction kit in accordance with the manufacturer's instructions (GeneAll, Seoul, Korea). The first-strand cDNA synthesis kit (Thermo Fisher Scientific, Waltham, MA, USA) was employed to produce cDNA in accordance with the manufacturer's instructions. To analyze the expression of BMPR1B and GAPDH as the reference gene, the cDNA products were kept at –20°C. Magnetic induction cycler was used to conduct the qRT-PCR experiment (Bio molecular Systems, Australia). The program used the device to multiply cDNA in three different temperatures. Initial 95°C for 15 min, the secondary 95°C for 15 s, a 60°C for 20 s, and a 72°C for 20 s. It should be mentioned that from the secondary 95°C onwards, the number of cycles has been 40. Microarray analysis was performed by R Studio software (4.0.2). The GEOquery,[ 28 ] limma,[ 29 ] ggplot2, and pheatmap packages were used to get the data from the GEO database, statistical analysis of microarray data, and drawing plots and heatmap. Three GSE31448by GPL570 ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array) were analyzed to find the significantly dysregulated genes. Quantile normalization method was performed to normalization of the raw microarray expression data. The Bioconductor and CRAN databases were used for downloading the R packages. About 353 BC samples and four control samples were analyzed. The miRNA interaction analysis was performed by miRWalk 2 online software. The RNA interaction network was visualized by the Cytoscape (3.8.0) software. Pathway enrichment analysis was performed by KEGG and Reactome online databases. Statistical analysis was performed by GraphPad Prism 8 software. The Student's t -test, Kolmogorov–Smirnov test, Welch's t -test, one-way ANOVA, and Tukey's multiple comparisons test were performed on the expression data to get the significance level of the experiment.