Differentially Expressed Genes Induced by Erythropoietin Receptor Overexpression in Rat Mammary Adenocarcinoma RAMA 37-28 Cells

The erythropoietin receptor (EPOR) is a transmembrane type I receptor with an essential role in the proliferation and differentiation of erythroid progenitors. Besides its function during erythropoiesis, EPOR is expressed and has protective effect in various non-hematopoietic tissues, including tumors. Currently, the advantageous aspect of EPOR related to different cellular events is still under scientific investigation. Besides its well-known effect on cell proliferation, apoptosis and differentiation, our integrative functional study revealed its possible associations with metabolic processes, transport of small molecules, signal transduction and tumorigenesis. Comparative transcriptome analysis (RNA-seq) identified 233 differentially expressed genes (DEGs) in EPOR overexpressed RAMA 37-28 cells compared to parental RAMA 37 cells, whereas 145 genes were downregulated and 88 upregulated. Of these, for example, GPC4, RAP2C, STK26, ZFP955A, KIT, GAS6, PTPRF and CXCR4 were downregulated and CDH13, NR0B1, OCM2, GPM6B, TM7SF3, PARVB, VEGFD and STAT5A were upregulated. Surprisingly, two ephrin receptors, EPHA4 and EPHB3, and EFNB1 ligand were found to be upregulated as well. Our study is the first demonstrating robust differentially expressed genes evoked by simple EPOR overexpression without the addition of erythropoietin ligand in a manner which remains to be elucidated.


Introduction
EPOR cytokine receptor occurs in diverse types of cells, mainly erythropoietic progenitors and several types of non-hematopoietic tissues as well as cancer cells. In early erythroblasts, EPOR holds the position of regulator of cell size and cell cycle duration [1], and induces transcriptional reprogramming throughout their maturation [2]. EPOR and its tissue-protective effect was described in nervous system, retina, heart, kidneys, endothelium, muscle and bone tissues [3][4][5][6][7][8][9]. EPO/EPOR axis is also implied in regulation of energy metabolism, obesity and insulin response [10,11]. As a typical member of the type I cytokine receptor superfamily with no intrinsic kinase activity, the EPOR molecule consists of extracellular, transmembrane and intracellular regions. Part of the EPOR intracellular region named cytoplasmic box1 is constitutively associated with JAK2 kinases [12] and after binding of the EPO ligand, a cascade of phosphorylations is activated. Activation of the EPO/EPOR complex triggers signaling cascades JAK2/STAT5, PI3K/AKT, RAS/MAP and PKC signaling pathways [13]. It was reported that for successful signal transmission, EPOR translocates to membrane lipid microdomains along with JAK2, STAT5 and LYN kinase [14]. After these events, the signal translocates to the nucleus, where transcription of EPO responsive target genes is turned on. To date, the overall understanding of EPOR effects on the biological processes in the cells is still missing. Global transcriptome and proteome studies were formerly applied to examine EPO/EPOR effects during embryonic and adult erythropoiesis [15,16]. Lately, phospho-proteome analysis of erythroid progenitors detected 121 novel EPO/EPOR target proteins, including core regulators of metabolic processes aldolase and pyruvate dehydrogenase-α1, and also phosphorylation of 19 novel cytoskeletal targets [17]. Although EPOR signaling is well described in a variety of tissues and cells [13,14] studies dealing with EPOR-induced DEGs in malignant cells are completely missing. EPOR plays an important role in the progression of numerous cancers and can be used as a prognostic marker [18], but the precise mechanism of EPOR action on cancer cells needs to be further investigated. In breast cancer, its expression on the surface of tumor cells is associated with better survival, proliferation and invasiveness [19].
The objective of the present study was to assess the effects of EPOR on the whole transcriptome and proteome changes and interpret these results with regard to cell biological processes. For our purpose, we used rat mammary adenocarcinoma RAMA 37-28 cell line with stable EPOR overexpression, which is a subclone of benign non-invasive mammary rat cell line RAMA 37 [3]. Using RNA-seq and proteomics, a global picture of signaling events in RAMA 37-28 vs. RAMA 37 cells has been obtained, in which 233 DEGs were affected. Subsequently, the accuracy of the expression profile of 13 DEGs from RNA-seq result was verified using qRT-PCR. The main goal of our study was to reveal new candidate genes associated with EPOR expression and to categorize them according to their biological processes.

