MiRNA Expression in Patients with Gaucher Disease Treated with Enzyme Replacement Therapy

Aims: The aim of the work was to establish potential biomarkers or drug targets by analysing changes in miRNA concentration among patients with Gaucher disease (GD) compared to in healthy subjects. Methods: This study was an observational, cross-sectional analysis of 30 adult participants: 10 controls and 20 adults with GD type 1. Patients with GD type 1 were treated with enzyme replacement therapy (ERT) for at least two years. The control group was composed of healthy volunteers, unrelated to the patients, adjusted with age, sex and body mass index (BMI). The miRNA alteration between these groups was examined. After obtaining preliminary results on a group of six GD patients by the high-output method (TaqMan low-density array (TLDA)), potential miRNAs were selected for confirming the results by using the qRT-PCR method. With Diane Tools, we analysed miRNAs of which differential expression is most significant and their potential role in GD pathophysiology. We also determined the essential pathways these miRNAs are involved in. Results: 266 dysregulated miRNAs were found among 753 tested. Seventy-eight miRNAs were downregulated, and 188 were upregulated. Thirty miRNAs were significantly altered; all of them were upregulated. The analysis of pathways regulated by the selected miRNAs showed an effect on bone development, inflammation or regulation of axonal transmission in association with Parkinson’s disease. Conclusions: We revealed few miRNAs, like miR-26-5p, which are highly altered and fit the GD pathophysiological model, might be considered as novel biomarkers of disease progression but need further evaluation.


Introduction
Among lysosomal storage diseases, Gaucher disease (GD) is one of the most common. It is estimated that around 1 per 57,000 newborns is affected worldwide [1]. The disease's phenotype is associated with a mutation in the lysosomal glucosylceramidase (GBA) gene; however, it differs remarkably even in patients with the same type of mutation [2][3][4]. Three types of GD are recognized upon central nervous system (CNS) involvement: nonneuronopathic (GD type 1: OMIM #230800), acute neuronopathic (GD type 2: OMIM #230900) and chronic neuronopathic (GD type 3: OMIM #231000). The most common form of GD, covering approximately 94% of cases, is type 1, which can manifest itself at any age [5]. It is characterised by haematological changes such as thrombocytopenia, anaemia and highly pronounced splenomegaly. Sometimes, hepatomegaly may coexist, and changes occurring in the skeletal system are also important [5]. Type 2 with a very

Material and Methods
This study was an observational, cross-sectional analysis of 30 adult participants: 20 adults with GD type 1 and 10 controls. First, samples from 8 random persons (4 from the GD group and 4 from the control one) were used for the TaqMan low-density array (TLDA), a high-throughput method for identifying alterations in miRNAs expression levels. Based on these preliminary results, potential miRNAs were selected for a confirmation on a whole group (30 persons) using the more precise qRT-PCR method. Patients with GD type 1 were treated in the Metabolic Diseases and Diabetes Department in the Jagiellonian University Hospital in Krakow, part of EU reference centres-MetabERN. The study included adult people with GD, treated in our centre, who agreed to participate in the study and signed up an informed consent. GD diagnosis was established on the positive genetic or/and enzymatic tests and, in two cases, additional histopathological examination. All these patients underwent ERT for at least two years. The exclusion factors for the analysis were as follows: neoplastic diseases, PD, autoimmune diseases, pregnancy, infection, fever or lack of patient's consent to participate in the study. The control group was composed of healthy volunteers, not related to the patients, adjusted with age, sex and BMI.
All the participants underwent a detailed physical examination (age, weight, height, BMI, age and any disease history). Fasting venous blood was used for molecular biology tests. Molecular tests were performed at the Clinical Biochemistry Department of the Jagiellonian University. The clinical characteristic of groups is presented in Table 1.

RNA Isolation
The purification of total RNA from human whole blood was performed using the PAXgene Blood RNA Kit (Qiagen, Germantown, MD, USA), following the manufacturer's protocol. The RNA quality was analysed using the Tapestation 2200 instrument (Agilent Technologies, Santa Clara, CA, USA) and quantified by spectrophotometry on the NanoDrop (Thermo Fisher Scientific, Wilmington, DE, USA).

