Identification ACTA2 and KDR as key proteins for prognosis of PD‐1/PD‐L1 blockade therapy in melanoma

Abstract Programmed cell death protein 1 (PD‐1) /programmed cell death ligand 1 (PD‐L1) blockade is an important therapeutic strategy for melanoma, despite its low clinical response. It is important to identify genes and pathways that may reflect the clinical outcomes of this therapy in patients. We analyzed clinical dataset GSE96619, which contains clinical information from five melanoma patients before and after anti‐PD‐1 therapy (five pairs of data). We identified 704 DEGs using these five pairs of data, and then the number of DEGs was narrowed down to 286 in patients who responded to treatment. Next, we performed KEGG pathway enrichment and constructed a DEG‐associated protein‐protein interaction network. Smooth muscle actin 2 (ACTA2) and tyrosine kinase growth factor receptor (KDR) were identified as the hub genes, which were significantly downregulated in the tumor tissue of the two patients who responded to treatment. To confirm our analysis, we demonstrated similar expression tendency to the clinical data for the two hub genes in a B16F10 subcutaneous xenograft model. This study demonstrates that ACTA2 and KDR are valuable responsive markers for PD‐1/PD‐L1 blockade therapy.

in patients with advanced or metastatic urothelial carcinoma in whom the disease has progressed after platinum-based chemotherapy. 15 Thus, it is important to explore biomarkers that can be used to assess the clinical outcomes of PD-1/PD-L1 blockade therapy.
Gene chip and differential gene profiles are efficient techniques that can probe the gene expression patterns in a sample at a time point. Gene expression levels and profiles obtained from microarrays have been widely used for identifying differentially expressed genes (DEGs). 16 The combination of gene expression microarrays and bioinformatics can be used to find potential biomarkers. 17 In this study, based on published clinical datasets from five melanoma patients treated with anti-PD-1, 18,19 we screened for DEGs.
Next, we analyzed the signaling pathways involved and constructed the protein-protein interaction (PPI) network for the DEGs. The expression of the key genes identified was measured in relevant animal models to validate these analyses. Our findings provide potential biomarkers for supporting the clinical use of PD-1/PD-L1 blockade therapy.

| Gene expression profile data
The Gene Expression Omnibus (GEO) datasets at the National Center for Biotechnology Information are an international public repository that contains microarray, next-generation sequencing and other formats of high-throughput functional genomic data. 20 The gene expression profile dataset, GSE96619, from GEO was used in the study. There were 10 biopsy specimens derived from five melanoma patients who responded (complete response/partial response, n = 2) or did not respond (progression, n = 3) to anti PD-1 therapy. Tumor samples were obtained from patients receiving anti-PD-1 therapy.
Samples were immediately fixed in formalin followed by paraffin embedding and processed for snap-freezing in liquid nitrogen. 19

