Differential microRNA Expression in USP8-Mutated and Wild-Type Corticotroph Pituitary Tumors Reflect the Difference in Protein Ubiquitination Processes

Background: USP8 mutations are the most common driver changes in corticotroph pituitary tumors. They have direct effect on cells’ proteome through disturbance of ubiquitination process and also influence gene expression. The aim of this study was to compare microRNA profiles in USP8-mutated and wild-type tumors and determine the probable role of differential microRNA expression by integrative microRNA and mRNA analysis. Methods: Patients with Cushing’s disease (n = 28) and silent corticotroph tumors (n = 20) were included. USP8 mutations were identified with Sanger sequencing. MicroRNA and gene expression was determined with next-generation sequencing. Results: USP8-mutated patients with Cushing’s disease showed higher rate of clinical remission and trend towards lower tumor volume than wild-type patients. Comparison of microRNA profiles of USP8-mutated and wild-type tumors revealed 68 differentially expressed microRNAs. Their target genes were determined by in silico prediction and microRNA/mRNA correlation analysis. GeneSet Enrichment analysis of putative targets showed that the most significantly overrepresented genes are involved in protein ubiquitination-related processes. Only few microRNAs influence the expression of genes differentially expressed between USP8-mutated and wild-type tumors. Conclusions: Differences in microRNA expression in corticotropinomas stratified according to USP8 status reflect disturbed ubiquitination processes, but do not correspond to differences in gene expression between these tumors.