TLDA
The initial expression analysis of human microRNAs (753 unique assays) was performed by a high-throughput method using the TLDAs Panel v3.0 (Applied Biosystems, Foster City, CA, USA), as described by the manufacturer. Total RNA enriched with microR-NAs was reverse-transcribed using stem-loop primers. In order to detect low abundant microRNAs, a preamplification step was performed. The preamplified product was loaded into a TLDA, and amplification signal detection was performed in the 7900 FAST real-time thermal cycler (ABI) (Thermo Fisher Scientific, Wilmington, DE, USA). Data were exported to the DataAssist software version 3.01 (Life Technologies, now Thermo Fisher Scientific, Waltham, MA, USA) and normalised using the small-nucleolar RNA U6snRNA. Mean relative quantity (RQ) was calculated, and microRNAs differentially expressed between groups were defined as those with a fold change of >1.5 and a p-value of <0.05. The differentially expressed miRNAs based on the TLDAs results were verified using microRNA assays (TaqMan Array Human MicroRNA Panel v3.0).

qRT-PCR
Subsequently, to validate the microarray experiment, qRT-PCR was conducted on 20 patients with GD type 1 and 10 volunteers from the control group. Total RNA was reverse-transcribed using a specific stem-loop RT primer for each microRNA and the MultiScribeTM reverse transcriptase (Thermo Fisher Scientific, Wilmington, DE, USA). Then, diluted retro-transcription reaction (1:15) was mixed with the Universal PCR Master Mix, No AmpErase R UNG (2X), in the presence of individual PCR microRNAs. The PCR reactions were done in a GeneAmp R PCR System 9700 (Applied Biosystems, Foster City, CA, USA), using the following program: 95 • C for 10 min, 40 cycles of 95 • C for 15 s and 60 • C for 1 min. The relative expression of microRNAs was measured by qRT-PCR using the comparative Ct (∆∆CT) method. A two-tailed Student's t-test comparing the 2 (-∆Ct) values of the two groups was performed, and the Benjamini-Hochberg-adjusted p-value was calculated. The U6 snoRNA was used as an internal control for data normalisation.
Briefly, the DIANA-mirPath utilised predicted miRNA targets provided by the miRNAgene interactions derived from experimentally supported DIANA-TarBase v.7.0 (Department of Electrical and Computer Engineering, 382 21 Volos, Greece) or microT-CDS target prediction, if there were no hits for DIANA-TarBase. The KEGG pathways union was used to identify the main pathways targeted by predicted genes (discovered miRNA targets).

Statistical Analysis of Patients' Baseline Characteristics
For the analysis of the characteristics of the study groups, IBM SPSS Statistic Software was used. The Shapiro-Wilk test was performed to check the normality. Differences between groups were assessed using the T-test or U-test for normally and non-normally distributed continuous variables, respectively, and by the chi-square or Fischer's exact test for categorical variables. Continuous variables are presented as mean ± SD or median (interquartile range). Statistical significance was set at p < 0.05.

Results
The analysed groups: control and GD patients were comparable in gender distribution, age and anthropometric parameters. To unify patients' group as much as possible, all of the patients with GD type 1 were treated with ERT.
Using the TLDA method, we were able to perform an accurate quantitation of 753 mi-croRNAs in the study group. U6 comparison data were used for normalization. Alterations in miRNAs expression were determined in 266 cases, among which 78 were downregulated and 188 were upregulated. Eventually, among all tested miRNAs, we identified 30 miRNAs, of which the expression significantly differed between GD patients and our control data (Table 2). Interestingly, all these miRNAs were upregulated in relation to the control group.
The most prominent alteration in the GD group included miR-26b-5p, miR-31-5p, miR-29a-3p, miR-454-3p, miR-660-5p and miR-148a-3p. The miRNAs were mapped for molecular pathways to identify those, which could be specifically modified in GD patients (Table 3). Among pathways, we chose the most noticeable ones according to the p-value, the proportion of miRNAs to controlled genes and actual knowledge about the pathophysiology of the lysosomal storage disease.
The analysis using a heatmap ( Figure A1) showed a strong relationship between the TGF-beta signal transmission pathway and miR-106a-5p. In addition, a significant association was found between pathways of fatty acid synthesis and metabolism, with miR-16-5p and miR-195-5p. Interestingly, these miRNAs played also a significant role in the regulation of the Akt-signalling pathway, focal adhesion or lysine degradation pathways. Figure 1 Additionally, we noted an important impact of differentially expressed miRNA on a very complex pathway of extracellular matrix (ECM)-receptor interaction, namely miR-29a-3p or miR-let-7g-5p. complex pathway of extracellular matrix (ECM)-receptor interaction, namely miR-29a-3p or miR-let-7g-5p.