| Identification of DEGs
We used MetaboAnalyst (https://www.metab oanal yst.ca/) to identify DEGs. The fold change of the gene expression level was calculated. Only genes with a P <.05 and fold change >1.5 or <0.67 were regarded as DEGs.

| Functional and pathway enrichment analysis of the DEGs
The significant GO terms and KEGG pathway enrichment analysis of the identified DEGs were performed by using DAVID with thresholds of significant functions and pathways of P <.05 and enrichment gene count of >5. 21,22 DAVID has been used for systematic and comprehensive analysis of massive lists of genes.

| Establishment of PPI network, modular analysis, and pathway identification
The STRING v10 online tool was used to construct and visualize the PPI network. 23 We mapped the DEGs into STRING, and interactions with the threshold of combined score of >0.4 were selected as significant. The network was captured and modified by Cytoscape (http://www.cytos cape.org/). 24 MCODE in Cytoscape was used for integrating the complex PPI network into unified conceptual frameworks and calculating the node degree (numbers of interconnections to filter hub genes). The hub genes were then selected with a cutoff degree of ≥10. 25

| GEPIA analysis of gene expression
GEPIA (http://gepia.cance r-pku.cn/) was used to validate the hub genes. GEPIA is an online tool based on the sequencing database for gene expression analysis between tumor and normal data from The Cancer Genome Atlas and the Genotype-Tissue Expression programs. We used GEPIA for the preliminary exploration of the differences in tyrosine kinase growth factor receptor (KDR) and smooth muscle actin 2 (ACTA2) expression level and the survival rate between skin cutaneous melanoma (SKCM) and normal samples. 26

| Establishment of mouse xenograft model
Eight-week-old male C57BL/6 mice were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing Vital River Laboratory Animal Technology Co., Ltd., Beijing, China). On day 0, B16F10 mouse melanoma cells were harvested in saline, and 1.5 × 10 6 cells in 0.2 mL saline were injected subcutaneously into the right flanks of each mouse. On day 1, treatment was initiated after the mice were assigned randomly to the control and experimental groups. In the control group, IgG (BE0093, Bio X Cell, West Lebanon, NH) was dissolved in saline for intraperitoneal treatment every 3 days. In the PD-L1 mAb treatment group, 10 mg/kg anti-mouse PD-L1 (BE0101, Bio X Cell, West Lebanon, NH) antibody was dissolved in saline for intraperitoneal treatment every 3 days. For the CTX group, 60 mg/kg CTX was dissolved in saline for intraperitoneal treatment every 7 days. The Tumor volume was calculated as Tumor volume = π × L × W 2 /6, in which L is the maximum length of the tumor and W is the maximum width of the tumor. Mice were sacrificed by CO2 asphyxiation. When the mice were sacrificed, the tumors were stripped and weighed. The

| RT-qPCR
Real-time quantitative PCR (RT-qPCR) with primers was conducted as previously described. 32 Total RNA from mouse tumor tissues

| IHC staining
Samples were processed for IHC by routine techniques, as previously described. 35,36 Xylene was used to dewax the paraffinembedded sections. The deparaffinized tissue sections were incubated with 3% H 2 O 2 for 10 minutes at 37°C to quench the activity of endogenous peroxidase. Proteinase K was used to digest the sections for antigen retrieval. Slides were incubated over-

| Statistical analysis
Data were expressed as the mean ± SD of at least three independent experiments. Statistical analysis was carried out with Prism (version 7.0, Graph Pad, San Diego, CA). Statistical significance was inferred at P <.05.

| Sample collection from the dataset
Clinical data from melanoma patients who received PD-1/PD-L1 blockade therapy was analyzed. The original microarray dataset GSE96619 from the GEO database contained five pairs of melanoma tissue samples (before and after receiving anti-PD-L1 therapy with atezolizumab) from five patients. Of these five patients, two responded (complete response/partial response, Pt1 and Pt2) to the treatment, and the other three did not respond (progression, Pt3, Pt4, and Pt5). Dataset GSE96619 was also analyzed in other publications with different aims. 18,19

| Identification of DEGs and analysis of gene expression
We identified 704 DEGs involved in antitumor immunity by comparing expression in the responding group and no responding group using a P-value of .05 (unpaired t test) as the criterion. Next, we set the criteria as P <.05 and fold change >1.5 or <.67, and 286 DEGs (232 upregulated and 54 downregulated genes) were identified in the responding group after the data from the treatment biopsy specimens were compared with their respective baselines (Table 1).
MetaboAnalyst software was used to integrate the 286 DEGs in an expression heat map of the significant DEGs' differential distribution ( Figure 1). The full figure is provided in Figure S1.

| Gene ontology analysis of DEGs
To understand the biological functions of our screened DEGs, we performed gene ontology (GO) analysis of the DEGs in DAVID (https://david.ncifc rf.gov/) 37 with the criteria of P <.05 and count of ≥5. The DEGs were summarized into the following GO categories: "biological process" (BP; 145 GO terms), "molecular function" (MF; 13 GO terms), and "cellular component" (CC; 1 GO terms). Detailed information is provided in Table S1.
In the BP group, the upregulated genes were mainly enriched in the "cellular response to chemical stimulus", "regulation of cell communication", and "regulation of signaling" GO terms. The downregulated genes were mainly concentrated in the "ncRNA metabolic process", "carboxylic acid metabolic process", and "oxoacid metabolic process" terms. In the MF group, the upregulated genes were mainly enriched in the "calcium ion binding", "substrate-specific channel activity", and "channel activity" terms. The downregulated genes were mainly enriched in the "structural molecule activity" term. In the CC group, the upregulated genes were mainly enriched in the "extracellular region", "extracellular region part", and "membrane-bounded vesicle" terms. The downregulated genes were mainly enriched in the "membrane-bounded vesicle", "extracellular region part", and "extracellular exosomes" terms. The 75 most significant GO terms are shown in Figure 2A.
The "extracellular region" (70 involved genes), "extracellular region part" (53 involved genes), and "membrane-bounded vesicle" (51 involved genes) GO terms in the CC group contained the majority of enriched DEGs. The 50 most significantly enriched GO terms are shown in Figure 2B, according to P-value. The top five GO terms were "blood vessel development", "angiogenesis", "vasculature development", "anatomical structure formation involved in morphogenesis", and "blood vessel morphogenesis". A total of 41 GO terms were located in the BP group, whereas only three GO terms were located in the MF group and six GO terms in the CC group.

| KEGG pathway enrichment analysis of DEGs
We performed KEGG pathway enrichment analysis to elucidate the important pathways of DEGs. Enriched pathways of DEGs and detailed information are listed in Table 2. The KEGG pathway enrichment analysis of upregulated DEGs indicated that the important pathways were the "Ras signaling pathway" (eight genes), "TNF signaling pathway" (five genes), "adipocytokine signaling pathway" (four genes), and "prolactin signaling pathway" (four genes). The downregulated DEGs were mainly enriched in the "phenylalanine metabolism pathway" (two genes).

| PPI network construction and hub gene identification
The data were imported into the STRING online database (http:// strin g-db.org), and proteins were linked by colored lines to indicate the different types of interaction evidence. 23 The potential interactions of the DEGs were obtained by mapping the upregulated and downregulated DEGs in the STRING database. After removing the isolated nodes, a complex PPI network comprising 170 edges and 249 nodes was established ( Figure 3). Each protein in the network is considered as node and the edges represent the predicted functional associations. 23 The degree of a node represents the number of interactions between two nodes. Proteins closely associated with others in the network were identified with a degree of ≥10, and the hub genes, including, were identified KDR and ACTA2.

| Modular analysis of the PPI network
Using Cytoscape Molecular Complex Detection (MCODE), the most significant module from the PPI network complex was selected, and the genes involved in the modules were analyzed ( Figure 4). 25 Enrichment analysis of the module showed that the genes were mainly associated with the "adipocytokine signaling pathway", "TNF signaling pathway", "Ras signaling pathway", "prolactin signaling pathway", "PI3K-Akt signaling pathway", "focal adhesion", "Rap1 signaling pathway", and "ECM-receptor interaction".

| Verification and survival analysis of hub genes
The expression level of these two genes were analyzed from dataset GSE96619. In general, the expression levels of ACTA2 and KDR in no responding group was clearly higher than that in responding group TA B L E 1 286 differentially expressed genes (DEGs) were identified from GSE96619 dataset, including 232 up-regulated genes and 54 down-regulated genes in responding group on-treatment (OnTx) biopsy specimens compared to their respective baselines (The up-regulated genes were listed from the largest to the smallest of fold changes and down-regulated genes were listed from the smallest to largest)

| In vivo anti-tumor study of PD-L1 monoclonal antibody treatment
The anti-tumor effect of PD-L1 monoclonal antibody (mAb) was evalu- The responding and no responding group of PD-L1 mAb treatment resulted in a 72.4% and 0.3% decrease in tumor weight compared with the control tumor after 18 days of treatment respectively ( Figure 6, Table 3).

| Validation of the expression of hub genes in the B16F10 subcutaneous xenograft model
To examine the quality of our hub gene exploration, transcriptional and immunohistochemical (IHC) analyses of the expression of hub genes in the responding and no responding group of B16F10 subcutaneous xenograft model were compared (Figure 7). RT-qPCR showed that the expression of the two hub genes was higher in the no responding group ( Figure 7A). We verified the overexpression of the two hub genes by IHC staining. The expression of the two hub genes were significantly elevated in no responding individuals compared with responding group (Figure 7B,C). Expression level of PD-L1 and CD4 + ; CD8 + positive cells were also detected by IHC staining.
Both the expression level of PD-L1 and number of CD4 + CD8 + positive cells were higher in the responding group ( Figure 7D,E).

F I G U R E 1
The heat map and cluster diagram of differential expression profiles of DEGs from dataset GSE96619 was developed. Each column represents an independent sample, and each row represents a gene. Blue represents the downregulated genes expression and red indicates the upregulated genes in each cell. The full figure with details is provided in Figure S1 F I G U R E 2 GO analysis and significantly enriched GO terms of DEGs. A, The 75 most significant GO terms for DEGs classified into three groups (ie, "molecular function", "biological process", and "cellular component"  to avelumab. 12 Lastly, the response to cemiplimab was observed no more than 50% in patients with advanced cutaneous squamous-cell carcinoma. 38 To understand treatment better, predictive biomarkers for evaluating the effectiveness of PD-1/PD-L1 blockade drugs should be identified.
Clinical research data is precious for providing a great deal of other important information that may have been irrelevant to the original analysis. Remining of these datasets of clinical research may help to identify the potential predictive biomarker. 39 The GSE96619 dataset used in this study contains a balanced responding and no responding cases to anti PD-1 therapy. 18 Given that the statistical analysis is properly performed, the data can be used to obtain highquality results.
Our analysis identified the 50 most significantly enriched GO terms from dataset GSE96619, and the top five were "blood vessel development", "angiogenesis", "vasculature development", "anatomical structure formation involved in morphogenesis", and "blood vessel morphogenesis". The terms are all related to ergistically. 45 In addition, PD-L1 expression is also significantly cor- In conclusion, we performed a bioinformatics analysis using clinical dataset GSE96619. We screened 704 DEGs that may be relevant to PD-1/PD-L1 blockade therapy. By analyzing the GO and KEGG pathways, we found that DEGs were mainly enriched in "angiogenesis" terms, which provides a theoretical basis for PD-1/PD-L1 blockade therapy. We constructed a PPI network of DEGs and identified two hub genes (KDR and ACTA2) that could be potential biomarkers for predicting prognosis. We confirmed the identification by observing similar expression of the two hub genes in a B16F10 subcutaneous xenograft model and clinical samples. Our study revealed that ACTA2 and KDR could be used as responsive markers for PD-1/PD-L1 blockade therapy in melanoma.

| ACKNOWLEDG EMENT
We appreciate the guidance on this project provided by Dr Yi Zhu and Dr Xu Zhang at Tianjin Medical University, Dr Qinghua Cui and Jiangcheng Shi at Peking University Health Science Center. We would like to thank Angel Garcia-Diaz and colleagues for uploading and sharing their dataset GSE96619.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

AUTH O R S' CO NTR I B UTI O N S
YW carried out the study and drafted the manuscript, ZL collected and analysed data, ZZ did animal experiment, XC designed the . Data are reported as mean ± SD from at least three independent experiments, which were performed in triplicate. (*P <.05, **P <.01, ***P <.001) project and revised the manuscript. All authors read and approved the final manuscript.

E TH I C S A PPROVA L A N D CO N S E NT TO PA RTI CI PATE
The study was reviewed and approved by the Experimental Animal

Management and Welfare Committee at the Institute of Materia
Medica, Peking Union Medical College.

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasets analysed during the current study are available from the GEO dataset GSE96619 repository, https://www.ncbi.nlm.nih.