{'fraction_non_alphanumeric': 0.028190806853347423, 'fraction_numerical': 0.01762397162049966, 'mean_word_length': 5.402271079971007, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deduction of Novel Genes Potentially Involved in Upper Tract Urothelial Carcinoma Using Next-Generation Sequencing and Bioinformatics Approaches

Upper tract urothelial carcinoma (UTUC) is a relatively uncommon cancer worldwide, however it accounts for approximately 30% of urothelial cancer in the Taiwanese population. The aim of the current study is to identify differential molecular signatures and novel miRNA regulations in UTUC, using next-generation sequencing and bioinformatics approaches. Two pairs of UTUC tumor and non-tumor tissues were collected during surgical resection, and RNAs extracted for deep sequencing. There were 317 differentially expressed genes identified in UTUC tissues, and the systematic bioinformatics analyses indicated dysregulated genes were enriched in biological processes related to aberration in cell cycle and matrisome-related genes. Additionally, 15 candidate genes with potential miRNA-mRNA interactions were identified. Using the clinical outcome prediction database, low expression of SLIT3 was found to be a prognostic predictor of poor survival in urothelial cancer, and a novel miRNA, miR-34a-5p, was a potential regulator of SLIT3, which may infer the potential role of miR-34a-5p-SLIT3 regulation in the altered tumor microenvironment in UTUC. Our findings suggested novel miRNA target with SLIT3 regulation exerts potential prognostic value in UTUC, and future investigation is necessary to explore the role of SLIT3 in the tumor development and progression of UTUC.


Introduction
Urothelial carcinoma (UC), arising from the urothelium of the urinary tract, including bladder, ureter and renal pelvis, is the most common cancer among genitourinary tract cancer types. Common risk factors of UC include cigarette smoking and exposure to aristolochic acid and arsenic [1,2]. Patients with chronic kidney disease and end-stage renal disease also have higher incidence of UC in the Taiwanese population [3][4][5]. According to the anatomical location of tumor, UC can be classified into lower tract UC and upper tract UC (UTUC). The differences in clinical, demographic and molecular features between bladder Ivyspring International Publisher UC and UTUC have been reported. For instance, UTUC tends to have higher stage and grade than bladder UC due to thinner smooth muscle layer in the renal pelvis, and certain risk factors have larger impact on UTUC than on bladder UC [6,7]. Worldwide, bladder UC accounts for most of the UC, while UTUC is approximately 5%, and 2/3 of the UTUC occur in the renal pelvis [8]. In the Taiwanese population, however, UTUC accounts for approximately 30% of UC, and a slight female predominance was reported, differing from that reported in the Western countries [7][8][9]. The distinct epidemiology suggests potential endemic and molecular characteristics among our UTUC population.
Although studies suggest that UTUC shares many similarities with bladder UC, and similar genetic alterations regarding cell cycle and proliferative tissue markers have been reported [10,11], distinct genetic and epigenetic differences and mutation frequencies between UTUC and bladder UC exist [12][13][14]. Other than clinical and pathological characteristics, the prognostic value of tissue-based molecular biomarkers in UTUC has evolved rapidly [15]. In addition to alteration in tumor cell genetics, the critical roles of non-tumor cells, adjacent stroma, extracellular matrix (ECM) and altered tumor microenvironment (TME) have gained much attention [16][17][18][19]. The emerging biomarkers and novel therapeutic targets shed light on the advance in cancer treatment and importance of precision oncology [20].
Taking into account the distinct genetic and epigenetic differences in UTUC, recent progress in high-throughput next-generation sequencing (NGS) of the whole genome can achieve good resolution in characterization of genome-wide variations and facilitate the advance of precision oncology [21]. The rapidly progressing NGS technologies and development of powerful bioinformatics tools has made large genomic studies and the discovery of novel oncotargets more feasible [22]. In the current study, we aimed to investigate the distinct molecular signatures and novel miRNA regulations in UTUC, combining the NGS technique and bioinformatics approaches. We hope to identify novel targets of clinical significance and potential prognostic value in patients with UTUC.

Materials and Methods
The aim of our current study was to identify differentially expressed genes between tumor part and non-tumor part of UTUC clinical specimen through deep sequencing and identify novel microRNAs potentially involved in UTUC through bioinformatics approaches. The study flowchart is illustrated in Figure 1.

Clinical specimen
The two pairs of tumor and non-tumor tissue specimens were obtained from two female patients with renal pelvis UC during surgical resection. The specimens were collected within 30 minutes after radical dissection, and immediately stored in liquid nitrogen container to ensure the quality of tissue preservation. The detailed clinical background of the two patients was listed in Table 1. The study was approved by the Institutional Review Board of our hospital (KMUHIRB-E(I)-20170018).

Figure 1. Flowchart of study design.
Clinical specimens were obtained from two patients of upper urinary tract urothelial carcinoma (UTUC) for RNA and small RNA deep sequencing. The differentially expressed genes between tumor and non-tumor tissues were selected for enrichment analyses using various bioinforatmics databases. Furthermore, putative targets of selected differentially expressed microRNAs were predicted by miRmap database. The expression pattern and outcome prediction of candidate genes were analyzed using Oncomine database and Prediction of Clinical Outcomes from Genomic Profiles (PRECOG) database. Candidate genes with significant prognostic prediction were then validated in related urothelial carcinoma arrays in the Gene Expression Omnibus (GEO) database.  [23].

Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources
Database for Annotation, Visualization and Integrated Discovery (DAVID) is one of the bioinformatics enrichment tools that integrates large public bioinformatics resources and provides powerful tools for enrichment analysis of large gene lists from genomic experiments or sequencing results. Through multiple pathway-mining tools within the database, researchers gain a general concept of the biological themes and concentrated biological networks among the gene list of interest [24]. The DAVID Bioinformatics Resources 6.8 version was used for enrichment analysis in this study.

Gene Set Enrichment Analysis (GSEA)
The Gene Set Enrichment Analysis (GSEA) tool extracts relevant biological functions of a given gene list through a computational method that determines if a pre-defined gene set of genes shows statistically significant difference between two states, such as tumor and non-tumor phenotypes, instead of single-gene analysis. A gene set is a group of genes that share common biological function or regulation. The GSEA also provides leading-edge subset analysis, which extracts a subset of genes in a gene set as the core that contributes mainly to the enrichment signal [25]. The GSEA Desktop v3.0 was used for analysis in this study.

Ingenuity Pathway Analysis (IPA)
The Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems Inc., Redwood City, CA, USA) provides search function and network building/ analysis of an uploaded data. In the "Core Analysis" result, IPA generates unique networks based on the highest-fold change in the uploaded data. A network of interest can be graphed and overlaid for canonical pathways or specific diseases and functions selected [26]. The upstream regulator analysis is also available, which identifies molecules upstream of the genes in a given data that potentially explain the altered gene expression [27].

