BIRC5 and PRLR Might Be Novel Immune-Related Prognostic Biomarkers in Adrenocortical Carcinoma (ACC)

Background: Adrenocortical carcinoma (ACC) is a highly aggressive tumor(cid:0)which occurs in adrenal cortex. This study aims to identify new immune-related molecular biological markers and therapeutic targets in ACC patients. Results: First, we identied 134 genes by constructing weighted gene co-expression network. Among them, two genes (BIRC5 and PRLR) were eventually regarded as immune-related prognosis biomarkers on the basis of The Immunology Database and Analysis Portal (ImmPort) database. The validation of mRNA level was meaningful (based on Oncomine database and TCGA-ACC). Four methods (overall survival (OS) analysis, disease-free survival (DFS) analysis, stage plot and expression value comparison) were used in GEPIA to verify that PRLR and BIRC5 have signicant correlation with the prognosis of ACC. Several other methods including survival analysis, analysis of variance, the spearman correlation analysis, distance correlation analysis and receiver operating characteristic curve were also used for the external validation of them in GEO datasets (GSE10279, GSE19750, GSE76019 and GSE76021). Furthermore, we evaluated the relationship between hub gene expression of and the degree of immune cell inltration in ACC by ESTIMATE algorithm and ssGSEA. BIRC5 expression was signicantly positively related to activated CD4 T cell and type 2 T helper cell. And PRLR expression was positively related to CD56dim NK cell and type 17 T helper cell. Moreover, Gene set enrichment analysis (GSEA) suggested that PRLR has strong correlation with neuroactive ligand receptor interaction and BIRC5 was signicantly enriched in cell cycle. In addition, according to Drug Signatures Database (DSigDB), BIRC5 were signicantly enriched in “resveratrol_mcf7_down”, “lucanthone_ctd_00006227”, “8-azaguanine_pc3_down”, “thioguanosine _mcf7_down”, “dasatinib_ctd_00004330”, “monobenzone_pc3_down”. Conclusions: In summary, BIRC5 and PRLR were

Furthermore, we evaluated the relationship between hub gene expression of and the degree of immune cell in ltration in ACC by ESTIMATE algorithm and ssGSEA. BIRC5 expression was signi cantly positively related to activated CD4 T cell and type 2 T helper cell. And PRLR expression was positively related to CD56dim NK cell and type 17 T helper cell. Moreover, Gene set enrichment analysis (GSEA) suggested that PRLR has strong correlation with neuroactive ligand receptor interaction and BIRC5 was signi cantly enriched in cell cycle. In addition, according to Drug Signatures Database (DSigDB), BIRC5 were signi cantly enriched in "resveratrol_mcf7_down", "lucanthone_ctd_00006227", "8-azaguanine_pc3_down", "thioguanosine _mcf7_down", "dasatinib_ctd_00004330", "monobenzone_pc3_down".
Conclusions: In summary, BIRC5 and PRLR were shown to be related to the prognosis of patients with ACC, indicating them as novel immune-related prognostic biomarkers and immunotherapy targets in ACC.