Introduction
Pituitary neuroendocrine tumors (PitNETs) represent about 10-20% of all intracranial neoplasms in adults. They may arise from various kinds of secretory cells of pituitary gland, including corticotroph cells, which produce adrenocorticotropic hormone (ACTH). Corticotroph PitNETs commonly cause ACTH-dependent Cushing's disease (CD); however, a notable proportion of tumors originated from corticotropic pituitary cells are endocrinologically non-functioning and classified as silent corticotroph tumors commonly referred to as silent corticotroph adenomas (SCA). Both active and silent corticotroph PitNETs share a similar molecular profile [1,2].
Recently, notable progress in the understanding of pathogenesis of CD has been made [3], including the discovery of recurrent USP8 mutations [4][5][6]. These mutations are observed in approximately 30-40% of patients suffering from Cushing's disease as well as in silent tumors [4,[6][7][8][9]. Patients with Cushing's disease with and without USP8 mutations have a slightly different clinical profile according to previously published data [4,6,8,[10][11][12][13][14][15]. The studies showed that USP8 mutation is related to lower tumor size [4,5,8] and clinical remission after surgery [8] [12]. Additionally, differences in expression of possible molecular predictive markers as MGMT or somatostatin receptor were also observed in this group of patients [2,8]. Perhaps testing the USP8 mutation status in patients could provide some kind of clinically useful information; however, the clinical results published up to today are insufficient, and molecular consequences of this mutation are only partially recognized.
USP8 gene encodes for deubiquitinase enzyme involved in the regulation of proteasomal protein degradation. USP8 mutations are small, single codon deletions or missense variants that occur in the region involved in binding 14-3-3 proteins family members. Thus, these changes impair interactions between USP8 and 14-3-3 proteins, which normally suppress deubiquitinase activity [16]. As a result of mutation, USP8 activity is enhanced and leads to preventing proteasomal degradation of particular proteins and dysregulation of natural protein turnover. This was clearly shown in in vitro experiments, providing the explanation of the sustained EGFR signaling in USP8-mutated corticotroph PitNETs [4,5]. Since USP8 deubiquitinase has many molecular substrates, its impaired functioning has potentially a wide effect on the protein level. The pleiotropic effect of the mutation is reflected by differences in gene expression between USP8-mutated (USP8mut) and USP8-wild type (USP8wt) tumors [1,2]. Accordingly, differences in expression of particular proteins related to corticotroph tumors growth were also found [15]. The understanding of the biological difference between wild type and mutated tumors appears important, since USP8 mutation may potentially serve as prognostic and predicting factor [2].
The aim of this study was to compare the profiles of microRNA (miRNA) expression in corticotroph tumors stratified according to USP8 mutational status and to determine the potential role of differential miRNA expression. Moreover, mutations of another deubiquitinase-encoding gene (USP48) contribute to pathogenesis of some USP8wt tumors [17,18]. Mutations of both genes were determined and taken into account in data analysis.

Patients and Samples
Pituitary tumor samples were collected during transsphenoidal surgery and fixed in formalin for routine diagnostic procedures, including immunohistochemical and ultrastructural evaluation. Archival formalin-fixed paraffin-embedded (FFPE) tissue samples from 48 patients, including 28 samples from patients with Cushing's disease and 20 samples from patients with SCA from years 2013-2017, were included. Patients were diagnosed according to WHO criteria applied during the time of tissue sampling [19]. Diagnosis was based on results of immunohistochemical staining for pituitary hormones and Ki-67 labeling as well as commonly accepted ultrastructural features of corticotroph tumors [20]. For this study, all the tumor samples were reevaluated histopathologically by one pathologist to confirm the diagnosis and determine tumor tissue content within each sample.
The diagnosis of Cushing's syndrome/hypercortisolism was based on standard hormonal criteria: increased urinary free cortisol (UFC) in three 24 h urine collections, disturbances of cortisol circadian rhythm, increased serum cortisol levels accompanied by increased or not suppressed plasma ACTH levels at 8 a.m., and a lack of suppression of serum cortisol levels to < 1.8 µg/dL during an overnight dexamethasone suppression test (1 mg at midnight). The pituitary etiology of Cushing's disease was confirmed based on the serum cortisol levels or UFC suppression < 50% with a high-dose dexamethasone suppression test (2 mg q.i.d. (lat. quater in die = four time a day) for 48 h) or a positive result of a corticotrophin-releasing hormone stimulation test (100 mg i.v. (intravenously)) and positive pituitary magnetic resonance imaging. In the group of SCAs, none of the patients had any evidence of hypercortisolism based on clinical signs and symptoms as well as basic laboratory tests. ACTH levels were assessed using IRMA (immunoradiometric assay) (ELSA-ACTH, CIS Bio International, Gif-sur-Yvette Cedex, France). The analytical sensitivity was 2 pg/mL (reference range: 10-60 pg/mL). Serum cortisol concentrations were determined by the Elecsys 2010 electrochemiluminescence immunoassay (Roche Diagnostics, Mannheim, Germany). Analytical sensitivity of the assay was 0.02 µg/dL (reference range: 6.2-19.4 µg/dL). UFC was determined after extraction (liquid/liquid with dichloromethane) by electrochemiluminescence immunoassay (Elecsys 2010, Roche Diagnostics, Mannheim, Germany)-reference range: 4.3-176 µg/24 h. Bilateral inferior petrosal sinus sampling was used as a routine investigation tool in any patient with proven ACTH-dependent Cushing's syndrome and negative or equivocal MRI findings (intrasellar lesion ≤ 6 mm) [21]). A macroadenoma was defined as a tumor with at least one diameter exceeding 10 mm, and the tumor volume was assessed with the diChiro Nelson formula (height × length × width × π/6). Invasive growth of the tumors was evaluated using Knosp grading [22]. Tumors with Knosp grades 0, 1 and 2 were considered non-invasive, while those with Knosp 3 and 4 were considered invasive.
Overall characteristics of the patients are presented in Table 1 while details are provided in Table S1. The content of tumor tissue in each FFPE sample ranged between 80 and 100% (median 99%) (details in Table S1). The study was approved by the local Ethics Committee of Maria Sklodowska-Curie Institute-Oncology Center in Warsaw, Poland. Each patient provided informed consent for the use of tissue samples for scientific purposes.
DNA and total RNA from FFPE samples was isolated using RecoverAll™ Total Nucleic Acid Isolation Kit for FFPE (Thermo Fisher Scientific, Waltham, MA, USA) and measured using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). The samples were stored at −70 • C.