MiRmap Database
The miRmap software library is an open-resource for the target prediction of a specific miRNA, and the repression strength of a miRNA target is indicated by the miRmap score, which was estimated through a comprehensive computational method. A higher miRmap score indicates higher repression strength [28]. In the current study, miRNA targets with miRmap scores ≥ 99.0 were selected for further analysis.

Oncomine Database
The Oncomine platform integrates more than 700 independent datasets, expert curated data, and standardized analysis. Users can select differential expression analysis for automatically computed differential expression profiles of a selected cancer type or subtype of interest. Raw data including clinical information of selected datasets can be extracted for further analysis [29]. Differential expression analysis results of the candidate genes in urothelial carcinoma and transitional cell carcinoma of bladder and renal pelvis origins were extracted in this study.

Prediction of Clinical Outcomes from Genomic Profiles (PRECOG)
PRECOG is a new resource integrating cancer gene expression profiles and clinical outcome data from public database. It contains approximately 30,000 expression profiles from various cancer expression datasets, and all data were curated according to related clinical parameters [30].

Gene Expression Omnibus (GEO)
The Gene Expression Omnibus (GEO) database is a publicly available resource created since 2000 that accumulates free-access high-throughput genomic datasets. A web-based tool with graphic gene expression and raw data extraction of the candidate gene for statistical analysis is also available [31,32].
Other than datasets available from the Oncomine database, urothelial carcinoma related arrays were searched in the GEO database, and a dataset (GSE32894) containing 308 bladder cancer samples of different stages and grades was selected for candidate gene analysis in this study.

