Expression profiling of muscle invasive and non-invasive bladder tumors for biomarkers identification related to drug resistance, sensitivity and tumor progression

Abstract

Many bladder cancer (BC) patients with early disease are asymptomatic and diagnosed at advanced stage when the therapeutic options are limited. This necessitates the development of reliable predictive molecular biomarkers that will ensure a positive therapeutic response in every patient. The aim of this study was to screen for alterations in gene expression levels related to drug sensitivity and resistance that may be further explored as potential predictive therapeutic biomarkers. Gene expression analysis of the 168 genes from two panels for Cancer drug resistance and metabolism (PAHS004) and Cancer Drug Targets (PAHS507z) was performed. A total of 47 transitorial cell bladder cancer samples of stage pTa, pT1, pT2 were investigated using the pooling method, which allows reducing the effect of biological variation and detecting only significant expression changes. Differential gene expression was calculated using the ΔΔCt method with GPDH as a housekeeping gene. The 4.0-fold change in gene expression was used as the cut-off threshold to determine upregulation or downregulation compared to normal bladder tissue (negative control). Significance of the differences in the expression profiles was assessed by nonparametric one-way analysis of variance (ANOVA) with Dunn’s multiple comparison tests and Mann-Witney test. We demonstrated a correlation of tumor invasion and several up-regulated genes related to chemotherapy resistance. For the first time, this study demonstrated overexpression of CDK8, CDK9, FIGF, HDAC11, IGF1 and PDGFRA genes in muscle-invasive bladder carcinomas. These genes and their proteins could be used as potential biomarkers for bladder cancer progression or prospective therapeutic targets.

Keywords:

• Bladder
• cancer
• biomarker
• muscle-invasive
• drug
• resistance

Introduction

The leading role of the environmental risk factors and occupational exposures in the initiation and development of bladder cancer (BC) is well known. BC is associated with exposure to workplace carcinogens and is linked to lifestyle, environmental factors, aging, male gender and Caucasian race [1]. The latter may be partly due to certain genetic polymorphisms like GSTM1 0/0 and NAT2 slow acetylators, which are more prevalent in Caucasian than non-Caucasian populations [2]. Even if tobacco smoking could be reduced through appropriate actions like healthy lifestyle promotion programs, and more rigorous control of the work-safety guidelines implementation in industries were to be introduced, it is still almost impossible to mitigate the role of occupational risk factors, ethnicity, gender and aging, and as a result the disease incidence levels remain steady [1].

Many bladder cancer patients with early disease are asymptomatic and are diagnosed at an advanced stage when the therapeutic options are limited. The initial treatment for bladder cancer consists of transurethral resection (TURB) followed by adjuvant chemotherapy. Although low grade pTa and pT1 non-muscle invasive bladder cancers (NMIBC) could be eradicated completely by TURB, it commonly reappears and progresses to muscle-invasive bladder cancer (MIBC) [3]. The recurrence phenomenon of non-invasive bladder cancer makes it one of the costliest among the common neoplasms to handle per patient from diagnosis to death [4]. Moreover, the timely diagnosis of NMIBC recurrences and MIBCs is crucial, because a delay in cancer treatment can be life-threatening. The existing treatment protocol for MIBC with bladder resection provides only a 50% five-year survival [1]. To improve these rates, prior to surgery or radiotherapy, a course of cisplatin-based neoadjuvant chemotherapy may be recommended [5]. Unfortunately, not all the patients respond to neoadjuvant chemotherapy. For responders, treatment has a major positive impact on overall survival, whereas for non-responders only side effects are reported [6]. This requires the development of reliable predictive molecular biomarkers to personalize the therapy.

The breakthroughs in genetics over the last 10 years brought the discovery of a number of associated biomarkers such as angiogenesis related factors, altered p53 expression, vascular endothelial growth factor, fibroblast growth factor, fibroblast growth factor receptor-3, thrombospondin-1, multidrug resistance gene expression etc. [1]. However, due to insufficient predictive power, none of them has been incorporated in the standard management of protocols. Molecular tumor profiling has potential value for tumor classification recurrence risk assessment, as well as tailoring the treatment to the specifics of every patient.

The aim of this study is to screen for alterations in gene expression levels related to drug sensitivity and resistance that may be further researched as potential predictive therapeutic biomarkers. For this screening, we applied the pooling method, because it allows reflection of general trends of the studied group at the expense of individual biological variation, which leads to the detection only of genes with the highest changes in the expression.

