Analysis of HSD11B2 as a prognostic marker in Melanoma via TCGA data mining

Background: Hydroxysteroid 11-Beta Dehydrogenase 2 (HSD11B2) expression has been reported to be present in melanoma. We investigated the association of HSD11B2 with melanoma using publicly available data from The Cancer Genome Atlas (TCGA). Methods: The relationship between clinical pathologic features and HSD11B2 were analysed via Wilcoxon signed-rank test and logistic regression. Clinicopathologic characteristics associated with overall survival in melanoma patients were calculated using Cox regression and the Kaplan-Meier method. Gene Set Enrichment Analysis (GSEA) and gene co-association of HSD11B2 were performed using TCGA data set. Results: Reduced HSD11B2 expression was signicantly lower in melanoma patients compared to normal patients (p value = 3.004e-122) and also associated with lower survivability. low HSD11B2 expression in melanoma was also signicantly associated with cancer stages T (p value = 0.002) and N (p value < 0.001) and age (p value = 0.003). Genes TFCP2L1, PRR15L, ATP6V1B1, C9orf152, AC009948.1, AL391244.1, WDCP, HNRNPCP2 and GTF2E1 were all shown to be co-associated with changes in HSD1B2 expression. Multiple signalling pathway including cytosolic DNA sensing pathway, JAK STAT signalling pathway, NOD like receptor signalling pathway, T cell receptor signalling pathway and Toll like receptor signalling pathway were differentially enriched in low HSD11B2 expression phenotype. Conclusion: Our study revealed that HSD11B2 expression is closely associated with melanoma development and age, as well as multiple cancer related genes and pathways, thus highlighting HSD11B2 as a potential therapeutic marker of melanoma.


Background
Melanoma is a type of cancer that occurs in the skin and on rare occasions in the mouth, intestines, or eye. Melanoma develops from the pigment-producing cells known as melanocytes located in the bottom layer of the skin's epidermis [1]. According to the World Health Organisation (WHO) melanoma is one of the most dangerous type of skin diseases [2] with a resulted death of 59,800 from 3.1 million cases in 2015 [3]. Although the incidence rate of melanoma is low, this cancer has a high degree of malignance and occurs mostly in adults. The development of melanoma is a contribution of genetical and environmental factors, or a combination of both. One of the main environmental causes of melanoma is exposure to ultraviolet light which results in DNA damage [4]. Indeed, countries such as New Zealand and Australia has the highest occurrence of skin melanoma due to higher and more direct exposure to ultraviolet light [2]. Rare mutations that are often inherited in families also greatly increases melanoma susceptibility. For example, mutation in the MC1R gene (a gene that causes red hair) are 2 to 4 times more likely to develop melanoma than those with two wild-type copies [5]. Other class of mutation such as alternative reading frame and acquired mutations can affect the genes such as CDKKN2A which leads to destabilization of p53, a transcription factor included in apoptosis and in 50% of human cancers [6].
Mutations in CDKKN2A can also results in a non-functional inhibitor of CDK4, a cyclin-dependent kinase that promotes cell division [7]. Mutations that cause the skin condition xeroderma pigmentosum (XP) also increase one's susceptibility to melanoma [8]. Because genetical factors play a key role in the occurrence of skin melanoma, the signi cant expansion of knowledge in the association of different genes with skin melanoma together with technical advances can create an opportunity to interrogate the genetic basis of cancer risk, response to therapy and outcome.
HSD11B2 gene encodes a human enzyme known as Corticosteroid 11-β-dehydrogenase isozyme 2 (11-βhydroxysteroid dehydrogenase) that is expressed in aldosterone-selective epithelial tissues such as the kidney, colon, salivary and sweat glands. HSD11B2 is involved in the conversion of 11-hydroxy cortisol to cortisone and modulates intracellular glucocorticoid (a vital lipid hormones regulator in stress responses, metabolism and immune homeostasis) levels, thus protecting the nonselective mineralocorticoid receptor from occupation by glucocorticoid. HSD11B2 participates in cancer [9,10] and promotes cancer progression. HSD11B2 is reported to be essential for colon-and lung carcinogenesis by relieving GCmediated cyclooxygenase-2 (COX-2) inhibition [10,11]. In addition, abnormal expression of HSD11B2 has been shown to result in signi cantly promoted migration and metastasis of colorectal carcinoma cells in vitro and in vivo [12]. HSD11B2 activity has also been reported to play a role in tumour proliferation in ovarian cancer via eliminating anti-in ammatory action of cortisol. Despite the role of HSD11B2 being well reported in cancer physiology [13][14][15][16], little is known about its role and expression in melanoma.
Thus, the objective of the current study was aimed to evaluate the prognostic value of HSD11B2 expression in human melanoma data obtained from TCGA. Indeed, gene set enrichment analysis (GSEA) was performed in combination to multiple bioinformatics and statistical analysis to gain a better understanding of clinicopathologic characteristics that are related to HSD11B2 as well as gain insight into the underlying pathophysiological pathway mechanisms associated with melanoma pathogenesis.