Genomic Mutation Testing
The presence of point mutation at the USP8 hotspot (exon 14) and USP48 hotspot (exon 10) was determined using Sanger sequencing. DNA was PCR amplified with Fast-Start Taq DNA Polymerase (Roche Diagnostics, Mannheim, Germany) using the GeneAmp 9700 PCR system (Applied Biosystems, Foster City, CA, USA). The PCR product was purified using ExoSAP-IT (USA Affymetrix, Cleveland, OH, USA), labeled with BigDye Terminator v.3.1 (Applied Biosytems, Foster City, CA, USA) according to the manufacturer's instructions and analyzed by capillary electrophoresis with the ABI PRISM 3300 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA), as described previously [1]. The following sequences of PCR primers were used: 5 -TCCACCCCTCCAACTCATAA and 5 -CTGACAGATTCAGAGTAGAAACT for USP8 mutation testing as well as 5 -GCCCGGCTAAAGAATAAACA and 5 -TGCCTGCTATAATCCTGGAAA for identification of USP48 variants.

Determining miRNA Expression Profile with Next Generation Sequencing (NGS)
The quality of small RNA fractions was assessed using Agilent 2100 Bioanalyzer with Small RNA Kit chip (Agilent, Santa Clara, CA, USA) and measured with Qubit RNA HS Assay Kits (Thermo Fisher Scientific, Waltham, MA, USA). One µg of total RNA was used to sequencing library construction with an Ion Total RNA-Seq Kit v2 (Thermo Fisher Scientific), according to the manufacturer's protocol. Ion Xpress™ RNA-Seq Barcode Kit, which allows for multiplexed sequencing, was used for hybridization and ligation of RNA adapters. RNA reverse transcription and subsequent cDNA purification and library size selection were performed using Nucleic Acid Binding Beads. cDNA was PCR-amplified, followed by DNA purification and size selection. Amount and size distribution of the amplified DNA was determined using Bioanalyzer 2100 using a High Sensitivity DNA Kit (Agilent, Santa Clara, CA, USA). The length of miRNA ligation products in barcoded libraries ranged between 94 and 114 bp. Template preparation for clonal amplification of up to four miRNA libraries at a concentration of 18 pM, and loading of the PI chip, were performed using Ion Chef instrument, with Ion PI™ Hi-Q™ Chef Kit (Thermo Fisher Scientific, Waltham, MA, USA). An Ion Proton sequencer (Thermo Fisher Scientific, Waltham, MA, USA) was used for sequencing. Unmapped bam files were converted into fastq files with a bamToFastq script from bedtools. Read mapping to known human miRNAs (according to miRBase release 22) and reads quantification were performed using miRDeep2.14. Data normalization and differential expression analysis were performed using DESeq2. Filtration for low-expression miRNAs and miRNAs genes with less than five sequencing reads in at least half of the samples were excluded. Fold change of expression (FC) calculated as ratio of the normalized read-count value in USP8-mutated and USP8-wt tumors was used as a measure of expression difference. Differentially expressed miRNAs were defined as those with |FC| > 2 and adjusted p-value < 0.05.

Gene Expression Profiling
Gene expression profiles were determined in 24 FFPE samples of corticotroph tumors by sequencing of amplicon-based library representing whole transcriptome, as described in detail previously [1]. Ion AmpliSeq™ Transcriptome Human Gene Expression Kit (Thermo Fisher Scientific, Waltham, MA, USA) was used for library preparation and semiconductive sequencing technology with Ion Proton instrument, PI chip, and the sequencing reagents included in Ion PI™ Hi-Q™ Sequencing 200 Kit according to the manufacturer's instructions (Thermo Fisher Scientific, Waltham, MA, USA). Data processing was performed using Bioconductor packages in R environment as described [1]. Differentially expressed genes (DEGs) were defined as those with adjusted p-value < 0.05.