Background
Adrenocortical carcinoma (ACC) is a highly aggressive tumor,which occurs in adrenal cortex. The overall 5-year survival rate is 35% due to highly aggressive biological behavior (1). It mostly occurs between the ages 40 and 60, with higher prevalence in female (up to 60%) (2). Since ACC has few symptoms in the early stage, many patients have already advanced in the late stage when they are diagnosed (3). At these stages, ACC has extensive invasiveness and metastasis with limited treatment options. Actually, the current treatments for ACC, including surgery, chemotherapy and radiotherapy, shows poor e cacy and prognosis (4). Under current circumstance, it is urgent to study the mechanism that promotes the progress of ACC and nd new molecular biological markers to predict the prognosis of ACC.
Many scholars are dedicated to studying the impact of the immune system on tumor progression. In recent years, accumulating studies indicated that the tumor immune microenvironment (TIM) plays a vital role in tumor progression, invasion and metastasis (5). The in ltration of various immune cells might affect the proliferation and metastasis of tumor cells(6). Papewalis C et al demonstrated that autologous DCs may serve as adjuvants for immunotherapy in ACC (7). Kanczkowski W et al revealed the signi cant reduction in the expression and signal transduction of TLR4 and CD14 in malignant but not benign adrenal tumors, which is a novel feature of adrenocortical tumors and may be used as diagnostic markers(8). Parise IZS et al showed that CD8 + T Lymphocytes could be the prognostic markers for ACC (9). To date, the development of immunotherapy for ACC was rarely documented because of the uncertainty in the relationship between tumor cells and the tumor microenvironment (TME). Therefore, in our study, we aimed to screen out novel diagnostic and prognostic biomarkers related to immune for ACC, which can also be used as therapeutic targets for immunotherapy meanwhile.
Based on the facts above, we identi ed modules associated with prognosis by weighted gene coexpression network. Then the hub genes in the key modules were selected. Following, we screened out two immune-related genes (baculoviral IAP repeat containing 5 (BIRC5) and prolactin receptor (PRLR)) from them. Then, we conducted internal and external validation, respectively, to further verify our results. It is worth mentioning that we regard survival time and survival status as important prognostic yardstick, so as to screen out prognostic-related modules. Previous studies chosen prognosis-related modules only by survival time or survival status. As far as we know, this study initially explore the potential relationship between immune cell in ltration and prognosis biomarkers by single-sample gene-set enrichment analysis (ssGSEA) in ACC. These two genes could be novel prognostic biomarkers and therapeutic targets in ACC. Moreover, based on Drug Signatures Database (DSigDB) (10), which is a new gene set resource that uses drugs/compounds and their target genes for gene set enrichment analysis (GSEA), we also studied the correlation between hub gene and six kinds of drugs in order to nd new targets of anticancer drugs.

Data Preprocessing
The expression matrix was displayed as count number. For subsequent analysis, the raw expression data of TCGA-ACC was preprocessed by R package "DEseq.2"(15), by which we took normalization and log2 transformation too. Moreover, normalization and log2 transformation of GEO datasets (GSE10927, GSE19750, GSE76019 and GSE76021) was conducted using the R package "affy"(16).

Co-expression network construction
Firstly, we use the gsg (goodSamplesGenes) method to check whether the top 5000 genes in the expression pro le are suitable for WGCNA. Moreover, we distinguished the outlying samples which were distant from other samples. All outliers were excluded from the ananlysis. Then, co-expression network was constructed by using the "WGCNA" R package (17). An appropriate β (soft threshold power beta) was selected according to the scale-free topology criterion. Next, we calculated topological overlap matrix (TOM) using the adjacency matrix to measure the network connectivity of genes. Then, we conducted three different branch cutting methods to divide genes into gene modules. The cut line of merging modules is based on the dissimilarity of module eigengenes (MEs).

Identi cation of Clinically Signi cant Modules
Firstly, the gene signi cance (GS) and module signi cance (MS, the average GS for genes in a module) were obtained to quantify the relationship between gene expression and clinical features. In general,, the module which has the largest absolute MS was regarded as the module related to clinical feature. In addition, we calculated the module membership (MM) to assess the correlation between ME and gene expression.

Functional and Pathway Enrichment Analysis
We performed Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis pathway by using R package "clusterPro ler"(18) to demonstrate the mechanism of genes in key modules. P < 0.05 were considered to be cutoff.

Identi cation of hub genes
After the key modules were chosen, hub genes were de ned under the threshold of |MM|≥0.80 and |GS| ≥0.20. Then, in order to identify immune-related prognosis biomarkers, we regarded the intersection between the hub genes in WGCNA and IRGs in Immport database as immune-related prognostic biomarkers.