Statistical Analysis
The between-group expression difference of candidate genes obtained from selected arrays of Oncomine and GEO databases were analyzed using student's t test or one-way ANOVA with Tukey test for post-hoc analysis. The IBM SPSS Statistics for Windows, version 19 (IBM Corp., Armonk, NY, USA) was used for statistical analysis. A p-value < 0.05 was determined as statistically significant between-group difference.

Identification of differentially expressed genes in UTUC
The sequencing results of differential expression pattern of the two UTUC specimen was plotted in Figure 2A. There were 326 significantly up-regulated genes and 834 significantly down-regulated genes in tumor part tissue of UTUC specimen from patient 1. In addition, 562 significantly up-regulated genes and 653 significantly down-regulated genes in tumor part tissue of UTUC specimen from patient 2 were identified. By overlapping these dysregulated genes from two pairs of clinical UTUC specimens, we identified 86 up-regulated genes and 231 down-regulated genes in tumor part tissues of UTUC patients ( Figure 2B).

The differentially expressed genes were involved in extracellular matrix organization and cell cycle related biological functions
To determine the biological functions involved in these 317 differentially expressed genes of UTUC specimen, these genes were uploaded into DAVID database for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis. The top 10 GO and KEGG terms were shown in Figure 3, indicating the involvement of dysregulated genes in ECM organization, cell adhesion, and cell cycle pathways.
The GSEA enrichment analysis was also performed for gene sets of hallmarks, canonical pathways, motif and oncogenic signatures. The gene sets enriched in UTUC tumor tissues included G2M checkpoint, E2F targets, mitotic spindle and cell cycle canonical pathway ( Figure 4A upper panel), whereas matrisome and ECM glycoprotein related canonical pathways and epithelial mesenchymal transition gene sets were enriched in UTUC non-tumor tissues ( Figure 4A lower panel). The expressions of genes in related gene sets were displayed as heat maps in Figure 4A. Additionally, the motif gene set analysis indicated nuclear factor Y (NFY) as transcriptional factor targeting the dysregulated genes in UTUC tumor tissues, and the oncogenic signature gene set analysis indicated the representative gene signatures in polycomb repressive complex 2 (PRC2)/enhancer of zeste homolog 2 (EZH2) ( Figure 4B).

Identification of candidate genes with potential miRNA regulations in UTUC
To explore differentially expressed miRNAs and candidate genes potentially involved in miRNA regulations, small RNA sequencing was simultaneously performed. There were total 80 dysregulated miRNAs identified. The putative targets of these 80 differentially expressed miRNAs with miRmap score ≥ 99.0 were obtained from miRmap database. The overlapping genes between miRNA putative targets and differentially expressed genes of our dataset were achieved by Venn diagram analysis. The heat maps with z-scores of the differentially expressed miRNAs and mRNAs along with Venn diagram were illustrated in Figure 5. A total of 14 down-regulated genes and 1 up-regulated gene were identified as potentially involved in miRNA regulations. The expression values and fold-changes of these 15 candidate genes in the two pairs of UTUC specimens were listed in Table 2.
Four datasets from Oncomine database containing specimens of normal bladder tissue and bladder UC tissue were selected for comparison of candidate gene expression patterns, including Lee, Dyrskjøt, Sanchez-Carbayo, and Blaveri datasets. The heat maps of these genes in each dataset were illustrated in Figure 6, indicating the similar molecular changes among different bladder cancer datasets. Additionally, searching for the histological term of transitional cell carcinoma for UC, a dataset (Jones renal dataset) containing 23 normal kidney tissues and 8 renal pelvis UC tissues was achieved. We therefore also compared the expression pattern of the candidate genes in this dataset. The heat maps of the candidate genes of our UTUC data and Jones renal dataset were shown in Figure 7. The Oncomine database analysis identified significantly downregulated genes, including LMOD1, PDE5A, IGFBP5, FAM107A, TNS1, NCALD, and SLIT3 (p value < 0.05) in the Jones renal dataset of renal pelvis UC, indicating the coinciding novel molecular signatures between bladder UC and UTUC.