Prediction of miRNA-mRNA Interactions
Analysis of the interactions between miRNAs and mRNAs was applied to determine the possible functional role of miRNAs, which are differentially expressed in USP8mut and USP8wt tumors. We used both mRNA target prediction and correlation analysis of the expression levels of particular miRNAs and their predicted mRNA targets in corticotroph tumor samples.
The MicroRNA Data Integration Portal (mirDIP) algorithm, which combines multiple sources for miRNA target prediction [23], was used for the identification of possible mRNA targets. Only mRNAs that were predicted as targets with a probability scored as VeryHigh, according to the mirDIP criterion, were taken into account and included into downstream analyses.
Then, the correlation between the expression levels of identified potentially interacting miRNAs and mRNAs was assessed using normalized read-count data from small RNA sequencing and matched data from gene expression profiling for the same tumor samples (GSE132982 dataset). Spearman rank correlation was calculated using data from expression profiling of 24 tumors in R environment. Unadjusted p < 0.01 was considered relevant for correlation results.

Statistical Analysis
Datasets of quantitative variables were tested for the normal distribution with Shapiro-Wilk test. Variables with normal distribution were analyzed with two-sided unpaired t-test, while a two-sided Mann-Whitney U-test was used when normal distribution was not verified. Exact Fisher's test was used for the analysis of proportions. Significance threshold of α = 0.05 was adopted. For the identification of differentially expressed miRNAs and genes, p-values were adjusted for multiple testing with the Benjamini-Hochberg method. The Spearman correlation method was used for correlation analysis. Data were analyzed using GraphPad Prism 6.07 (GraphPad Software). Hierarchical clustering analysis was conducted with Cluster 3.0, and the results were visualized using TreeView 1.6 software (Stanford University School of Medicine, Stanford, CA, USA).

USP8 and USP48 Mutations
The incidence of hotspot mutations in USP8 and USP48 was determined with Sanger sequencing in 48 patients. USP8 mutations were identified in 15/48 (35.4%) patients. USP48 mutations were identified in 2/48 (4.2%) patients; both were female patients suffering from Cushing's disease with a diagnosis of densely granulated corticotroph tumors. One of the tumors was microadenoma while the other was macroadenoma. Both USP48 mutations had p.Met415Ile substitution. USP8 and USP48 mutations were mutually exclusive. Details of the results are presented in Table S1. The possible relationship between the incidence of USP8 mutations and demographic/clinical parameters was investigated in groups of Cushing's disease patients and SCA patients separately. Since only two patients with USP48 mutation were identified, they were excluded from the analysis. No difference in age of diagnosis was observed between USP8-mutated (USP8mut) and USP8-wild type (USP8wt) patients, both in the Cushing's disease group and in SCA group. Except for one male patient, all USP8 mutations were identified in females; however, differences of proportions did not reach statistical significance in the Cushing's disease group or SCA group.
All USP8mut patients suffering from Cushing's disease were in clinical remission after surgery, while clinical remission was observed in 9/15 of USP8wt patients. Difference of proportions of patients with/without remission was significant (11/0 vs. 9/6; p = 0.0237). A trend towards lower tumor volume was observed in USP8mut patients vs. USP8wt patients in both the Cushing's disease group (median 445 mm 3 vs. 2730 mm 3 , respectively; p = 0.0798) and SCA group (median 1844 mm 3 vs. 3893 mm 3 , respectively; p = 0.1707), but no significant difference was observed. Patients suffering from Cushing's disease stratified according to USP8 mutations status did not differ in terms of preoperative clinical parameters: morning serum ACTH level, morning serum cortisol level, or 24 h UFC.
Among patients with silent corticotroph tumors, significantly higher 24 h UFC level was observed in patients with USP8 mutations than in USP8wt patients (median 124.4 vs. 66.32, respectively; p = 0.0334). No difference in morning serum ACTH level, morning serum cortisol level, or midnight serum cortisol level were observed between these patients. We did not observe any difference between USP8mut and USP8wt patients in invasive growth status as determined with Knosp grading, proliferation index, or histological subtype (sparsely vs. densely granulated) either in the group of Cushing's disease patients or in those with silent tumors. The results are presented in Table 2.