RNA-sequencing patient data and bioinformatics analysis
Unlike liver and kidney tumors, it is di cult to obtain normal controls for melanoma. Therefore, the TCGA database regarded control cases as only one control. In order to perform statistical analysis, GTEx database was used as a supplementary control. High throughput sequencing of gene expression data (HTSeq-counts) and clinical information of 323 control cases and 471 cases of melanoma cases were downloaded from TCGA o cial website. Expression of HSD11B2 differences for discrete variables were visualized using boxplots and whiskers plot [17].
Gene set enrichment analysis GSEA is bioinformatics method aimed to identify whether prior sets of genes or proteins are signi cantly different between two phenotypes [18]. Our study applied GSEA to generate an order list of all genes according to their correlation with HSD11B2 expression, signi cant survival differences observed between high and low HSD11B2 groups were elucidated. Gene set permutations were performed 1000 times for each analysis. The expression level of HSD11B2 was used as a phenotype label. The nominal p value and normalized enrichment score (NES) were used to sort the pathways enriched in each phenotype. In addition, to investigate associated genes in performing different molecular function and biological pathway, gene correlation analysis was performed for the HSD11B2 gene. Top 9 signi cant genes (p value < 0.05) that are correlated with expression alteration were ranked based on their Pearson's correlation factor. R (v3.4.3) was used to construct a circus plot displaying these Pearson's correlation between these genes and HSD11B2.
Statistical analysis R (v3.4.3) was used to perform all statistical related analysis. Relationship between clinical pathological features and HSD11B2 expression were analysed via Wilcoxon signed-rank test and logistic regression. Nomogram construction was performed according to the guidelines proposed by Iasonos [19]. We randomly assigned two-thirds of patient information to the development of the cohort model (n = 314) and one-third of patient information to the validation cohort model (n = 157).
To identify independent prognostic predictors, we used a Cox proportional hazards regression model for univariable and multivariable analyses by the "Enter" method. The nomogram was developed to predict the 1-, 3-and 5-year prognosis mainly based on the results of the multivariable Cox regression model.
The performance of the nomogram was estimated regarding discrimination and calibration. The C-index was applied to evaluate discrimination [20], which refers to the models' ability to accurately distinguish the outcomes. A higher C-index indicates more precise model predictions [21]. Calibration curves were performed by comparing the means of the nomogram-predicted outcomes with the actual outcomes estimated with Kaplan-Meier. The bootstrapping (1000 repetitions) method was applied to reduce the estimate bias. In addition, model validations were performed using the data of the validation melanoma cases as follows. First, we calculated the total points of the patients in the validation group using the established nomogram. Next, we used the total points as a factor to perform Cox regression analysis.
Finally, the C-index and calibration curves were developed with the results of regression analysis. Receiver operating characteristics (ROCs) curve was used for the sensitivity and speci city of nomogram.