SLIT3 as a potential prognostic biomarker in UTUC
To predict the prognostic value of these candidate genes, the PRECOG database was used to determine the meta-Z score of each of the following genes: LMOD1, PDE5A, IGFBP5, FAM107A, TNS1, NCALD, and SLIT3. The higher meta-Z scores were observed in SLIT3 (meta-Z score = -2.47) and FAM107A (meta-Z score = -1.24) for bladder cancer, and no data available for UTUC in the database. The Kaplan-Meier plots for SLIT3 expression in two datasets (GSE5287 and GSE13507) were extracted from the PRECOG database, as displayed in Figure  8A, indicating lower survival probability in bladder UC patients with low SLIT3 expression. The expression pattern of SLIT3 between superficial and infiltrating bladder UC in Oncomine database was further investigated. There were four datasets comparing the expression pattern of SLIT3 between superficial and infiltrating bladder UC available, including Lee, Dyrskjøt, Sanchez-Carbayo, and Stransky (superficial bladder UC=25, infiltrating bladder UC=32). The results indicated higher ranking of SLIT3 under-expression in infiltrating bladder UC across each dataset ( Figure 8B).
We also searched in the GEO database for related UC arrays, and a bladder cancer array (GSE32894) containing 308 samples was selected for analysis of SLIT3 expression among different tumor stages and tumor grades. The expressions of SLIT3 were significantly lower in higher stages ( Figure 9A, upper panel) and grades ( Figure 9A, lower panel). The expression pattern of SLIT3 in the four datasets from Oncomine database also revealed lower expression of SLIT3 in infiltrating UC than in superficial UC ( Figure  9, B-E). The expression level of SLIT3 was also lower in renal pelvis UC than in normal kidney tissue ( Figure 9F).       Genomic Profiles (PRECOG) database was used for outcome prediction of candidate genes. Urothelial carcinoma patients with higher expression of SLIT3 had better survival rate. (B) Comparison between infiltrating and superficial urothelial carcinoma among four bladder cancer datasets from Oncomine database revealed higher ranking of SLIT3 under-expression in infiltrating type across each dataset. The heatmap analysis result was extracted from Oncomine database. The rank for a gene indicated the median rank for that gene across each analysis. The p-value for a specific gene indicated p-value for the median-ranked analysis. Red color represented ranking of over-expression genes, and blue color represented ranking of under-expression genes.

SLIT3 participates in cancer, renal and urological disease
To clarify the role of SLIT3 in dysregulated genes of UTUC, the 317 differentially expressed genes in UTUC were uploaded into IPA software for core analysis. The top three networks associated with differentially expressed genes in UTUC were listed in Table 3, where SLIT3 was clustered in network 1 related to diseases and functions of "Cancer, Organismal Injury and Abnormalities, Renal and Urological Disease". The hierarchical layout of network 1 was shown in Figure 10, with SLIT3 interconnecting to LMNB1. The overlay disease and function tool in IPA indicated the involvement of HGF, RRM2, TP63, KRT7, CDC6, MKI67, and SLIT3 in UC.