Subjects and methods

Ethics statement

The study protocol was approved by the Ethics Committee of the Medical University of Sofia, Bulgaria (protocol No 76/23.04.2012). Each patient provided written informed consent and the relevant clinical data.

Patients and samples

Bladder tissue samples were collected from 52 Caucasian patients who underwent surgical treatment for bladder cancer at the Clinic of Urology of Aleksandrovska University Hospital, Sofia, Bulgaria. Data about lifetime use of tobacco products and occupational history were collected. All the samples were immediately frozen in liquid nitrogen and stored at −80 °C. All samples underwent histological assessment. Normal bladder tissue was taken from patients suffering from benign prostatic hyperplasia. Complete staging procedures were carried out to determine precisely the primary tumor size - T, nodal involvement - N and presence of distant metastases - M, according to the TNM International Staging system for bladder cancer. None of the patients had received any prior radiotherapy or chemotherapy.

RNA isolation, pool preparation and reverse transcription

Total RNA was isolated using the ready-to-use TRIzol® LS Reagent (Invitrogen, Cat. No 10296-010) according to the manufacturer’s protocol. The double-stranded cDNA was prepared from 1 µg of extracted total RNA using the RT2 First Strand Kit for cDNA synthesis (Qiagen, Cat. No 330401), which contains a step for genomic DNA elimination. RNA yield and purity were determined using NanoDrop® ND-1000, and RNA integrity, by 1% gel electrophoresis.

Five different pools were prepared from the isolated RNA. Each pool was composed of tumor samples from the same tumor stage: pTa; pT1; pT2; chronic inflammation (precancerous group); or normal bladder tissue (negative control). RNA samples comprising each pool were diluted to 100 ng/μL. One microgram of RNA from each sample was added to the respective pool.

Quantitative real time PCR (qRT-PCR) analysis

The relative mRNA expression of genes involved in the regulation of cell cycle, DNA repair, drug resistance, drug metabolism, growth factor receptors and transcription factors were analyzed with Cancer drug resistance & metabolism (PAHS-004) and Cancer drug targets PAHS-507Z) RT² Profiler PCR Arrays from (SA Biosciences, Qiagen) according to the manufacturer’s instructions. Quantitative PCR was performed with the usage of RT² qPCR SYBR Green MasterMix in a 96-well plate format using the ABI 7500 instrument (Applied Biosystems). The instrument program was as follows: initial denaturation at 95 °C for 10 min and 40 cycles of 95 °C for 15 s followed by 60 °C for 1 min.

The gene expression assays on PAHS004 and PAHS507z PCR array were performed in five pools from pTa, pT1 and pT2 bladder cancer samples, precancerous group and negative control (normal bladder tissue). To confirm pool results, three randomly chosen individual tumor samples from each tumor stage (pTa, pT1 and pT2) were tested for PAHS004 pathway.

Statistical analysis

Differential gene expression in the experimental samples relative to the control samples was calculated using the ΔΔCt method with GPDH as a housekeeping gene. According to the manufacturer’s protocol, a 4.0-fold change in gene expression was used as the cut-off threshold to determine upregulation or downregulation compared to the negative control. Data analysis was performed by web-based SuperArray software - RT² Profiler PCR Arrays Data Analysis version 3.5. Only results that meet the aforementioned quality criteria were used for further data analysis.

Statistical significance testing was assessed by nonparametric one-way analysis of variance (ANOVA) with Dunn’s multiple comparison tests and Mann-Witney test. Nonpaired t-test was performed to analyse the difference in regard to muscle invasion in the tested cancer samples. p < 0.05 was considered to indicate a statistically significant difference; ‘*’ denotes p < 0.05, ‘**’ - p < 0.01 and ‘***’ - p < 0.001.

Results and discussion

Characteristics of bladder cancer samples

Of the 52 patients who underwent transurethral resection, 41 had histologically confirmed primary transitional cell carcinoma; four patients had prostatic metastasis in the bladder; five ones showed chronic bladder inflammation, one, nonmalignant bladder lesion and one case was not investigated because of improper handling and damaged sample. Demographic, clinical and pathological characteristics are presented in Table 1.