Differentially Expressed Genes and Validation
We used RNA-seq analysis to understand the molecular changes occurring after the introduction and overexpression of human EPOR in the mammary adenocarcinoma cell line RAMA 37. RNA isolated from both the EPOR overexpressed clone RAMA 37-28 and the parental line RAMA 37 was evaluated for its quality. cDNA libraries were prepared from biological replicates that had an optimal fragment size between 150-300 nt. Sequencing yielded 10 million raw reads per sample of both the RAMA 37 and RAMA 37-28 lines. A total of 15,177 genes were mapped for each sample. Genes with a minimum mean logCPM (counts per million) of 3 were considered differentially expressed and included in the differential expression analysis. Genes with logFC (fold change) in the range above ±1.2 were included in the final list of differentially expressed genes (DEGs). Additionally, we checked the p-value for shortlisted DEGs and removed any gene with p-value > 0.01 (Supplementary Figure S1).
A total of 233 genes were differentially expressed in RAMA 37-28 versus RAMA 37 cells (Supplementary Table S1). Among them, 145 were downregulated and 88 were upregulated. To validate the results obtained from RNA-seq, 13 DEGs were analyzed by qRT-PCR. The correlation of both results was carried out by calculating the Pearson correlation coefficient, the value of which was equal 0.943 (p-value = 1.3 × 10 −7 ) (Figure 1). DEGs were divided based on GO biological processes using the peer-reviewed server Omnalysis (http://lbmi.uvlf.sk/omnalysis.html, accessed on 15 February 2023) ( Figure 2).

DEGs Involved in Erythropoietin Signaling
Three genes are associated with erythropoietin signaling, of which two are upregulated and one is downregulated (GO biological processes "erythropoietin signaling"). EPOR overexpression evoked the upregulation of STAT5A and LYN and downregulation of KIT.

DEGs Involved in Erythropoietin Signaling
Three genes are associated with erythropoietin signaling, of which two are upregulated and one is downregulated (GO biological processes "erythropoietin signaling"). EPOR overexpression evoked the upregulation of STAT5A and LYN and downregulation of KIT.

DEGs Involved in Erythropoietin Signaling
Three genes are associated with erythropoietin signaling, of which two are upregulated and one is downregulated (GO biological processes "erythropoietin signaling"). EPOR overexpression evoked the upregulation of STAT5A and LYN and downregulation of KIT.

DEGs Involved in Proliferation
Forty genes associate with cell proliferation (GO biological processes "cell proliferation"). In this regard, 25 genes were downregulated and 15 genes were upregulated as a result of EPOR overexpression in RAMA 37-28 cells.

DEGs Involved in Proliferation
Forty genes associate with cell proliferation (GO biological processes "cell proliferation"). In this regard, 25 genes were downregulated and 15 genes were upregulated as a result of EPOR overexpression in RAMA 37-28 cells.

DEGs Involved in Apoptosis
Twenty-nine genes were categorized either in the regulation of apoptotic signaling pathway or simply involved in the apoptotic process. Of these, 15 were downregulated and 14 were upregulated (GO biological process "apoptosis"). The most downregulated PTPRF (logFC −5.67), SFRP2 (logFC −5.07), EGLN3 (logFC −2.91), XDH (logFC −2.71) and NR4A1 (logFC −2.45) associate with negative regulation of both intrinsic and extrinsic apoptotic signaling pathways. Moreover, downregulated EGLN3, NR4A1 and SMAD3 can participate also in the inhibition of cysteine-type endopeptidase activity (involved in apoptotic process) and downregulated ANGPT1 and RGCC correlate with diminishing of both epithelial and endothelial apoptosis. In spite of the proapoptotic (tumor suppressor) potential of upregulated genes AKAP12 (logFC 3.36), TGFB2 and BTG2, the group of overexpressed genes CXCL12 (logFC 2.28), ZFP385A and CDKN1A associate with negative regulation of the intrinsic apoptotic signaling pathway in response to DNA damage by p53 class mediator (Figure 3a).

DEGs Involved in Metabolic Processes
Sixty-two genes associate with metabolic processes, of which 38 are downregulated and 24 upregulated (GO biological processes "metabolic processes"). Many metabolic processes could be regulated as a result of EPOR-induced differentially expressed genes in RAMA 37-28 cells. The most significant seems to be altered lipid metabolism, in which downregulated KIT (logFC −7.93), PMP22 and CLN6 and upregulated NR0B1 (logFC 7.70), IGFBP7, THRB, STAT5A, ADM and CYP1B1 associate with steroid metabolic process (GO biological processes "metabolic processes"). While upregulated LYN and RAB38 associate with the regulation of phospholipid and lipid metabolic processes, upregulated THRB could be directly involved in the regulation of triglyceride or cholesterol metabolic processes. In addition, PTGS1, GSTA1, PDPN, DAGLA and CYP1B1 genes associate with icosanoid, unsaturated fatty acid and long-chain fatty acid metabolism. In relation to apoptosis, the following associations between the metabolic processes and DEGs appear to be important. The first one is between the metabolism of reactive oxygen species and both downregulated XDH and SMAD3 and upregulated DDAH1 (logFC 3.67), CDKN1A and CYP1B1 genes. The second is between the metabolism of both xenobiotic and glutathione and the upregulated GSTA1 gene, and the last one is between diminished phosphatidylserine metabolism and downregulated PLSCR1 and SERINC5 genes ( Figure 3b).

Proteomic Analysis and Validation by Western Blot
Proteome analysis of more than 5000 peptides confirmed the altered expression of OCM2 (logFC 7.43) and SNCG (logFC 2.08) genes from the transcriptome analysis, but also revealed new DEGs and/or proteins such as ALDH6 (logFC 5.51), SPRR1B (logFC 2.64), ZYX (logFC 2.28), CTPB (logFC 4.62), HPRT (logFC −5.31) and hypothetical protein MJ0443 (logFC −3.51) (Supplementary Figure S2). To validate the results obtained from proteomic analysis, five representative proteins were analyzed using Western blot. Results obtained from both techniques were consistent (Supplementary Figure S3). EPOR-induced changes were confirmed by using specific siRNA designed against EPOR mRNA. In addition, 48 h silencing of EPOR did not abolish the expression of selected STAT5A and OCM2 genes; however, it significantly reduced their protein levels. On the contrary, EPOR silencing did not increase the expression of downregulated HPRT (Figure 4).

Proteomic Analysis and Validation by Western Blot
Proteome analysis of more than 5000 peptides confirmed the altered expression of OCM2 (logFC 7.43) and SNCG (logFC 2.08) genes from the transcriptome analysis, but also revealed new DEGs and/or proteins such as ALDH6 (logFC 5.51), SPRR1B (logFC 2.64), ZYX (logFC 2.28), CTPB (logFC 4.62), HPRT (logFC −5.31) and hypothetical protein MJ0443 (logFC −3.51) (Supplementary Figure S2). To validate the results obtained from proteomic analysis, five representative proteins were analyzed using Western blot. Results obtained from both techniques were consistent (Supplementary Figure S3). EPOR-induced changes were confirmed by using specific siRNA designed against EPOR mRNA. In addition, 48 h silencing of EPOR did not abolish the expression of selected STAT5A and OCM2 genes; however, it significantly reduced their protein levels. On the contrary, EPOR silencing did not increase the expression of downregulated HPRT (Figure 4).

Discussion
EPOR was discovered in 1989 [20] and has been intensively studied for over 30 years. Activation of EPO-EPOR cascades during erythropoiesis results in transcriptional activation of proliferation and differentiation [2], while in non-hematopoietic tissues, anti-apoptotic and protective mechanism of EPOR signaling is described [21][22][23]. The expression of EPOR has been demonstrated in a panel of 29 tumor cell lines [24] and may be related to increased resistance of cancer cells to various therapies. In the current study, we describe EPOR as a causative factor of extensive expression changes referred to 233 DEGs in RAMA 37-28 cell line, affecting numerous biological functions and resulting in the origin of a new cell identity. The genes obtained from RNA-seq have been processed and clustered into

Discussion
EPOR was discovered in 1989 [20] and has been intensively studied for over 30 years. Activation of EPO-EPOR cascades during erythropoiesis results in transcriptional activation of proliferation and differentiation [2], while in non-hematopoietic tissues, anti-apoptotic and protective mechanism of EPOR signaling is described [21][22][23]. The expression of EPOR has been demonstrated in a panel of 29 tumor cell lines [24] and may be related to increased resistance of cancer cells to various therapies. In the current study, we describe EPOR as a causative factor of extensive expression changes referred to 233 DEGs in RAMA 37-28 cell line, affecting numerous biological functions and resulting in the origin of a new cell identity. The genes obtained from RNA-seq have been processed and clustered into particular biological categories (GO biological processes). In our previous studies, we have shown that RAMA 37-28 cells have lower proliferation capacity compared to RAMA 37 cells, as a result of EPOR overexpression [25,26] (Supplementary Figure S4). Currently, we used xCelligence screening to show faster proliferation of RAMA 37-28 cells under conditions of paclitaxel added, whereas EPOR silencing suppressed such an effect (Supplementary Figure S5).
According to our transcriptome analysis, EPOR overexpression in RAMA 37-28 cell line results in 40 DEGs associated with cell proliferation. Indeed, 12 downregulated genes such as KIT, NPPB, SFRP2, SLURP1, XDH, NR4A1, CLDN1, SMAD3, IFT80, CDH3, ECM1 and RGCC and seven upregulated genes, CDH13, CXCL12, TGFB2, FST, FGFR2, STAT5A, SIX1, are involved in the inhibition of epithelial cell proliferation. In this regard, STAT5A, as a member of JAK2/STAT5 signaling cascade regulated by EPO/EPOR activation, is identified as a chemoresistance inducer and the regulator of ABCB1 transporter protein [27] with the ability to stabilize also heterochromatin structure followed by the inhibition of cell growth [28]. We demonstrate EPOR-induced overexpression of STAT5A, whose direct relationship (without EPO addition) using siRNA against EPOR was confirmed. Upregulated TGFB itself encodes the transforming growth factor beta family of cytokines which functions in proliferation, differentiation, adhesion and migration in many cell types [29]. TGFB receptors, moreover, activate both SMAD-dependent and independent pathways that not only regulate SMAD signaling, but also allow SMAD-independent TGFB responses [30]. The data of Livitsanou et al. [31] demonstrated a novel mechanism of TGFB/SMAD signaling modulation by the small GTPase RHOB and show that TGFB/RHOB signaling cross talk affects the nuclear and cytoplasmic responses to TGFB in opposite ways. The slower proliferation of our RAMA 37-28 cells could therefore be explained by the downregulation of SMAD3 and the upregulation of TGBF2 and RHOB genes. Furthermore, Marlow et al. [32] showed that RHOB has divergent downstream signaling partners, which are dependent on the HDAC isoform that is inhibited. When RHOB upregulates only P21 using a class I HDACi (romidepsin), cells undergo cytostasis. When RHOB upregulates BIM using class II/(I) HDACi (belinostat or vorinostat), apoptosis occurs. Combinatorial synergy with paclitaxel is dependent upon RHOB and BIM while upregulation of RHOB and only P21 blocks synergy. Interestingly, RAMA 37-28 cells which demonstrate paclitaxel resistance compared to RAMA 37 cells [26] reveal overexpression of both RHOB and CDKN1A genes. In search of mechanisms responsible for the proliferative changes of RAMA 37-28 cells, we emphasize the overexpression of IGFBP7 and simultaneous downregulation of NPPB genes. Indeed, it has been reported that IGFBP7, involved in p53-dependent growth suppression of lung and colorectal tumors [33,34], decreased the production of NPPB during viral tumorigenesis [35], which points to the conditional expression of some DEGs. Moreover, overexpressed P21 activated PAK1 may also be involved in the negative regulation of cell proliferation involved in the contact inhibition together with downregulated CLDN-1, whose lower expression is either associated with cancer progression (invasion) or improved survival of cancer patients [36]. Highly downregulated proto-oncogene and receptor tyrosine kinase KIT itself might be another important factor contributing to both neoplastic breast epithelium transformation and inhibition of cell migration which were already demonstrated by Janostiak et al. [37].
Furthermore, we demonstrate the upregulation of ephrin receptors (EPHA4, EPHB3) and ligand (EFNB1) which can associate with a slowdown in RAMA 37-28 proliferation. While EPHA4 signaling promotes mesenchymal-to-epithelial transition during morphogenesis in zebra fish embryo development [38], it was reported that overexpression of EPHB3 enhances cell-cell contacts and inhibits the growth in HT-29 human colon cancer cells [39]. EFNB1 plays important role in cell adhesion in healthy tissues [40], while in tumors, it inhibits cell proliferation and adhesion [41].
Interestingly, overexpression of EPOR is associated with the emergence of a malignant phenotype in RAMA 37-28 cell line [42]. In our previous studies, overexpression of EPOR in RAMA 37-28 breast cancer cells altered the sensitivity of these cells to tamoxifen and paclitaxel via mechanism related to prolonged activation of AKT and ERK1/2 pathways, respectively [25,26]. Acquisition of chemotherapy-resistant cancer phenotype of RAMA 37-28 cells may also associate with the overexpression of widely described cell cycle regulator CDKN1A (P21 Waf1/Cip1 ) [43], which can induce cell cycle arrest, followed by reduction of cancer cell growth rate. Indeed, cytoplasmic P21 phosphorylated by AKT increases cell survival and contributes to taxol resistance in glioblastoma cells [44], cisplatin resistance in testicular and ovarian cancer [45,46], doxorubicin resistance in triple-negative breast cancer cells SUM159 [47] and failure of paclitaxel treatment in human nasal squamous carcinoma RPMI-2650 [48]. In addition, P21 mediates 5-fluorouracil resistance in colorectal cancer cells [49] and is associated with poor response to tamoxifen in MCF7 breast cancer cells [50]. CDKN1A overexpression drives cells to acquire a more aggressive phenotype that is capable of escaping cell block, senescence and apoptosis [43]. Equally important in the resistance of RAMA 37-28 cells to tamoxifen alone is PAK1, whose overexpression and the role in tamoxifen resistance of breast cancer patients has already been described [51]. Besides this, ID4, SEMA3F, FOXD, KCNK5, WNK4 and ECM1 associate with chemoresistance origin and SOX5, ITM2A, SNCG, EPHA4, NTS, FGFR2 and ENPEP associate moreover with breast cancer tumorigenesis. The most significantly upregulated DEGs such as SNCG, CPNE8, GRIA3, SOX5 and CRISP3 can be involved in metastasis and invasivity of cancer cells as well. In this regard, SNCG is overexpressed in human infiltrating breast carcinomas and promotes metastasis [52], while it is undetectable in normal and benign breast tissues [53]. Tian et al. [54] indicated a possible relationship between SNCG and P21 and their relation to radioresistance. We show that EPOR evoked upregulation of SNCG (both mRNA and protein levels) in RAMA 37-28 cells and suppose that the overexpression of both SNCG and P21 genes could play a significant role in the chemoresistance phenotype of RAMA 37-28 cells.
The sensitivity of a tumor cell to therapy can also be affected by changes in the cell's metabolism [55]. Our transcriptome analysis reveals 62 DEGs associated with such a biological property and genes such as XDH, CYP1B1, GSTA1, PLSCR1 and DDAH1 are definitely genes worth mentioning. In addition, proteomic analysis enriched this category with ALDH6 and HPRT1 proteins. The expression of XDH inversely correlates with the expression of cancer stem cell-related genes, such as CD44 or CD133, and their downregulation promotes TGFβ signaling in hepatocellular carcinoma [56]. It is noteworthy that knocking-down or inhibiting XDH results in development and progression of hepatocellular HCC and promotes migration and invasion but not proliferation of these cells, while the changes themselves are dependent on increased activation of TGFB2/SMAD2/3 signaling pathway [56]. It seems that low activity of XDH provides a selective privilege to cancer cells and may also contribute to neoplastic differentiation of RAMA 37-28 cells. However, the main biological effects of TGFβ are inhibition or promotion of proliferation, differentiation, apoptosis, cell dormancy, autophagy and cellular senescence; reduced levels of SMAD3 in RAMA 37-28 cells can suppress the function of TGF-β-induced expression of tumor suppressor genes. So, downregulation of TGFB2 results in the expression of anti-apoptotic proteins BCL2 and BCL-W, and enhanced cancer cell survival, which confers platinum resistance in NSCLC and 5-FU resistance in CRC cells, respectively [57,58]. A similar relationship to resistance to cisplatin and paclitaxel chemotherapy was described in the case of GSTA1 and CYP1B1 genes, respectively [59,60]. Although both genes are overexpressed in RAMA 37-28 cells, the CYP1B1 could be one of the mechanisms explaining the resistance phenotype of RAMA 37-28 cells to paclitaxel [26]. The sensitivity and/or proliferative capacity of tumor cells can also be affected by phospholipid scramblase 1 (PLSCR1) which is involved in several processes including phosphatidylserine exposure on an apoptotic cell surface. Its upregulation promoted cell proliferation, invasion and migration, while its downregulation inhibited these effects [61]. On the contrary, downregulation of PLSCR1 expression, which is also seen in RAMA 37-28 cells, significantly inhibited the proliferation, adhesion, migration and invasion of cancer cells [62]. Another important gene that is involved in cellular metabolism and at the same time plays an important role in the altered phenotype of the tumor cell is DDAH1. The encoded enzyme plays a role in NO generation by regulating cellular concentrations of methylarginines, which in turn inhibits NO synthase activity. While endothelium-derived NO stimulates angiogenesis through the inhibition of apoptosis [63], in cancer cells, the roles of NO are diverse, and might have dual pro-and anti-tumor effects depending on local concentration [64]. Finally, the upregulated ALDH6 enzyme confirmed in RAMA 37-28 cells by proteome analysis may also play a significant role in the resistance [65] of our cells, while the role of our downregulated HPRT1 in relation to the published results [66] is questionable.
It is well known that the characteristic attribute of cancer cells is alteration of cell adhesion and/or extracellular matrix organization (ECM). According to DEGs identified by transcriptome analysis of RAMA 37-28 cells compared to RAMA 37 control, we suppose the loss of cell adhesion ability of RAMA 37-28 cells via plasma membrane adhesion molecules which are closely related to strong downregulation of morphogen and tumor suppressor glypican 4 (GPC4). The effect of GPC4 downregulation was earlier demonstrated in relation to breast tumor progression [67], as well as in healthy tissues, to disruption of epithelial integrity and tight junction organization [68]. GPCs have been shown as well-known and accepted cell surface coreceptors for growth factors such as Wnt/β-catenin, FGF, IGF, VEGF and TGFB and matrix modifying enzymes in many cancer cells, thus being involved in control of tumor growth, metastasis, angiogenesis and ECM [69][70][71]. Moreover, Varma et al. [72] showed that GPC4 mRNA level was downregulated in oxaliplatin-resistant ovarian carcinoma cell line A2780/C10. We demonstrate expression changes also in Frizzledrelated protein 2 (SFRP2), which is secreted and incorporated into ECM of the normal and tumor cells. Its association with the fibronectin-integrin protein complex, promotion of cell adhesion and/or transformation of normal mammary epithelial cells into a tumor was already described earlier [73]. SFRP2 was downregulated in radiotherapy-treated glioma patients, and low SFRP2 expression was correlated with advanced tumor stage and poor prognosis. CRISP/Cas9-mediated SFRP2 knockdown promoted soft agar colony formation, cancer stemness and radioresistance of glioma cells, while increased SFRP2 expression exhibited opposite effects. Moreover, SFRP2 knockdown activated Wnt/β-catenin signaling in glioma cell lines, while overexpression of SFRP2 inhibited Wnt/β-catenin activation [74]. In the case of ECM assembly, downregulated LAMB3, SMAD3, RGCC and GAS6 and upregulated GPM6B were categorized. On the other hand, downregulated MMP13 and CST3 and upregulated TGFB2 and PDPN correlate with ECM disassembly. Interestingly, upregulated FBN1 associates with sequestering of TGFB in ECM. Many more DEGs of RAMA 37-28 cells are related to ECM remodeling and changes of cell adhesion, such as upregulated CYR61, MMP10, NOV, CTGF, TMEM47, CDH13, CTNND2, PAK1 and LOXL1 and downregulated FBLN1, CLDN1, ECM1, TENM2, RGD1561161, COL16a1, COLA3, SGCE, THBS2, PTPRF, MEPE, C1QTNF1, PLEKHA2, PCSK6, FLRT2 and GPM6A.
Not only alterations in cell growth and proliferation but also different morphology in RAMA 37-28 cells compared to parental RAMA 37 ones was observed [25,26]. Although the molecular basis of such difference has not been fully described, actually, 84 genes out of 233 associate with cell differentiation of RAMA 37-28 cells. The most significantly overexpressed gene in this category, both at the level of mRNA and the protein itself, is high-affinity calcium ion-binding protein (OCM2), which belongs to the superfamily of calmodulin proteins. We conclude and show for the first time that the expression of the OCM2 gene which associates with positive regulation of cell morphogenesis involved in differentiation (GO biological processes "cell differentiation") is directly related to the expression of EPOR. This effect was confirmed by siRNA against EPOR. Indeed, high expression of OCM2 was demonstrated earlier in rat cancer cell lines undergoing neoplastic transformation [75,76]. Besides OCM2, 19 DEGs participate directly in the morphogenesis process. Thus, we demonstrate several upregulated DEGs related to columnar or cuboidal epithelial differentiation, namely, P21, SHROOM3 and FGFR2. A large body of evidence suggests that P21 also plays an important role in the differentiation of various normal or malignant cells and tissues, but its effect is dependent on the cell type and the stage of differentiation. Interestingly, P21 has a positive role in differentiation in most studies, partly by inhibiting apoptosis, which promotes cell survival [77]. The STAT5A protein itself plays a key role both in the positive regulation of myeloid cell differentiation and gamma-delta T cell differentiation, as well as in the negative regulation of erythrocyte differentiation, while its role in the differentiation of RAMA 37-28 cells remains unclear.

Preparation of the Library
A quantity of 250 ng of RNA was reverse transcribed to first strand cDNA using oligodT primers and QuantSeq. 3 mRNA-Seq Library Prep Kit (Lexogen, Wien, Austria). RNA template was removed using RNA removal solution (RS buffer) (Lexogen, Wien, Austria) and the second strand was synthesized using random hexamer primer with Illumina-compatible linker sequences at its 5 end. The dsDNA libraries were purified using magnetic beads and then amplified by PCR using specific single indexing i7 primers with the aim of adding adapter sequences important for cluster generation and to generate DNA for quality control and sequencing. PCR Add-on kit for Illumina (Lexogen, Wien, Austria) was used for the determination of 20 PCR cycles required for RAMA cells. Later, magnetic beads of the kit were used for purification of amplified libraries, which were checked for the length of the fragments by fragment analyzer.

Sequencing and Data Analysis
cDNA libraries were sequenced using Illumina NextSeq, single-end 75 bp, with 8 million reads per sample. Fastq files were processed and aligned to reference genome (Rattus norvegicus, mRatBN7.2) using STAR aligner (STAR V 2.5.2b, (https://github.com/ alexdobin/STAR/releases/tag/2.5.2b, accessed on 20 February 2023). The pre-processing included adaptor trimming and removal of recommended initial 10 bases. Reads were counted in STAR V 2.5.2b. to perform differential gene expression analysis edgeR (open source R package version 3.12 was used: https://bioconductor.org/packages/release/ bioc/html/edgeR.html, accessed on 20 February 2023). The low read count with less than 3 CPM (count per million) was filtered out using the filterByExp function of edgeR package. The identification of differentially expressed gene (DEGs) was accomplished by using glmTreat and glmQLFit (quasi-likelihood, QL) functions of edgeR in R package considering log fold change (logFC) values beyond ±1.2 and FDR less than 0.05. The logical relation of DEGs between the challenged RAMA cells (RAMA 37-28 vs. RAMA 37) was calculated using Excel (MS office) and Venn diagrams were constructed. To categorize DEGs groups into GO biological processes and to construct heat maps https://reactome.org/ (accessed on 20 February 2023) and http://www.heatmapper.ca/expression/ (accessed on 20 February 2023), respectively, were used. ECL Western Blotting Substrate (Pierce, Thermo Fisher Scientific, Waltham, MA, USA) and visualised by bio-imaging system (DNR MF Chemibis 2.0, Neve Yamin, Israel). Equal sample loading was verified by immunodetection of anti-β-actin Monoclonal Antibody AC-15 (MA1-91399). The pictures were scanned with GS-800 Calibrated Densitometer and the quantification was performed using Image J software version 1.52 (NIH; National Institutes of Health, Bethesda, MD, USA).

Proteomics Analysis
Total proteins (100 µg) were reduced with 10 mM dithiothreitol in 100 mM ammonium bicarbonate at 56 • C for 30 min. Alkylation was performed with 100 µL of 50 mM iodoacetamide in 100 mM ammonium bicarbonate in the dark for 30 min. Proteins were digested with trypsin in ratio 1:100 at 37 • C, overnight. Aliquots of purified complex peptide mixtures of 100 ng were separated using Acquity M-Class UHPLC (Waters, Etten-Leur, The Netherlands). Samples were loaded onto the nanoEase Symmetry C18 trap column (25 mm length, 180 µm diameter, 5 µm particles size). After 2 min of desalting/concentration by 1% acetonitrile containing 0.1% formic acid at a flow rate 8 µL/min, peptides were introduced to the nanoEase HSS T3 C18 analytical column (100 mm length, 75 µm diameter, 1.8 µm particle size). For the thorough separation, a 90 min gradient of 5-35% acetonitrile with 0.1% formic acid was applied at a flow rate of 300 nL/min. The samples were nanosprayed (3.1 kV capillary voltage) to the quadrupole time-of-flight mass spectrometer Synapt G2-Si with ion mobility option (Waters, Etten-Leur, The Netherlands). Spectra were recorded in a data-independent manner in high definition MSE mode. Ions with 50-2000 m/z were detected in both channels, with a 1 s spectral acquisition scan rate. Spectra were preprocessed with the Compression and Archival Tool 1.0 (Waters, Etten-Leur, The Netherlands) to reduce noise, removing ion counts below 15.

Silencing of EPOR Gene Expression by siRNA
For silence of EPOR expression in RAMA 37-28 cells, we used siRNAs (siRNA Dharma SMARTpool ON-TARGETplus human EPOR siRNA, GE HealthCare, PerkinElmer Holdings Inc., Lafayette, CO, USA), ntRNAs (negative controls to siRNA Dharma ON-TARGETplus-non-targeting) and transfect reagent DharmaFECT 1 (GE HealthCare, PerkinElmer Holdings Inc., Lafayette, CO, USA) following the manufacturer's instructions. Transfection reagents were stable and none of the tested conditions significantly affected cell toxicity. RAMA 37-28 cells were cultured in 6-well plates (500,000 cells/well) in RPMI media without antibiotics, supplemented by 10% of FBS and incubated for 24 h at 37 • C in the presence of 5% CO 2 . Subsequently, cells were cultured in plates with changed fresh media without antibiotics, and were supplemented by 10% of FCS; the experimental groups were enriched by siRNA (14 µM) or ntRNA (5 µM) for the next 48 h. After that, the media was changed and followed by isolation and purification of proteins. Subsequently, Western blot analysis was carried out according to the protocol described above.

Agilent xCELLigence Real-Time Cell Analysis
RAMA 38-28 cells (8 × 10 3 cells/well) were seeded in 96-well plates (RTCA E-Plates 96) in xCELLigence RTCA systems (Agilent Technologies, Santa Clara, CA, USA) for 24 h followed by the treatment using siRNA (14 nM), or ntRNA (5 nM) for 48 h. Subsequently, cells were treated by PTX (200 nM) for 60 h. The cell adhesion and spread of the cells were continuously monitored in 60 min intervals over the course of monitoring period using the xCELLigence RTCA system.

Statistical Analysis of xCELLigence Analysis
Experiments under all conditions were performed in at least three independent measurements. Mean value and standard deviation were calculated using descriptive statistics. The data were analyzed by using the RTCA software Pro 1.2.1 (ACEA Bioscience, Santa Clara, CA, USA). Statistical analysis was carried out by a non-parametric method, one-way ANOVA using SigmaPlot (Ver. 12.0); p < 0.05 was considered significant (see Supplementary Figure S5).

Conclusions
We demonstrate for the first time robust 233 DEGs evoked by human EPOR overexpression in rat mammary adenocarcinoma RAMA 37-28 cells. Identified DEGs are associated with many biological processes such as proliferation, apoptosis, tumorigenesis, cellular metabolism, differentiation and others. It is obvious that EPOR overexpression affects also signaling pathways either by attenuation of RAS signaling, small GTPase, cytokine-mediated transduction, intrinsic apoptotic pathway in response to DNA damage and cAMP-mediated signaling or through the activation of STAT5, NOTCH signaling and certain Ephrin and WNT signaling mediators. We are aware of certain limitations of this work, which is based on the comparison of rat EPOR overexpressed RAMA 37-28 cells and their parental control only. However, we believe that our newly prepared EPORoverexpressed human cell lines will soon confirm EPOR-induced gene expression and its potential tissue-specific origin.