Discussion
The current study identified the differentially expressed genes in UTUC tissues were associated with aberration in cell cycle and ECM-related genes, analyzed by systematic bioinformatics approach. In addition, low tissue expression of SLIT3 in invasive UC was potentially a prognostic predictor of poor survival rate. Among the 15 candidate genes with potential miRNA-mRNA interactions, novel miR-34a-5p was a potential regulator of SLIT3, which may infer the potential role of miR-34a-5p-SLIT3 regulation in the altered TME in UC. A schematic figure summarizing the proposed molecular signatures of UTUC is displayed in Figure 11.
Alteration in cell cycle checkpoint pathways increases the risk of carcinogenesis [33]. The characteristic feature of mature urothelium is quiescent with low mitotic index and turnover rate, and the homeostasis of urothelial cell regeneration upon injury relies on epithelial-mesenchymal crosstalk, local secreted growth factors, and epigenetic regulation [34]. Similar molecular signatures of cell cycle and proliferative tissue markers between UC of bladder and upper urinary tract origin has been reported, and Ki-67, pRb, p53, and CDCA5 are of prognostic values [10,15,[35][36][37][38]. Using DAVID and GSEA for pathway enrichment analysis, MKI67, a gene encoding nuclear protein Ki-67, was identified as one of the top dysregulated genes related to cell proliferation and cell cycle checkpoint pathways in our UTUC dataset. In addition, enrichment in transcription factor NFY related genes was also identified in our UTUC dataset. NFY is a key regulator of cell proliferation, transcribing cell cycle regulatory genes, and the NFY targets are up-regulated in various cancer types [39].
A closer look into the dysregulated genes related to cell cycle checkpoint pathways, LMNB1 molecule was consistently involved, and identified in IPA database to be interconnected to SLIT3 (Figure 10). Amodeo et al. reported promyelocytic leukemia protein controls cell migration via PRC2-mediated SLIT repression in neoplastic brain cells [40]. Interestingly, the bioinformatics analysis also identified the dysregulated genes in our UTUC dataset were enriched in representative PRC2/EZH2 oncogenic gene signatures ( Figure 4B), with an average 7.1-fold increased expression of EZH2 in UTUC tumor tissue. PRC2 exerts epigenetic regulatory role during transcription and leads to target gene silencing [41], and ensures proliferative and regenerative potential of urothelial progenitor cells in response to injury [34]. EZH2 is one of the polycomb group genes, and the overexpression of EZH2 occur in prostate cancer and bladder cancer, and is associated with poor outcome in high grade UTUC cohort [42][43][44][45]. Taken together, epigenetic regulation of PRC2/ EZH2 on cell cycle checkpoint related genes and downstream SLIT3 suppression may have potential contribution to cancer development. However, the role of this regulatory axis in UTUC pathogenesis remains to be elucidated. SLIT3 (Slit guidance ligand 3) is one of the SLIT gene family members encoding secreted glycoprotein, Slits, that is highly expressed in mammalian tissues, and serves as ligands for roundabout (Robo) receptor, propagating many downstream signaling responses related to kinases and microtubule cytoskeleton [46,47]. SLIT3 is a tumor suppressor gene, where decreased expression of SLIT3 was observed in many human cancers through hypermethylation of SLIT3 promoter region or synergistic inhibitory effect with intronic miR-218, which also revealed an average 8.6-fold decreased expression in our dataset, leading to tumor invasion and progression [47][48][49][50][51]. The Expression Atlas, a database including gene and transcript expression from microarray-and sequencing-based functional genomics experiments [52], for RNA-seq mRNA differential expression of SLIT3 also revealed low SLIT3 expression in breast cancer, skin cancer, small cell lung cancer, and glioblastoma tissues, while higher SLIT3 expression was observed in esophageal adenocarcinoma and Barrett's esophagus tissues. The role of SLIT3 in UC has not been reported in the literature. Using bioinformatics approach, we identified low expression level of SLIT3 in high grade, invasive bladder UC and renal pelvis UC, and poorer survival rate in bladder UC patients with lower SLIT3 expression. This implicated the potential role of SLIT3 in UC tumor progression.
The crosstalk between tumor and adjacent non-tumor cells can be mediated by exosomes containing miRNAs, and change the tumor microenvironment to an invasion-promoting environment [53]. The role of secreted miRNA exosome in muscle-invasive bladder cancer has recently been reported, indicating the potential role of liquid biopsy for urinary exosomes as molecular marker in UC [54,55]. We identified up-regulated miR-34a-5p was the potential miRNA regulating SLIT3. Literature supported miR-34a as a tumor suppressor miRNA targeting key regulators of cell cycle and apoptosis through p53-dependent and p53-independent pathways [56][57][58][59][60], being silenced due to CpG methylation of its promoter in many cancer cell lines, including UC [61]. High expression of miR-34a was also reported to have reduced risk of bladder cancer recurrence [62]. Several studies also proposed the promoting effect of miR-34a in chemoresistance [63,64]. In bladder cancer cell lines, miR-34a activation can reduce cancer stem cell properties and sensitize these cells to gemcitabine and cisplatin [65,66]. These results suggested miR-34a may exert both oncogenic and oncosuppressive properties in different tumors and tumor environments. The role of miR-34a in UTUC has not been reported. Our current findings identified an average of 2.8-fold increased expression of miR-34a-5p in UTUC tumor tissue compared to adjacent normal urothelial tissue. Whether miR-34a-5p-SLIT3 regulation exerts regulatory role in UTUC development or progression merits further investigation.
The ECM is a major component of TME that regulates cell behavior [67]. The study of ECM proteome, namely matrisome, in cancer is evolving, emphasizing the important roles of tumor surroundding microenvironment and other non-tumor cells in cancer [17,19]. The gene ontology enrichment analysis of our UTUC dataset identified the dysregulated genes in UTUC were most differentially enriched in cellular component of ECM region (135 of the 317 dysregulated genes), as shown in Figure 3C. Among the 135 genes, 57 molecules were ECM-associated, as identified from MatrisomeDB 2.0 introduced by Naba et al. [68,69], where SLIT3 was grouped as one of the core matrisome genes. Excessive ECM deposition is one of the characteristic features of cancer, and several core matrisome genes predicted cancer outcome in various cancer types [70]. However, our GSEA enrichment analysis results indicated gene sets related to matrisome were enriched in UTUC adjacent non-tumor tissue instead ( Figure 4A, lower panel). The concept of TME holds true that tumor cells react to altered TME and affect gene expression of normal tissue adjacent to the tumor (NAT), and unique characteristics of NAT from healthy tissue has recently been reported [71]. The mechanisms underlying altered gene expression in NAT remain to be validated, and changed microenvironment in the adjacent stroma may possibly be regulated by the tumor [72]. Our result may implicate the dysregulated genes in UTUC are potentially involved in altered ECM tumor microenvironment, while the mechanistic link between matrisome-associated genes and tumor related biological processes necessitate further investigation.
There are several potential weaknesses to be addressed. Firstly, the current findings of potential SLIT3 regulation and prognostic prediction for UTUC were based on two patients of UTUC and validated by systematic bioinformatics analyses. Further investigation and longitudinal follow-up in a larger UTUC cohort is necessary to confirm its prognostic value. Additionally, the UTUC specimens were collected from Taiwanese population, where epidemiology of UTUC was different from reported worldwide, suggesting the potential influence of environmental and/or genetic factors, thus may limit the generalizability of the current findings in different populations.

Conclusions
The current study identified novel miRNA target with potential SLIT3 regulation may have potential prognostic value in UTUC, possibly related to aberrant cell cycle progression and altered tumor microenvironment. Further investigation to confirm the role of SLIT3 and its miRNA regulation in UTUC tumor initiation and progression is of potential clinical significance.