Table 1. Demographic, clinical and pathological characteristics of bladder cancer patients and samples included in the study.

	Female (n = 7, 17.1%)	Male (n = 34, 82.9%)	Total (n = 41, 100%)
Patients’ characteristics
Mean age at diagnosis	65.4 (range 52–80)	66 (range 42–80)	65.9 (range 42–80)
Smoking status (10–20 cigarettes per day)	0	10 (29.4%)	10 (24.3%)
Professional risk factors	0	1 (2.9%)	1 (2.4%)
Mytomicin-C treatment	3 (42.85%)	4 (11.76%)	7 (17.07%)
Tumor characteristic
Primary tumors	2 (28.5%)	14 (41.17%)	16 (39.02%)
Recurrent tumors	5 (71.4%)	20 (58.8%)	25 (60.97%)
Tumor stage:
pTa	1 (14.28%)	9 (26.47%)	10 (24.4%)
pT1	2 (28.57%)	11 (32.35%)	13 (31.7%)
pT2	4 (57.14%)	14 (41.18%)	18 (43.9%)
Tumor grade:
G1	3 (42.85%)	17 (50%)	20 (48.8%)
G2	3 (42.85%)	7 (20.6%)	10 (24.4%)
G3	1 (14.3%)	10 (29.4%)	11 (26.8%)

Changes in gene expression

The results of the gene expression assay on PAHS004 and PAHS507z PCR array show statistically significant change (p < 0.05) for the cell cycle, growth factors and receptors, PI-3 kinases and phosphatases, drug resistance and metabolism, DNA repair, hormone receptors, and transcription factors functional groups mainly between muscle-invasive pT2 and precancerous, and non-invasive bladder tumors (pTa and pT1) (Figure 1). Statistically significant difference (p < 0.01) was also found between the precancerous pool and pTa for the genes related to drug resistance and metabolism. Although there was no significant difference in the expression levels of the genes coding for Histone Deacetylases between precancerous, pT1 and pT2 bladder tumors, the level of gene expression was 5 to 8 times higher than those in normal bladder tissue (Figure 1).

Figure 1. Relative expression level of genes from PAHS-507 Z and PAHS-004 PCR arrays of precancerous, non-invasive (pTa, pT1), and invasive bladder pT2 tumors.

Note: 4.0-fold change in gene expression was used as cut-off threshold to determine upregulation or downregulation compared to negative control; ‘*’ denotes p < 0.05, ‘**’ - p < 0.01, and ‘***’ – p < 0.001 significant difference in gene expression between pT2 and other bladder cancer samples; ‘^^’ denotes p < 0.01 significant difference in gene expression between pTa and precancerous samples.

When analysing each pool separately, in the precancerous pool upregulation was found for genes: ABCC1, TOP2A, TOP2B, TXN and ATF2 (drug resistance and metabolism); EGFR, MDM2, PARP1 and FLT1 (growth factor and hormone receptors); AKT1, AKT2 (receptor tyrosine kinase signalling); CTSB, CTSD, CTSS (cathepsins); HDAC1, HDAC7 (histone deacetylases); AURKA and PIK3C3 (PI-3 kinases and phosphatases); as well as genes coding for the cell cycle kinases CDK1, CDK4, CDK5, CDK7 and MDM2 (Figure 2). Among these genes, AKT2, CTSS, MDM2 and CDK7 had more than 20 times higher expression compared to the negative control.

Figure 2. Relative expression levels of genes related to Cancer drug resistance and metabolism; Growth factors and Hormonal receptors; Transcription factors, Cathepsins and Histone deacetylases; PI-3 kinases and phosphatases and Cell cycle proteins of precancerous pools, non-invasive (pTA, pT1) and invasive bladder pT2 tumors.

Note: 4.0-fold change in gene expression was used as cut-off threshold to determine upregulation or downregulation compared to negative control.

In the pTa urinary bladder tumors, the majority of the tested genes had expression level less than the ±4.0-fold change cut-off threshold. However, ABCC1, ABCC3, CYP3A5 and AHR (drug resistance and metabolism), EGFR, ERBB2, HRAS and PPARG (growth factor and hormone receptors); AKT1 and AKT2 (receptor tyrosine kinase signalling); CTSD and CTSS (cathepsins) were up-regulated. No genes related to the cell cycle regulation were found to be overexpressed (Figure 2).