Comparison to Available Data
The biological purpose of miRNAs is related to their conservative function, common for many lifeforms, to control groups of genes by muting them when needed. They influence biochemical pathways on many levels to create a complex and tangled network of dependencies and controlling apparatus.
Nowadays, research to recognise miRNA in-depth role in the ability to regulate genes expression is blooming. As for GD, only a few miRNAs analyses have been described so far. In 2014, Siebert reported data for two miRNAs: 195-5p and 16-5p concentration changes in GD and linked them to the regulation (increase) of GCase activity [10]. Our data showed the upregulation of these two. In the presence of GCase insufficiency, the results suggested a secondary hyperconcentration of regulators and may confirm in vivo the importance of miR-195-5p and miR-16-5p in GCase secretion.

GD and Malignancies
miR-195-5p was also recently found to be downregulated in a broad spectrum of solid cancers. It was described as a cancer suppressor [14], which mediates its effect by complex and still vogue pathways. Similarly, miR-16-5p was described as a protective factor in breast carcinoma [15], and again, we showed significant upregulation of this miRNA.
Highly elevated levels of miR-195-5p and 16-5p may provide a protective role in GD in terms of solid tumours [16][17][18]. Until now, GD patients were found to be more susceptible primarily for haematologic neoplasms [19]. However, the data of breast cancer in GD are inconsistent. In 2015, European Registry did not report any incidence of breast cancer; neither did Shiran in 1993 [19,20]. In 2009, Taddei reported nine cases of breast cancer, with a risk ratio of 1.84 [21]. Therefore, these observations require further studies.
For miR-148b-3p, the connection between high concentration and poor survival in breast cancer was established [22]. Our data showed a 27.5-fold increase in this miRNA in fold change of miRNAs in the studied group

Comparison to Available Data
The biological purpose of miRNAs is related to their conservative function, common for many lifeforms, to control groups of genes by muting them when needed. They influence biochemical pathways on many levels to create a complex and tangled network of dependencies and controlling apparatus.
Nowadays, research to recognise miRNA in-depth role in the ability to regulate genes expression is blooming. As for GD, only a few miRNAs analyses have been described so far. In 2014, Siebert reported data for two miRNAs: 195-5p and 16-5p concentration changes in GD and linked them to the regulation (increase) of GCase activity [10]. Our data showed the upregulation of these two. In the presence of GCase insufficiency, the results suggested a secondary hyperconcentration of regulators and may confirm in vivo the importance of miR-195-5p and miR-16-5p in GCase secretion.