HSD11B2 expression
Box whiskers plot (Fig. 1A) reveals that HSD11B2 is similarly expressed in the skin of both male and female. In addition, expression of HSD11B2 is signi cantly lower in patients with tumour (p value = 3.004e-122) (Fig. 1B). The Kaplan-Meier survival analysis (Fig. 1C) showed that melanoma patients with low expression of HSD11B2 had lower survival probability compared to patients with high expression of HSD11B2 (p value = 9e-5).

Survival outcomes and multivariate analysis
The univariate analysis revealed that age (p value = 0.003), size or direct extent of primary tumour (T) (p value = 0.002), degree of spread to regional lymph nodes (p value < 0.001) and HSD11B2 expression (p value < 0.01) are the most highly correlated factor to poor overall melanoma survivability ( Fig. 2A). The multivariate analysis revealed similar results including age, T, N and HSD11B2 expression as associated with overall survival (Fig. 2B). Classical univariate ROC curve analyses showed that HSD11B2 expression, T and N had reasonably high AUC of 0.67, 0.70 and 0.67 respectively (Fig. 2C). A point ranking system was also developed to rank the association of each factor with survivability (Fig. 2D). The result represented via the nomogram that the higher the point for a given factor, the lower the survivability. Stages T (T0, T1, T2, T3, T4, TX) and higher HSD11B2 expression level are major point contributors to low survivability (1year-, 2-years-and 5 years-). In addition, higher age, higher stages of N and the Asian ethic group were also major point contributors to low survivability. Finally, AUC curves on both the predicted and validated nomograph model veri es (all AUC values > 0.7) the 1 year-, 2-years-and 5-years survivability prediction and highlights its statistical signi cance (Fig. 3).

GSEA identi es HSD11B2-related pathways
To identify pathways that are differentially activated in melanoma, we conducted Gene Set Enrichment Analysis (GSEA) between low and high HSD11B2 expression data sets. GSEA reveal signi cant differences (FDR > 0.05, NOM p-val > 0.05) in enrichment of MSigDB Collection (c2.cp.biocarta and h.all. v6.1. symbols). We selected the most signi cantly enriched signalling pathways based on their normalized enrichment score (NES) (Fig. 5). The results showed that cytosolic DNA sensing pathway, JAK STAT signalling pathway, NOD like receptor signalling pathway, T cell receptor signalling pathway and Toll like receptor signalling pathway are differentially enriched in high HSD11B2 expression phenotype.