Among the pT1 bladder cancer samples, up-regulated genes from the drug resistance and metabolism group were ABCC1, ABCC3, CYP3A5, TOP2B, TXN, AHR and ATF2; growth factor and hormone receptors - EGFR, ERBB2, HRAS, MDM2, PPARG and RXRA; transcription factors, cathepsins and histone deacetylases - AKT1, AKT2, CTSB, CTSD, CTSS, HDAC1, HDAC7 and HDAC8; PI-3 kinases & phosphatases - AURKA, AURKC, PIK3CA, PIK3C3 and PRKCA; as well as genes related to the cell cycle progression - CDK1, CDK7 and MDM2 (Figure 2).

In the muscle-invasive urinary bladder tumors from pT2 stage, there were up-regulated genes in the cancer drug resistance and metabolism functional group (ABCC1, ABCC3, CYP3A5, TOP2A, TOP2B, TXN, AHR and ATF2); genes coding for growth factor and hormone receptors (EGFR, ERBB2, HRAS, MDM2, PARP1, PPARG, RXRA, RXRB and FLT1); genes coding for transcription factors, cathepsins and histone deacetylases (AKT1, AKT2, CTSB, CTSD, CTSS, HDAC1, HDAC7 and HDAC8); genes coding for PI-3 kinases & phosphatases (AURKA, AURKC, PIK3CA, PIK3C2A, PIK3C3, PRKCD and MTOR); genes for the cell cycle proteins (CDK1, CDK2, CDK4, CDK5, CDK7, CDK8, CDK9, CDKN2D, MDM2 and MDM4). The highest gene expression level was observed for CYP3A5, MDM2, PPARG, CDK1, ESR2 and FLT4. The upregulation of the other genes was within the range between 4- and 25-fold (Figure 2).

Gene expression of muscle-invasive versus non-invasive bladder tumors

When comparing the pT2 muscle-invasive tumors with pT1 non-invasive ones, we observed genes that are up-regulated in pT2 muscle-invasive tumors only. These were ARNT, BCL2, CYP1A1, EPHX1, GSTP1, MVP, PTGS2, RB1, SOD1, CDK2, CDK4, CDK5, CDK8, CDK9, MDM4, FIGF, FLT1, FLT4, IGF1, IGF1R, KIT, MET, PDGFRA, PDGFRB, ERBB4, KRAS, NRAS, HDAC3, HDAC6 and HDAC11 (Figure 3). The gene expression level was within the 5–40 range in comparison to the negative control. These results showed statistically significant difference between muscle invasive (pT2) and non-invasive (pT1) bladder tumors (p < 0.0001, Mann-Witney test).

Figure 3. Relative expression level of genes of PAHS-507 Z and PAHS-004 PCR chips of non-invasive (pT1) and muscle-invasive (pT2) bladder tumors.

Note: 4.0-fold change in gene expression was used as cut-off threshold to determine upregulation or downregulation compared to negative control.

The applied confirmation assays show that, in general, the results of the expression analysis in the pools match those of the individual samples. Both in the individual sample and in the pool analysis, the expression of genes involved in functional pathways related to drug sensitivity and metabolism among pTa tumors was higher than that among pT1 tumors, although some genes from these pathways showed different expression levels (data not shown). In the pool analysis, genes with reduced expression that were not detected in the individual tumors were identified. This finding is expected due to the difference in the number of samples in the two types of assays: between 8 and 12 tumors for the pools and 3 tumors for the individual studies, respectively.

Pool approach and qRT-PCR in cancer studies

In order to assess the bladder cancer samples for alterations in drug resistance/sensitivity pathways and possible drug targets it is important to use a safe and reliable method that is widely applied in clinical practice. That is why we applied qRT-PCR, which is considered to be both an accurate and safe method routinely used in medicine on a daily basis. To identify common differentially expressed genes, we applied the pooling method, which allows reducing the effect of biological variation and detecting only significant expression changes [7]. The same pooling approach has been used previously in other gene expression research for nasopharyngeal carcinomas [8], thyroid adenomas [9], aldosterone-producing adenomas [10], pancreatic tumors [11] as well as sporadic colorectal carcinomas [12].

Possible role of inflammation in bladder cancer