GD and Malignancies
miR-195-5p was also recently found to be downregulated in a broad spectrum of solid cancers. It was described as a cancer suppressor [14], which mediates its effect by complex and still vogue pathways. Similarly, miR-16-5p was described as a protective factor in breast carcinoma [15], and again, we showed significant upregulation of this miRNA.
Highly elevated levels of miR-195-5p and 16-5p may provide a protective role in GD in terms of solid tumours [16][17][18]. Until now, GD patients were found to be more susceptible primarily for haematologic neoplasms [19]. However, the data of breast cancer in GD are inconsistent. In 2015, European Registry did not report any incidence of breast cancer; neither did Shiran in 1993 [19,20]. In 2009, Taddei reported nine cases of breast cancer, with a risk ratio of 1.84 [21]. Therefore, these observations require further studies.
For miR-148b-3p, the connection between high concentration and poor survival in breast cancer was established [22]. Our data showed a 27.5-fold increase in this miRNA in the study group, but again, there are no available studies that investigated the mentioned issue.
Few papers described an increased number of cases with hepatocellular carcinoma (HCC) [23]. Zhao in 2017 revealed that the expression of miR-31-5p was significantly upregulated in HCC tissues [24]. This information is consistent with our analysis showing highly altered miR-31-5p that was more than 50 times higher than in the comparison group. In our tested groups, there was no HCC incidence. The question remains: is miR-31-5p just a marker for HCC or does it present the properties for the highly sensitive biomarker of the liver lesion which could reveal the destruction of hepatic cells, even when classical tests such as alanine aminotransferase (ALT) and aspartate aminotransferase (AST) are within a normal range?
All the patients in the study group were on ERT for at least two years. This treatment could influence the concentrations of many miRNAs and change the course of the disease, including complication occurrence. However, the data from European Registry [20] did not reveal differences in cancer incidence between patients ever treated with ERT and those never treated with ERT, but still, the impact of ERT on miRNAs, which are associated with malignancies, is limited as the data are missing. There are also no data on the dependence of the miRNA concentration on the ERT presence.

Other Highly Upregulated miRNA
The most prominent upregulated miRNA was found to be miR-26b-5p, which was linked with TRPS1 gene, which is responsible for regulating genes involved in the growth of bone and cartilage. Therefore TRPS1 controls chondrocyte proliferation and differentiation. MiR-26b-5p controls also bone morphogenetic protein 2 (BMP-2) that affects osteoblasts and regulates proliferation, differentiation and apoptosis [25]. In patients with GD, bone involvement is one of the main symptoms and consists of distortion of the vertebra and distension of long bones, known as Erlenmeyer flask bone deformity, with infiltration of bone marrow. It has been demonstrated that miR-26b-5p with a 66-fold increase of expression in our study group is a potential biomarker of bone involvement.

miRNA Pathways Analysis
To put altered miRNAs into the bigger picture, we used the DIANA-mirPath v3 tool for analysing miRNA-predicted targets with KEGG pathways. Nineteen pathways were significantly altered. Those, most noticeable influenced and possibly important in GD pathophysiology are "TGF-beta signalling pathway", "focal adhesion", "FoxO signalling pathway", "axon guidance" and "ECM-receptor interaction".
The TGF-beta signalling pathway is found in numerous human disorders. As it controls the growth inhibitors, the alteration of its activity is observed in cancers [26], which are more frequent in GD. In addition, TGF-beta stimulates ECM deposition that may lead to fibrosis and scarring [27]. The higher prevalence of liver fibrosis in GD was also shown [28].
The focal adhesion pathway focuses on the regulation of call-matrix interaction. The central point of that pathway-focal adhesion kinase (FAK) was roughly studied in terms of cancer progression [29][30][31]. Likewise, FAK was found to promote proinflammatory gene expression via TNF-α or IL-1β [32].
The FoxO signalling pathway, on the other hand, is involved in oxidation stress and may play an important role in neurodegenerative disease like Parkinson's [33,34].
These data confirmed the previously described role of inflammatory processes and the disrupted pathways of cell migration and interaction with the ECM, that consists of a mixture of structural and macromolecules and serves an important role not only in tissue morphogenesis but also in the regulation of cell and tissue structure and function.

miRNA and PD
Since the relationship between PD and GD was revealed, many studies were trying to establish the underlying pathophysiology. It was shown that patients with PD had a decreased level of GCase in the brain tissue. Moreover, even carries of the single GBA gene mutation present tendency towards PD [36]. This finding showed that GD patients should be more susceptible to the occurrence of PD, which was confirmed in [36]. Research on the expression of genes related to neuroinflammation and synaptic plasticity are ongoing [12]. MiR-26b and miR-106a-5p acting through HSPA8 [37], as well as miR-7 and miRNA-153, are found to increase the concentration of alpha-synuclein [38], the primary accumulation material in PD. We noted a 14-fold increase of miR-106-5p in our study group, which is consistent with the known data.