Discussion
HSD11B2 gene encodes the type II isoform of 11-beta-hydroxysteroid dehydrogenase, a microsomal enzyme complex responsible for the interconversion of biologically active cortisol and inactive cortisone. This gene is highly expressed in kidney (highest expression) [22], colon [23], pancreas [24] and placenta [25]. According to the enzyme level expression pro le of HSD11B2 from Protein Atlas (http://www.proteinatlas.org), HSD11B2 is downregulated in human skin squamous cell carcinoma and melanoma as well as in eight other cancer types [26]. To our knowledge, the underlying pathophysiological mechanism and clinicopathologic characteristics of HSD11B2 in melanoma has not been closely examined. Therefore, this study applied bioinformatics analysis using high-throughput RNAsequencing data from TCGA to examine HSD11B2 expression in melanoma patients to demonstrate that low expression of HSD11B2 is mostly associated with direct extent of the primary tumour, degree of spread to regional lymph nodes, and age of melanoma patients. All these contributing factors signi cantly results in a low survivability. In addition, we discovered that changes HSD11B expression is associated with numerous receptors signalling pathways and multiple genes.
Anatomical extent of melanoma tumour and degree of spread to regional lymph nodes alongside increase of HSD11B2 gene expression are major contributing factors in lower predicted survivability of melanoma. Indeed, the thickness of the primary tumour is an important prognostic factor in predicting the risk that the cancer will spread. Melanoma skin cancer with tumour size less than 1 mm thick has a low chance of spreading, whereas tumour with more than 4 mm in thickness has a higher risk of spreading and recurrence [27]. Furthermore, prognosis is poorer if the cancer has spread to nearby lymph nodes which is an indication of metastasis stage which directly contributes to survivability. Interestingly, The Asian ethnic group has a higher contributing factor in lower survival rate. Though people of colour are less likely to became a icted with melanoma, they are more likely to die from it due to delay in detection and presentation. This may be associated with lack of awareness, diagnosis at more advanced stage and socioeconomic barriers hindering access to care [28].
Our study revealed that patients with melanoma tumour has a lower HSD11B2 expression, as well as a lower predicted survivability. Type I and Type II 11-beta-hydroxysteroid dehydrogenase encoded by the HSD11B1 and HSD11B2 gene respectably governs the glucocorticoid action on target tissues.
Glucocorticoid are a class of steroid hormones that is vital in regulation of stress responses, metabolism and immune homeostasis [29]. HSD11B1, functions as a reductase and catalyses the regeneration of active glucocorticoids, thus amplifying cellular action. On the other hand, HSD11B2 is a high-a nity dehydrogenase that inactivates glucocorticoids [30]. Therefore, 11-beta-hydroxysteroid dehydrogenase are important enzymes in the tissue-speci c regulation of glucocorticoids, with reports suggesting that alteration in regulation of 11-beta-hydroxysteroid dehydrogenase is related to numerous pathological process in the skin such as wound healing, cell motility, scar formation and remodelling and promotion of terminal epidermal differentiation [31][32][33][34]. A study by Cirillo et al. (2017) demonstrated that HSD11B2 is synthesised by human malignant keratinocytes (epidermis cells that forms barriers against environmental damages such as UV radiation, heat, and viruses [26]. The same report also revealed that a signi cant downregulation of HSD11B2 was observed in human skin squamous cell carcinoma and demonstrated that skin cancer cells have the ability to regulate cortisol levels through expression of both HSD11B1 and HSD11B2 levels. Similarly, downregulation of HSD11B2 mRNA has also been reported in skin lesions of leprosy patients compared to normal skin [35]. Multiple types of cancers such as colorectal cancer [36] and oral cancer [33] also displays a trend of downregulated HSD11B2 expression. These evidences agree with our results and highlights the importance of HSD11B2 in the glucocorticoid system and its potential signi cance in melanoma. As shown by our results, changes in the expression of HSD11B2 is co-associated with nine other genes. Among these, ATP6V1B1 (positively correlated with HSD11B2) is a gene that encodes a component of vacuolar ATPase, an enzyme which mediates acidi cation of eukaryotic intracellular organelles [37]. This acidi cation in uences a wide range of cellular process, many of which are dysregulated in cancers [38]. Dysregulation of vacuolar ATPase and its association to numerous cancers was a conclusion made based on altered expression of speci c vacuolar ATPase subunits in malignant cells and clinical features of cancer. For Instance, a study concluded that the downregulation of vacuolar-ATPase subunit ATP6V1B1 leads to resistant mechanism of trastuzumab mediated antibody dependent cellular cytotoxicity in human breast cancer cell line (HER2 overexpressing) [39]. Analysis of the TCGA database also revealed that ampli cation or overexpression of ATP6VC1 (gene that codes for the C subunit of vacuolar ATPase) was observed in 34% of human breast cancer resulting in low survivability [40]. Other studies revealed expression of vacuolar ATPase speci c subunits in other forms of cancers such as pancreatic cancer [41], lung cancer [42], oral squamous cell carcinoma [43,44], ovarian cancer [45] and cervical cancer [46]. Therefore, it is likely that ATP61B1 also plays a role in melanoma and its association with HSD11B2 requires further investigation. Another gene that was positively associated to HSD11B2 expression changes was KRT19. KRT19 referred to Keratin 19, is a member of the acidic type 1 cytokeratin family protein that is expressed predominantly in various cells types in the epithelium [47][48][49]. Since 80% of cancers are of epithelial cells in origin [50,51], Keratins are widely used as diagnostic markers to detect tumours in both primary and distal sites and to determine the tissue of origin of the tumour for the purpose of assisting in treatment strategies [52]. In this context, KRT19 has been suggested to be one of the more sensitive diagnostic markers across a broad range of cancer types. Indeed, studies indicated that KRT19 is highly expressed in breast, colon, liver, and intestine cancer, and also associated with poor clinical outcomes in patients [53][54][55][56]. Since activated keratinocytes expresses various forms of keratin proteins (KRT6 KRT16 and KRT17), the association of KRT19 expression to HSD11B2 is an underlying molecular mechanism that needs to be examined in future melanoma researches.
Numerous pathways were revealed to perturbate between HSD11B2 low expression and high expression phenotype via the GESA pathway analysis of TCGA data. Among these pathways, the cytosolic DNA sensing pathway (cGAS-STING pathway) is part of the innate immune system which recognises presence of cytosolic DNA, trigger expression of in ammatory genes (leading to senescence) and activates defence mechanisms [57]. Numerous studies have indicated that tumour cells developed strategies to inhibit DNA sensing pathway activation, likely for immune evasion during carcinogenesis in colorectal carcinoma and melanoma [58,59]. The JAK-STAT signalling pathway is another pathway closely related in the process such as immunity, cell death and tumour formation. Disruption in JAK-STAT signalling may lead skin conditions, cancers, and disputed immune systems [60]. Indeed, a study by Pansky et al. (2000) hypothesizes that defective signal transduction through the JAK-STAT pathway could contribute to the poor response rate of interferons (antiviral, antiproliferative and immunomodulatory protein) treatment in patients with advanced melanoma [61]. Both the NOD like receptor signalling pathway and toll-like receptor signalling pathway were detected to be differentially expressed in relationship to changes in HSD11B2. These pathways belongs to families of the innate immune system. Activation of these immune pathways leads to a broad range of pro-and/or antiin ammatory signals, including the secretion of interferons, tumour necrosis factors and cytokines [62]. Disruption in these signals may result in chronic in ammatory states that directly affect cell cycle progression and apoptosis which subsequently creates an environment for diseases, such as cancer [62,63]. Based on the presented evidence, all the identi ed pathways revealed by our results plays major roles in cancer and highlights the importance of their association with HSD11B2 expressional changes in melanoma for future researches.

Conclusion
HSD11B2 expression may be a potential prognostic molecular marker of poor survival in melanoma.
Changes HSD11B2 expression in melanoma patients is closely related with the size of the primary tumour, degree of spread to regional lymph nodes and age which in turn affects survival outcomes. The HSD11B2 gene is also closely co-associated with multiple other genes and signalling pathways that plays critical roles in various cancer. Our study highlights HSD11B2 as a possible target for therapeutic target in melanoma cases and further experimental validations.

Declarations
Ethics approval and consent to participate The research was approved by the Ethics Committee of Zhuhai People's Hospital.

Consent for publication
Not applicable.

Availability of data and materials
All data generated or analyzed during this study are included in this published article.

Funding
Not applicable.

Competing interests
The authors declare that they had no competing interests.
Authors' contributions Qixian Chen wrote the paper, performed the experiments and analyzed the data, Zhaohui Liu , Lintao Zhong, Zhike Liang performed the experiments and analyzed the data, Haoyu Song and Lijun Chen analyzed the data and wrote the paper, Jun Li and Jian Yang designed the research, Yuqing Weng designed the research and wrote the paper. All authors read and approved the nal manuscript.     Plot displaying the result of GESA. Multiple pathways were shown to have signi cant changes in enrichment scores in relation to high expression verses low expression of HSD11B2 gene.