It is noteworthy that the TXN, EGFR, MDM2, PARP1, FLT1, PIK3C3, CDK4, CDK5 and CDK7 genes were significantly expressed among the precancerous pool and pT2, but not to the same extent in the pTa and pT1 bladder tumors. We assume that these genes are associated with the inflammatory processes, because on the one hand, the precancerous pool is composed of nonmalignant lesions of the bladder - chronic inflammation (cystitis), and on the other hand, inflammation can be seen in advanced tumors as a secondary process. Some scientific reports link bladder cancer initiation to chronic inflammation and dysregulated expression and activity of certain chemokines [13–15]. This chronic inflammatory milieu associated with macrophage and T cell accumulation, chronic activation of macrophages, abnormal angiogenesis and DNA damage due to the presence of reactive oxygen species has been implicated in the potential progression from an inflammatory environment to cancer initiation and is a target of novel therapeutic approaches in bladder cancer treatment [16–18].

Early stage cancer and abnormalities in drug metabolism pathways

In the pTa urinary bladder tumors, the majority of the tested genes had expression levels close to normal bladder epithelia and no genes related to the cell cycle regulation were found to be overexpressed. However, there was overexpression of genes responsible for drug sensitivity, resistance and metabolism in both pTa and pT1 non-invasive bladder cancer - ABCC1, ABCC3, CYP3A5, AHR, EGFR and ERBB2 (Figure 2).

The proteins encoded by the ABC genes (ABCC1 and ABCC3) are members of the superfamily of ATP-binding cassettes, or ABC transporters. These proteins transport a variety of xenobiotics and physiological substrates through intra- and extra-cellular membranes, thus contributing to chemotherapeutic resistance to gemcitabine in bladder cancer patients [19]. Similar to lung carcinoma, pancreatic carcinoma and glioblastoma multiforme – we demonstrate that the expression of ABCC3 is also elevated in bladder cancer and is associated with tumor grade. The highest expression is found in pT2 bladder tumors.

Concomitant up-regulation of ABCC3, EGFR and ERBB2 is in accordance with reports on HER2 positive mammary tumors, where increased expression of ABCC3 was observed as a result of ABCC3 gene amplification [20, 21]. Indeed, a number of studies suggest a possible novel therapeutic approach by inhibiting these molecules [22].

CYP3A5 and AHR are involved in the metabolism of a number of exogenous and endogenous substances, anti-tumor drugs and steroid hormones. AHR (Aryl Hydrocarbon Receptors) regulates the function of xenobiotic-metabolizing enzymes of the cytochrome P450 superfamily. Dysregulation of the AHR gene leads to tumor initiation, promotion and progression [23]. AHR directly binds polycyclic aromatic hydrocarbons (PAHs) and dioxin and activates xenobiotic metabolism, histone modification and tumorigenesis [24, 25]. Expression of the CYP3A5 gene varies greatly across populations and in the urinary bladder it contributes to metabolitic Phase I activation of bladder carcinogens [26]. Upregulation of this gene has not yet been reported in bladder tumors.

Muscle invasive disease and deregulation of cell cycle control, neovascularization and chromatin remodelling systems

Upregulation of actin, cathepsins and histone deacetylases (Figure 2) is associated with tumor development. The AKT1 and AKT2 genes are putative oncogenes that suppress apoptosis and are researched as potential drug targets [27]. Cathepsins are proteases that play a leading role in the progression and invasion of malignancies. Cysteine cathepsins (CTSB, CTSS) break down the extracellular matrix, thereby facilitating growth and invasion into surrounding tissues, as well as vascularization. Increased CTSD concentrations in urothelial carcinomas have been already demonstrated, and increase of CTSD concentrations in serum and urine can be used as a tumor marker in urothelial carcinoma [28, 29].

Genes coding for histone acetylases are expressed in bladder tumors, colorectal and prostate carcinomas [30]. Increased expression of HDAC1 and HDAC2 is observed in non-invasive bladder papillary tumors and is associated with worse prognosis in higher-grade tumors [31].

In this study, we compared the gene expression levels between the muscle invasive and non-invasive bladder cancers (Figure 3). We demonstrated the upregulation of numerous genes in pT2 samples in comparison to pT1 tumors: ARNT, BCL2, CDK2, CDK4, CDK5, CDK8, CDK9, CYP1A1, EPHX1, FIGF, FLT1, FLT4, GSTP1, HDAC3, HDAC6, HDAC 11, IGF1, IGF1R, KIT, KRAS, NRAS, MDM4, MET, MVP, PDGFRA, PDGFRB, PTGS2, RB1, SOD1 and ERBB4. Among them BCL2, CDK2, CDK4, CDK5, CDK8, CYP1A1, FLT1, FLT4, GSTP1, HDAC3, HDAC6, IGF1R, KIT, KRAS, MDM4, MET, MVP, PDGFRB, RB1, SOD1 and ERBB4 encode various receptors or proteins involved in drug resistance, cell cycle activation, vascularization and invasion, and are already reported to be involved in bladder cancer progression or to influence therapeutic response [32–39].