Comparing miRNA Expression in USP8mut and USP8wt Corticotroph Tumors
The entire collection of 48 corticotroph tumors was subjected to miRNA expression profiling with next-generation sequencing of a small RNA fraction. Sequencing of small RNA libraries produced approximately 2,497,367 reads per sample, which were mapped to the human genome (hg19) and used for quantification of expression levels of known miRNAs, according to miRBase 22 release. Sequencing reads were annotated to 1917 miRNAs. Measurements of 1902 mature miRNAs were included in the analysis, after filtering out the miRNAs with low expression.
The overall analysis of the entire dataset with Principal Component Analysis (PCA) and hierarchical clustering methods did not show a clear separation between the groups of tumor samples stratified according to the mutation status, which indicates that the differences are not as pronounced as previously observed differences in gene expression profiles of USP8mut and USP8wt tumors [1,2]. Principal components 1 and 2 are presented in Figure S1a, while a dendrogram showing similarity of miRNA expression between the samples is shown in Figure S1b. Forty-six tumor samples were used for identification miRNAs differentially expressed in USP8mut and USP8wt tumors. Two samples with USP48 mutations were excluded from differential analysis to avoid bias resulting from possibly different molecular features of these tumors. A total of 250 differentially expressed miRNAs were found (adjusted p-value < 0.05), including 68 miRNAs that met the criterion of |FC| > 2, as shown in Figure 1a,b. Most of them (57/68) were miRNAs with the expression higher in USP8mut than in wild-type tumors (Figure 1b).

Putative mRNA Targets for Differentially Expressed miRNAs
To identify the mRNA targets of 68 differentially expressed miRNAs, a two-step procedure was applied. First, miRNA-mRNA interactions were predicted with the use of an mirDIP tool [23], and subsequently, the correlation analysis of matched miRNA and gene expression profiles was applied. This analysis included 24 tumor samples. For 49 out of 68 miRNAs differentially expressed in USP8mut vs. USP8wt tumors, significant correlation with predicted target mRNA was observed. A total of 442 miRNA-mRNAs interactions were identified with a median of four putative target mRNAs per single miRNA particle (ranging from 1 to 38 target mRNAs). Mostly negative correlation between miRNA and gene expression was observed as found for 303 miRNA/mRNA pairs (range of Spearman R coefficient −0.575 to −0.7922; median −0.6182). Positive correlation was observed for 139 miRNA/mRNA pairs (range of Spearman R coefficient: 0.5753 to 0.8361 median: 0.6215). Results are presented in detail in Table S3.
These analyses indicated 400 putative target genes that were identified as regulated by differentially expressed miRNAs. For the evaluation of potential functional significance of these genes, a subsequent gene set enrichment analysis (GSEA) was applied with the use of three gene ontology catalogs: KEGG Pathways, Gene Ontology (GO) Molecular Function and GO Biological Processes. Four KEGG pathways were found as significantly enriched for the putative target genes, including "Ubiquitin mediated proteolysis" as the most significantly enriched (according to p-value). The analysis with GO Molecular Function showed three protein ubiquitination related pathways as being in the top nine significantly enriched functions: Ubiquitin-like protein ligase activity (GO:0061659), Ubiquitin protein ligase activity (GO:0061630), and ubiquitin-protein transferase activity (GO:0004842). GO Biological process database indicates two processes related to the regulation of transcriptional activity and protein ubiquitination (GO:0016567) as the three most significantly enriched process (Figure 2a). The details of GSEA results are presented in Table S3. The genes with a clear ubiquitination function, which are common for the ubiquitination-related processes and pathways, are listed in Table 3 with details of the miRNA/mRNA correlation analysis. Since a difference in interaction patterns between miRNAs with negative and positive miRNA-gene correlation were described previously [24], GSE analyses were also performed separately for putative gene targets, the expression of which is positively or negatively correlated with DEMs. This showed no significant overrepresentation for genes with positive miRNA-mRNA correlation, while we observed a clear enrichment of protein ubiquitination-related pathways and processes for genes characterized by negative miRNA-mRNA correlation. The analysis with GO Biological process indicated a proteasome-mediated ubiquitin-dependent protein catabolic process ( Figure 2b and in Table  S3).   Table 3. The list of protein ubiquitination-related genes regulated by differentially expressed miRNAs (according to target prediction and correlation analysis) that were commonly found in the significantly enriched process in GSE analysis. The results of the correlation analysis between miRNA and predicted target mRNA levels and the results of miRNA differential analysis of USP8mut and USP8wt tumors.