Validation of hub genes
Firstly, based on the Gene Expression Pro ling Interactive Analysis (GEPIA) database (http://gepia.cancerpku.cn/), the relationship between hub genes and prognosis was veri ed by four methods including survival analysis (overall survival (OS), disease-free survival (DFS)), stage plot and expression values comparison.
Then, four additional independent validation GEO datasets (GSE10927, GSE19750, GSE76019 and GSE76021) were obtained to further test the prognostic role of hub genes by six methods. The analysis of variance (ANOVA) test and the Spearman correlation analysis were carried out by using SPSS (version 21.0), which was visualized by R package "ggstatsplot" (19). Meanwhile, the distance correlation analysis was completed by R package "energy" (20). Survival analysis was carried out by R package "survival" (21).
We divide ACCs into high gene expression group and low gene expression group and generate Kaplan-Meier survival curve. In addition, by virtue of R package "plotROC"(22), we performed receiver operating characteristic curve (ROC) analysis and calculated the area under the curve (AUC) to measure the resolving power for distinguishing the adrenocortical tumor from normal adrenal glands. Finally, the mRNA-level expression of hub genes between ACC and normal tissues was compared on the basis of Oncomine database(23) (https://www.oncomine.org/).

Relationship between gene expression and TME
We performed ssGSEA (24) to assess the in ltration of various immune cell in the TME using R package "GSVA" (25). The characteristic genes for each immune cell were obtained from other research(26). (27) was applied to value stromal and immune microenvironment in ltration. ESTIMATE score is de ned as tumor purity, which is the sum of stromal and immune score.

Gene Set Enrichment Analysis (GSEA)
To reveal the underlying function of potential prognostic biomarkers for ACC, we performed gene set enrichment analysis (GSEA). We divide ACCs into high gene expression group and low gene expression group. In our study, "DSigDB.ALL.gmt" obtained from Molecular Signatures Database (MSigDB) and "c2.cp.kegg.v6.0.symbols.gmt" were selected as the reference sets. Gene size ≥ 100, FDR < 25%, and |enrichment score (ES)|> 0.60 were considered as the cut-off criteria for signi cantly enriched KEGG pathways. Figure 1 provided the research process of this study Co-expression network construction to identify clinically signi cant modules TCGA_PK_A4HP and TCGA_OR_A5JB were excluded from the subsequent study ( Fig. 2(a)). Seventyseven samples with were included in the analysis (Fig. 2(b)). In this research, the power of β = 4 (scalefree R^2 = 0.92) was determined as the soft thresholding to construct a scale-free network (Fig. 3).

GO Biological Processes and KEGG Analysis
We conducted GO and KEGG pathway enrichment analysis to explore the underlying function of genes in prognosis-related modules. On the basis of GO analysis, genes in the blue module were enriched in four BPs, including pattern speci cation process, regulation of membrane potential, embryonic organ morphogenesis and embryonic skeletal system development ( Fig. 5(a)). As for the brown module, genes were enriched in 209 BPs (Supplementary Table S1). The top 10 terms were organelle ssion, nuclear division, chromosome segregation, nuclear chromosome segregation, mitotic nuclear division, DNA replication, DNA conformation change, sister chromatid segregation, mitotic sister chromatid segregation and regulation of chromosome segregation (Fig. 5(c)). While in KEGG pathway enrichment analysis, genes in blue module were enriched in two KEGG pathways, including neuroactive ligand-receptor interaction and axon guidance (Fig. 5(b)). Meanwhile, genes in the brown module were enriched in six pathways, including cell cycle, oocyte meiosis, fanconi anemia pathway, progesterone-mediated oocyte maturation, DNA replication and homologous recombination (Fig. 5(d)).

Identi cation of hub genes
First, 1770 genes that reached the cutoff criterion (|GS|≥0.2) for both OS time and OS were selected ( Fig. 6(a)). Then we got 16 genes from the blue module based on the screening criterion (|GS|≥0.2, |MM| ≥0.8) (Fig. 6(b)). In the same way, we got 118 genes from the brown module (Fig. 6(c)). In all, 134 hub genes in key modules were screened out. Moreover, we tried to nd hub genes associated with immunity. Therefore, two genes including BIRC5 (baculoviral IAP repeat containing 5) and PRLR (prolactin receptor), which were overlapping in WGCNA and IRGs gene lists, were eventually identi ed as immune-related prognosis biomarkers for further validation (Fig. 6(d)).

Hub Gene Validation
After screening out two biomarkers, we performed internal validation. At the beginning, according to the GEPIA, higher BIRC5 expression was signi cantly correlated with worse OS (Hazard ratio [HR] = 6.5, p = 6.5e-06, Fig. 7(a)) and DFS (HR = 3.3, p = 0.00054, Fig. 7(b)). Furthermore, we found that patients in high PRLR expression group had better OS (HR = 0.19, p = 7.2e-05, Fig. 7(e)) and DFS (HR = 0.29, p = 0.00039, Fig. 7(f)). Then we made a comparison of the mRNA expression of selected genes in ACC tumor tissue with normal tissue, the results revealed higher expression of BIRC5 (p < 0.05, Fig. 7(c)) in ACC samples compared with normal samples. It should be pointed out that although it was not signi cantly that PRLR expression in normal adrenal glands was higher than that in ACC tissues, there was still such a trend that the median expression of PRLR in normal adrenal glands was higher than that in ACC samples ( Fig. 7(g)). Furthermore, the results also showed that patients with high BIRC5 expression (F = 6.06, p = 0.000981; Fig. 7(d)) had higher tumor stage. On the contrary, low expression of PRLR (F = 3.53, p = 0.0191; Fig. 7(h)) was signi cantly correlated to higher tumor stage.

Relationship Between Gene Expression And TME
Varieties of immune cells in TME play an vital role in regulating anti-tumor response. Therefore, we conducted ssGSEA to determine the degree of alteration in the particular component of tumor in ltrating lymphocytes (TILs). Our study involved 28 immunophenotypes, including B cells, T cells, dendritic cells (DC) and other tumor-in ltrating lymphocytes. The results revealed that activated CD4 T cell and type 2 T helper cell were positively related with the expression levels of BIRC5 ( Fig. 12(a)). And the level of CD56dim natural killer cell and type 17 T helper cell was positively related with PRLR expression (Fig. 12(c)). Pearson's correlation analysis suggested that the abundance of anti-tumor immune cells and pro-tumor immune cells was positively correlated in local environment (BIRC5: R = 0.9008, p < 0.001, Fig. 12(b); PRLR: R = 0.9008, p < 0.001, Fig. 12(d)).
On the other hand, we also used DSigDB to explore the drug pathways enriched in BIRC5 high expression  (Fig. 13(c)).

Discussion
Adrenocortical carcinoma is prone to local invasion and metastasis, with poor prognosis and limited treatment(28). Currently, the preferred treatment for ACC is radical resection of the solid tumor (29).
However, almost half of the patients present disseminated metastasis at the rst diagnosis, and locoregional metastases or recurrence occur in approximately one-third of patients after surgery (30). For these patients, mitotane is the only approved treatment by the Food and Drug Administration and is generally poorly tolerated with limited e cacy (31). Recently, the development of immunotherapy has changed the treatment model of various cancers (32), and increasing researches have focused on immune-related molecules. However, some studies had shown that current immunotherapy did not show the expected effect on ACC. Fiorentiniv C et al pointed out that in numerous cancers including ACC, the preliminary clinical ndings about anti-PD-L1 drugs seemed to be unsatisfactory (33). ACC was even considered by some scholars to be a tumor with low immune response (34). Therefore, it is urgently to explore novel immune-related prognosis biomarkers for ACC.
In this study, we constructed a weight gene co-expression network by using TCGA-ACC data. Tow clinical features (OS time, OS) were considered as important prognostic indicators. We also tried to nd the key modules which were most related to the OS Time or OS of ACC. Among 23 modules, the blue and brown modules were the most associated with OS Time and OS. One hundred and thirty-four genes which reached the cutoff were considered as hub genes in key modules. To determine the immune-related genes among the hub genes, we used the ImmPort database. Finally, BIRC5 and PRLR were identi ed as novel immune-related prognosis biomarkers for ACC.
We used two databases to perform the validations of hub genes. GEPIA analysis con rmed that PRLR expression was lower in ACC than that in normal adrenal glands based on the TCGA-ACC data. In addition, Oncomine analysis showed that PRLR expression was higher in normal adrenal glands than in ACC samples. And BIRC5 showed extremely opposite results. Furthermore, whether on the basis of GEPIA analysis, GSE10927, GSE76019, or GSE76021 database, PRLR and BIRC5 were signi cantly related with the prognosis of ACC, demonstrating that PRLR and BIRC5 were powerful immune-related prognostic biomarkers in ACC.
Since the level of immune in ltration has a strong correlation with the prognosis of tumors(35), we used the ESTIMATE algorithm and ssGSEA to study the relationship between the hub gene expression and the level of immune cell in ltration. BIRC5 expression was signi cantly positively correlated with activated CD4 T cell and type 2 T helper cell. And PRLR expression was signi cantly positively related to CD56dim natural killer cell and type 17 T helper cell. Thus, BIRC5 and PRLR might be potential therapeutic targets for immunotherapy in ACC. The analysis of immune cell that associated with hub genes further demonstrated the possible action pathway of immune-related genes. BIRC5 and PRLR may in uence the prognosis of ACC by different immune pathways. According to DSigDB, we screened six potential drug targets for BIRC5 with GSEA analysis, which may become the targets of immunotherapy.
BIRC5 (also called survivin) is one of the inhibitor of apoptosis protein(36), which can inhibit apoptosisrelated signaling pathways and promote cell proliferation (37). It has been previously identi ed as one of the prognostic biomarkers for malignant tumor (38-41). Functionally, BIRC5 can stabilize the mitotic apparatus, ensure proper segregation of chromosomes and maintain microtubule integrity. It contributes to complete cell division stably and effectively (42)(43)(44)(45)(46), which is indispensable in the tumor cell growth. In addition to being a contributing factor of almost all cancer type, increasing data indicates a critical role of BIRC5 in regulating the proliferation of non-malignant, physiological adult T cells (40,47,48). This is in line with our study, in which the expression of activated CD4 T cell and type 2 T helper cell was signi cantly higher in high expression BIRC5 samples than that in low expression BIRC5 samples.
Although BIRC5 is considered as a vital cancer gene, only limited progress has been made in nding novel treatments targeting this pathway so far (49). Under these circumstances, scholars have studied the survivin inhibitor YM155, which is expected to shut down the transcription of the BIRC5 gene (50). But further study of YM155 has been less successful (51). Therefore, whether BIRC5 could be immune-related therapeutic target still requires more researches.
PRLR (prolactin receptor) is one of the cytokine receptor and its role in the genesis and progression of tumors is still controversial. Prolactin (PRL) is not only expressed on the surface of breast epithelial cells, but also on the surface of many other cells (52). Knockout of the PRLR locus suppresses the progression of breast tumors (53). In humans, PRLR expression higher in malignant epithelium than adjacent normal tissue (54). But various research on the relationship between PRLR and PRL and breast tumors remains inconclusive (55). On the other hand, the expression of prolactin in prostate cancer is correlated with high tumor grade and aggressive course(56, 57). However, the relationship between PRLR and adrenal tissue is uncertain. Glasow's research showed that PRL inhibits the proliferation of adrenal cells and PRLR is expressed in adrenocortical tumors but not in pheochromocytoma(58). Our results actually suggested that PRLR contributed to a favorable outcome in the prognosis of ACC, which is in accordance with Glasow's results. These differences may be due to differences in research methods, or differences in the biological characteristics of PRLR in different tissues.
There are several limitations of our study. First, there are no experimental data (in vivo and in vitro validation) to support these ndings at present. Secondly, the relationship between the two hub genes has not been further explored. Finally, the mechanism of the in uence of the two biomarkers on the prognosis of ACC needs to be further studied. Therefore, more molecular biology experiments are supposed to be carried out to study the functions of the hub genes.

Conclusions
In conclusion, by using various bioinformatics methods, we identi ed two immune-related prognostic biomarkers (BIRC5, PRLR) of ACC. Furthermore, the underlying function of genes related to prognosis were predicted. These biomarkers have great potential in translational medicine and may become immunotherapy targets for ACC in the future..

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. TableS1.xlsx