In this study, we demonstrated for the first time, the overexpression of CDK8, CDK9, FIGF, HDAC11, IGF1 and PDGFRA genes in muscle-invasive bladder carcinomas in comparison to non-invasive bladder carcinomas. The CDK8 and CDK9 genes code for cyclin-dependent kinases. The expression levels of CDK8 are considered as potential biomarkers in breast cancer [40]. Dysregulation of the antiapoptotic Cdk9 signalling pathway leads to malignant transformation of lymphocytes and monocytes [41, 42]. The involvement of CDK9 in the genesis of advanced prostate cancer has been previously demonstrated by Mohaparta et al. [43].

FIGF (c-fos induced growth factor, VEGFD) promotes angiogenesis, lymphangiogenesis and endothelial cell growth. Increased FIGF expression is associated with lung cancer [44]. It has been reported that the expression level of VEGFD (FIGF) is significantly reduced after andrographolide treatment, which inhibits tumor growth of hepatocellular carcinoma in vivo [45].

The HDAC11 gene encodes a class IV histone deacetylase. HDAC11, together with HDAC6, is involved in the construction of a protein complex [46]. In our study, we found that the expression of HDAC6 and HDAC11 in muscle-invasive carcinomas was much higher than that in non-invasive carcinomas (Figure 3). Increased expression of HDAC11 has also been reported in other malignancies, which is in line with our results. Moreover, HDAC11 depletion is sufficient to cause cell death and inhibition of metabolic activity in colon cancer (HCT-116), prostate cancer (PC-3) and mammary gland (MCF-7) cell lines. At the same time, according to Deubzer et al., a decrease in the level of HDAC11 in healthy tissue does not alter the metabolic activity and cell viability, making HDAC11 suitable drug target for various oncological diseases [47]. It has been suggested that inhibitors that block histone deacetylases can be used in antitumor therapy. HDAC inhibitors can modulate the expression of genes related to cell cycle, apoptosis and angiogenesis. These inhibitors, in combination with DNA-demethylating agents and chemoprophylaxis are potential novel therapeutic approaches to bladder cancer treatment.

IGF1 gene encodes a protein similar in function to insulin and involved in mediating growth and development. There is data that upregulation of IGF1 by tumor-associated macrophages promotes the proliferation and migration of epithelial ovarian cancer cells, but the potential mechanism of action in bladder cancer is still unclear [48].

PDGFRA is a cell surface tyrosine kinase receptor for members of the platelet-derived growth factor family. These growth factors are mitogens for cells of mesenchymal origin. PDGFRA overexpression has been reported in patients with oral cancer with significant trend of overpresentation among tobacco users [49].

In conclusion, we have applied robust and widely available method, to screen bladder cancer samples for alterations in drug resistance/sensitivity pathways and for possible biomarkers related to cancer progression. We have found a correlation of tumor invasion and several up-regulated genes involved in resistance to common chemotherapeutics, which after careful consideration and replication of the results in larger sample sets may provide a rationale for therapeutic choice on the basis of the molecular characteristics of the patient’s tumor.

For the first time, it has been demonstrated overexpression of CDK8, CDK9, FIGF, HDAC11, IGF1 and PDGFRA genes in muscle-invasive bladder carcinomas. These genes and coding proteins could be used as potential biomarkers for bladder cancer progression or therapeutic targets in the future. This necessitates examination of their expression level, methylation and mutation status among larger patient’s cohort.

Disclosure statement

No potential conflict of interest was reported by the authors.

Availability of data and materials

The data that support the findings of this study are available from the corresponding author, upon reasonable request.

Additional information

Funding

This work was supported by Bulgarian Ministry of Education and Science, NSF under Grants: DM13/4, 2017 and KP-06-OPR01/3-2018, the National Program “Young Scientists and Postdoctoral Students” and Medical University of Sofia, Grant 13-D, 2012.