Difference in miRNA Profile and Differential Gene Expression
The difference in gene expression profiles between USP8mut and USP8wt tumors was determined using sequencing data for 24 tumor samples, which were also included in miRNA-mRNA correlation analysis. The results of differential analysis indicated 1648 DEGs that met the criterion of adjusted p < 0.05. In order to identify DEGs with expression differences resulting from distinct miRNA profile in corticotroph tumors, stratified according to USP8 mutation status, we compared the results of three analyses: differential gene expression analysis, differential miRNA expression analysis, and identification of putative miRNA-mRNA interaction. DEGs with a direction of expression fold change that is concordant with the sign of correlation and the expression change of the corresponding miRNA were considered relevant. In case of genes with a negative miRNA-mRNA correlation, we looked for genes with opposite fold change signs, while in the case of genes with positive correlation, we looked for genes with the same fold change signs.

4.
Genes with distinct expression levels in USP8mut and USP8wt tumors with an expression difference related to els of differentially expressed miRNAs.

Discussion
Hotspot mutations in the USP8 gene, encoding ubiquitin carboxyl-terminal hydrolase 8, are the most common driver changes in corticotroph PitNETs. They have been detected in 30-40% of patients in previous studies [4,[6][7][8], as also observed in our cohort.
The clinical relevance of USP8 mutations was examined in few previous studies [4,6,8,[10][11][12][13][14][15]. Mutations were identified predominantly in women [4,5,12,13,25] and more frequently in younger patients according to some reports [6,14]. We found mutations nearly exceptionally in women (with only one male patient); however, due to a high overrepresentation of female patients in the study group, we cannot conclude that there is a sex-related difference. We also did not observe a relationship between mutation and age of onset.
In patients with Cushing's disease, these mutations appear related to lower tumor size and clinical remission after surgery [8,12], which could allow to consider USP8 mutations a possible favorable prognosis marker. However, they are also related to a higher risk of recurrence [14]. In concordance with these previous observations, we found a significantly higher rate of clinical remission in USP8mut patients suffering from Cushing's disease and a trend of lower tumor volume in mutated patients; however, no statistically significant difference was determined. We did not find differences in preoperative 24-h urinary free cortisol, which was observed to be significantly higher in patients with mutations in some previous studies [6,14]. Moreover, we did not find the difference in the proportion of tumors with invasive growth reported in other published results [5].
The analysis of the role of USP8 mutations in silent corticotroph tumors has not been reported previously. The mutations are less frequent in this group of patients than in the case of Cushing's disease. We identified the mutations in 20% of patients, while the mutation rate was 10% in the other study by Castellnoum, which included 20 SCAs [9]. No SCA with mutation was found in a study that included 11 such patients [6]. Because only four mutated SCAs were included in the analysis, the results may only be treated as preliminary. The only clinical parameter that we found significantly different in SCA patients with and without mutation was 24 h UFC, which was significantly higher in USP8mut patients; however, it was still within the reference range.
The generalization of our results should be done cautiously, and we have to emphasize an important limitation of our clinical data analysis. The numbers of patients included in the analysis are probably too low to draw a firm conclusion and the group is not representative of the general population, especially in the case of patients suffering from Cushing's disease. Since the primary aim of our study was molecular profiling of tumor tissue, we intentionally preselected large tumors, which allowed us to have enough tissue for DNA/RNA isolation and successful molecular procedures.
From a molecular biology point of view, USP8 mutations cause deregulation of the protein polyubiquitination/deubiquitination balance and impair the normal proteasomal degradation process [16]. In corticotroph tumors, a sustained EGFR signaling was found to be a consequence of the mutation [4,5] and the USP8 changes probably affect normal turnover of many other proteins regulated by ubiquitin carboxyl-terminal hydrolase 8 [26]. The pleiotropic effect of the mutation is observed not only at the level of protein degradation, but it is also manifested at the level of gene expression [1,2]. In this study, we aimed to characterize the difference in miRNA profile in tumors with USP8 mutations and wild-type PitNETs.
Recently, a novel driver mutation in another deubiquitinating enzyme i.e., USP48, was found in patients negative for USP8 changes. USP48 pathogenic variants cause increased activity of encoded deubiquitinase against its substrates Gli1 and H2A [17]. Because USP48 mutations affect the processes of protein degradation similarly to mutations of USP8, we screened the tumor samples for both the mutations. Two samples with USP48 Met/Val variant were identified and excluded from differential miRNA analysis to avoid a bias resulting from similarities in pathogenic mechanism. Due to a very low number of USP48-mutated tumors, these patients were also excluded from the analysis of clinical data.
We found that USP8mut and USP8wt tumors differ in miRNA expression; however, the differences are less pronounced than previously reported differences in mRNA expression [1,2]. We did not observe a clear distinction of USP8mut and USP8wt tumors in overall analysis, including PCA or hierarchical clustering based on the entire set of miRNA sequencing data, while a pronounced difference in mRNA expression profiles between USP8mut and USP8wt was reported previously [1,2].
USP8mut and USP8wt tumors differ in levels of relatively small proportion of all miRNAs that were included in differential analysis (approximately 3.5%). Conversely, a much higher proportion of differentially expressed protein-coding mRNAs was identified in corticotroph tumors stratified according USP8 mutational status [1,2].
Since miRNAs play a regulatory role in expression and translation [27], the consequence of differential miRNAs' expression depends on particular mRNA targets. To determine the possible functional role of DEMs in a high-throughput approach covering multiple miRNA-mRNA interactions, we predicted the putative mRNA targets for each DEM, followed by calculating the correlation coefficient for the expression levels of matched miRNAs and mRNAs. With this procedure, we mostly identified miRNA-mRNA pairs with negative correlation between expression levels where a high level of a particular miRNA corresponded with a low expression of its target gene. This relationship is concordant with a generally accepted concept that miRNAs are negative regulators of gene expression. Still, a notable part of the identified putative target mRNAs showed positive miRNA-mRNA correlation, indicating an activating role of miRNA. Activating action was previously reported for many miRNAs [28]. Recently published pan-cancer analyses [24,29] reported that many of miRNAs dysregulated in human cancer are positively correlated with their target genes.
GSE analysis was applied to identify the pathways where the identified target genes are overrepresented. The results showed that target genes of DEMs, especially those with negative miRNA-mRNA correlation, are related mainly to pathways and processes of protein ubiquitination. Since direct effects of USP8 mutations are the changes at protein ubiquitination level [4][5][6], we believe that the different expression of miRNAs that are involved in editing ubiquitin marks may reflect this major biological difference between USP8mut and USP8wt corticotroph PitNETs. Protein ubiquitination processes are directly orchestrated by a high number of proteins belonging mainly to three classes: ubiquitinactivating enzymes (E1), ubiquitin conjugating enzymes (E2), and ubiquitin ligases (E3) [30]. The enzymes of each class catalyze the subsequent stages of protein ubiquitination and the reverse reaction of protein deubiquitination is conducted by deubiquitinating enzymes [30]. Each class of the enzymes involved in editing ubiquitin marks includes multiple proteins, and genes encoding for proteins belonging to each class were found as putative targets of DEMs.
Distinct expression of miRNAs appears to have a limited effect on differential gene expression in USP8mut and USP8wt tumors. Less than 10% of predicted target genes with correlated miRNA-mRNA expression levels have significantly different expression in corticotroph tumors with and without mutation. The mRNA level of only 25 out of over 1600 DEGs could be considered as related to a different miRNA expression. This means that factors other than miRNA are responsible for the previously described highly different gene expression profile in USP8mut and USP8wt PitNETs.
None of the ubiquitination processes-related genes that were identified as putative targets of DEMs have significantly distinct expression levels. However, some of DEGs that were identified as targets of DEMs may have an interesting role in the biology of USP8mut corticotroph PitNETs. For example, our results indicate that hsa-mir-382-5p, which has higher expression in USP8mut tumors, may regulate genes involved in transcriptional regulation: KMT2C, KDM5A, and AFF1. KMT2C encodes for lysine methyltransferase that introduces mono-methylation mark at histone H3K4 [31], while KDM5A is lysine demethylase that converts di-and trimethylated H3K4 into mono-methylated form [32]. H3K4 mono-methylation plays a role in regulation of gene enhancers activation [33]. In turn, AFF1 functions as a regulator of transcription elongation and chromatin remodeling [34]. This suggests that hsa-mir-382-5p may contribute to a large difference in gene expression levels between USP8mut and USP8wt corticotrophinomas.
Our results also suggest that hsa-miR-655-3p, which has higher expression in mutated tumors, may affect expression levels of GAS1. The protein encoded by this gene is a negative regulator of the cell cycle and is considered a tumor suppressor in gastric and colorectal cancer [35,36]. It is known that corticotroph PitNETs with a distinct USP8 status differ in the expression status of cell cycle regulators at gene and protein level [1,2,15].
It is worth emphasizing that miRNA-mRNA interaction analysis results based on target prediction and calculation of expression correlation should be treated as preliminary. The commonly used methodology of detail validation of miRNA-mRNA interactions utilize laborious experimenting in vitro to confirm the impact of miRNA level on target gene expression and confirmation of physical miRNA-mRNA interaction with luciferase assay. This wet-lab approach is practically unfeasible for simultaneous investigation of many target genes of multiple miRNAs that we attempted to perform in our study. Additionally, this approach requires an appropriate cell model, but no human cell line of corticotroph cells is available. The only stable line of corticotroph cells are mouse AtT-20 cells and its usefulness in investigation of miRNA-mRNA interaction in human is limited due to evolutionary differences between species [37]. Some data on miRNA function in corticotroph cells based on mouse cell line were published [38][39][40]; however, it must be taken into account that approximately 46% miRNAs are considered primate-specific, while 14% are human-specific [37].
In summary, in our study we compared miRNA profiles of USP8mut and USP8wt corticotroph PitNETs and determined miRNAs with different expression levels. With target prediction and comprehensive miRNA and mRNA expression data analysis, we found that putative targets of DEMs are mainly the genes involved in processes and pathways of protein ubiquitination. However, differences in only a few miRNAs appear to affect the levels of genes with significantly diverse expression in corticotrophinomas with and without USP8 mutations. Thus, the difference in miRNA levels is not the cause of a pronounced differences in the gene expression between these tumors.

Supplementary Materials:
The following are available online at https://www.mdpi.com/2077-038 3/10/3/375/s1. Table S1: Detail patients' characteristics, Table S2: MiRNAs differentially expressed in USP8mut and USP8wt tumors, Table S3: Identification of putative mRNA targets of miRNAs differentially expressed in USP8 mutated and wild type tumors by target prediction and correlation analysis, Table S4: Top 10 pathways enriched for the genes that are targets of differentially expressed miRNAs (based on in silico target prediction and miRNA/mRNA correlation analysis); Figure S1: The overall analysis of corticotroph tumor samples stratified according to USP8 mutation status based on the entire miRNA dataset. a) Principal Component Analysis; b) Hierarchical clustering.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.