Study of LY9 as a potential biomarker for prognosis and prediction of immunotherapy efficacy in lung adenocarcinoma

Differential expression and prognostic value analysis of LY9 in pan-cancer

First, we retrieved the transcriptome and clinical data of 33 types of cancer from the UCSC Xena database (https://xena.ucsc.edu/) and obtained LY9 expression levels for all tumor samples using the “limma” R package. Considering that too few normal samples might lead to errors in the results, we deleted tumors with less than five normal samples, then used the “ggpubr” R package to draw a box plot of the LY9 expression level in the remaining tumors (n = 18, *p < 0.05, **p < 0.01, ***p < 0.001). Then we conducted Cox analysis through Xenashiny website (https://shiny.zhoulab.ac.cn/UCSCXenaShiny/) to measure the prognostic value (overall survival, OS) of LY9 in 33 cancers from the TCGA database. We used the Kaplan–Meier plotter (https://kmplot.com/analysis/) (Györffy et al., 2010; Győrffy, 2021) and the Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) (Tang et al., 2017) databases to explore the prognostic value of LY9 in the abovementioned tumors (n = 18). The purpose of using these two databases for prognostic analysis is to serve as a cross-validation of results. Sources of the Kaplan–Meier plotter include GEO, EGA, and TCGA databases. When we performed survival analysis in the GEPIA database, the statistical difference between the curves was measured by the log-rank test. Portions of this text were previously published as part of a preprint (Deng et al., 2022).

Gene Set Enrichment Analysis and correlation of LY9 expression in LUAD with clinical features from TCGA

To explore the functions and pathways that enriched for LY9, we used the Gene Set Enrichment Analysis (GSEA) software (v. 4.1.0) to perform functional enrichment analysis based on the LY9 expression in LUAD. First, we classified patients into high- and low-expression groups based on the median value of the LY9 expression levels in LUAD patients. Then, we explored the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway in the two LY9 expression groups in LUAD patients. If the conditions of —NES—>1, NOM p-value <0.05, and FDR q-value <0.25 were satisfied, we considered that the gene set under the pathway was meaningful. Before analyzing the LUAD dataset from the TCGA database, we first performed filtering according to fdrFilter <0.05. We also used the “ggpubr” R package to draw boxplots to demonstrate the association of LY9 expression in LUAD patients in the TCGA database with other clinical features, including age, gender, clinical stage, and TMN stage. In this study, we removed the samples with incomplete clinical information, and the clinical information data are summarized in Table 1. In addition, to verify the differential expression of LY9 in LUAD, we downloaded the corresponding dataset from the GEO dataset (GSE68465) for analysis. Moreover, we explored the protein-protein interaction (PPI) network associated with LY9 using the Genecards database (https://www.genecards.org/). Finally, we explored the prognostic value of LY9 expression levels in LUAD patients by univariate and multifactorial independent prognostic analyses. M stage was not included because the sample with an unknown M stage accounted for a large proportion of the total sample (Table 1, 26.6%) and would have introduced a large error in the results.

Construction of the lncRNA-miRNA-LY9 ceRNA network

The miRNA expression data was downloaded from TCGA data. Then we retrieved miRNAs predicted to target LY9 from the TargetScan 7.2 database (http://www.targetscan.org/vert_72/) for correlation and differential expression analyses and filtered the predicted upstream miRNAs of LY9 according to correlation >0.2 and p < 0.05. The statistical method was Spearman correlation analysis. Also, we constructed correlation networks for miRNAs with correlations >0.1 and performed correlation, differential, and prognostic survival analyses for hsa-miR-141-3p that met the required conditions and the results were statistically significant. Thereafter, we downloaded lncRNAs predicted to target hsa-miR-141-3p from the DIANA-tools database LncBase Experimental v.2 (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=lncExperimental/index) and performed correlation and differential expression analyses. We found that RNF216P1 and LINC00943 correlations were >0.2 and p < 0.05. Finally, we constructed the LY9 lncRNA-miRNA network and performed correlation, differential expression, and survival analyses of RNF216P1 and LINC00943 with LY9 expression.

Correlation analyses of LY9 expression in LUAD with immune cells infiltration and immune checkpoints

To explore whether LY9 could be an immunotherapy target for LUAD, we analyzed the correlation between LY9 expression levels and the infiltration of various immune cells in LUAD from the TIMER 2.0 database (http://timer.cistrome.org/), including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. Next, we calculated the expression of immune cells for each LUAD patient in the TCGA dataset using the CIBERSOFT algorithm and calculated the correlation between various immune cells. Results were plotted as histograms and correlation heat maps using the “corrplot” R package. Then, we used differential expression analysis and the Spearman correlation test to explore which immune cell infiltration was associated with LY9 expression in LUAD patients. Finally, we used the intersection of immune cells obtained by both methods to obtain the final results. Moreover, we used these two databases (TIMER 2.0 and GEPIA 2.0) to analyze the correlation between the LY9 expression in LUAD and three main immune checkpoints (CD274, PDCD1, and CTLA4). Further, we calculated the correlation of LY9 with various immune cell marker genes by Spearman analysis.

Validation of GEO dataset, qRT-PCR and IHC experiments

GEO database LUAD RNA-seq dataset GSE68465 for validation of differential expression of LY9. We collected 26 pairs of LUAD tissues and their paraneoplastic normal tissues from the First Affiliated Hospital of Guangxi Medical University from January 2015 to September 2021 to perform qRT-PCR experiments to explore the mRNA expression levels and locations of LY9, CD4, and CD8, and calculated their correlation. Written informed consent was obtained from each subject and our study was approved by the Medical Ethics Committee of the First Affiliated Hospital of Guangxi Medical University (NO. 2022-KY-E-063). All experiments in this study were in accordance with the Declaration of Helsinki. First, we extracted total RNA and synthesized cDNA using random hmers and SuperScript III reverse transcriptase, repeating this step three times, and taking the mean value for subsequent qRT-PCR. The final results were estimated by the 2^−ΔΔCt method (Cycle threshold value, Ct). A p < 0.05 was considered statistically significant. The following primers were used: (5′–3′): forward sequence of LY9: TCTCTCCTCTTCCTGCTCATG and reverse sequence of LY9: ACGTTCTCAATCTCTGTGTCTAC; forward sequence of CD8: GGACTTCGCCTGTGATATCTAC and reverse sequence of CD8: TTGTCTCCCGATTTGACCAC; forward sequence of CD4: GCCCTTGAAGCGAAAACAG and reverse sequence of CD4: CTCCTTGTTCTCCAGTTTCAAAC; GAPDH was used as an internal control: forward sequence (5′–3′): ACATCGCTCAGACACCATG and reverse sequence (5′–3′): TGTAGTTGAGGTCAATGAAGGG.

When performing IHC experiments, anti-LY9 was rabbit antibody purchased from Solarbio-Beijing, and the goat anti-rabbit secondary antibody was purchased from Beyotime. All experimental reagents were used in strict accordance with the instructions, and the experimental procedures were not special.

Statistical analyses

The version of R software used in our study was 4.1.1 (R Core Team, 2021). We used the “limma” R package to perform Wilcoxon tests to analyze variances. In the survival analysis using the Kaplan–Meier database, the log-rank test was used. Spearman’s correlations were performed to analyze correlations. A p < 0.05 was considered statistically significant.

The mRNA expression levels of LY9 in pan-cancer

The results from the Xenashiny website were presented in a forest plot, from which we learned that LY9 expression could be a risk factor in Uveal Melanoma (UVM). Meanwhile, LY9 may play a protective role in skin cutaneous melanoma (SKCM), sarcoma (SARC), lung adenocarcinoma (LUAD), head and neck squamous cell carcinoma (HNSC), cervical squamous cell carcinoma, and endocervical adenocarcinoma (CESC), breast invasive carcinoma (BRCA), and adrenocortical carcinoma (ACC) (Fig. 2A). To determine whether LY9 expression levels differed between cancerous and paraneoplastic tissues in various tumors, we retrieved transcriptomic and clinical data of 33 tumors from the UCSC database for further analysis. Results showed that tumors with higher LY9 expression in normal tissues than in cancerous tissues were bladder urothelial carcinoma (BLCA, p < 0.05), colon adenocarcinoma (COAD, p < 0.001), lung squamous cell carcinoma (LUSC, p < 0.05), rectum adenocarcinoma (READ, p < 0.001), and thyroid carcinoma (THCA, p < 0.001). However, the tumors with higher LY9 expression, compared to normal tissues, were glioblastoma multiforme (GBM, p < 0.05), kidney renal clear cell carcinoma (KIRC, p < 0.001), kidney renal papillary cell carcinoma (KIRP, p < 0.05), and lung adenocarcinoma (LUAD, p < 0.05) (Fig. 2B). The above results indicated that the mRNA expression level of LY9 was up-regulated in a variety of solid tumors, including LUAD, and acted as a protective factor.

Prognostic value of LY9 in pan-cancer

We explored the association of LY9 with the overall survival (OS) in patients with various tumors using the Kaplan–Meier plotter (Figs. 3A–3F) and GEPIA (Figs. 3G–3L) databases. Results showed six tumors significantly associated in both databases, including BRCA, HNSC, LUAD, CESC, LIHC, and OV. Among these tumors, patients with high LY9 expression had longer OS. These results suggested that LY9 was a protective factor for these six cancers.

GSEA and correlation analysis of clinical characteristics in LUAD

Further, the Wilcoxon Tests for LY9 expression in samples from the LUAD dataset (TCGA database) indicated that LY9 was indeed expressed at a higher level in cancerous tissues (Fig. 4A, p = 0.013; B, p = 0.0022). Also, the survival analysis confirmed that patients with high LY9 expression had longer OS (Fig. 4C, p = 0.005). The results of correlation analysis with other clinical characteristics showed that, in LUAD patients, the LY9 expression tended to increase with age (Fig. 4D, p = 0.003). Also, the LY9 expression level was higher in females than in males (Fig. 4E, p = 0.02). Additionally, the overall trend of clinical and TMN staging showed a decreasing LY9 expression trend in early to advanced patients (Figs. 4F–4I), consistent with the result that patients with higher LY9 expression in tumors have better prognoses. The univariate independent prognostic analysis results showed that the clinical stage, T-stage, and N-stage were high-risk factors for LUAD patients, while the LY9 expression level was a low-risk factor (Fig. 4J). Similarly, the results of the multifactorial independent prognostic analysis suggested that the expression level of LY9 in LUAD is a protective factor (Fig. 4K). Moreover, we classified the LUAD transcriptome, retrieved from the TCGA database, into high- and low-expression groups, based on the expression of LY9, then performed the GSEA. Results indicated that LY9-related gene sets were involved in many immune-related and non-small cell lung cancer pathways in the high-expression group. We selected the top 11 pathways to draw enrichment maps (Fig. 4L). In the low-expression group, LY9-related gene sets were mainly involved in the biosynthesis of unsaturated fatty acids. Together, these results implied that the prognosis of LUAD patients can be assessed by LY9 expression levels. LY9 was primarily involved in immune-related pathways, which inspired us to further explore the value of LY9 immunotherapy.

Construction of the lncRNA-miRNA-LY9 ceRNA network

To explore the potential lncRNA-miRNA-mRNA ceRNA action network of LY9, we conducted an analysis using online databases. First, we retrieved miRNAs predicted to target LY9 from the TargetScan 7.2 database for correlation and differential expression analyses and filtered the predicted upstream miRNAs of LY9 according to correlation >0.1 and p < 0.05. Results showed that seven miRNAs, including hsa-miR-141-3p, hsa-miR-4746-5p, hsa-miR-151a-5p, hsa-miR-17-3p, hsa-miR-301b-5p, hsa-miR-4786-3p, and hsa-miR-629-5p were predicted to target LY9 (Table 2). Also, we constructed correlation networks for these seven miRNAs (Fig. 5L). Considering that the miRNA with the strongest correlation with LY9 was hsa-miR-141-3p (R = −0.2, p = 9.8e −06), we further performed differential and prognostic survival analyses for hsa-miR-141-3p. Results showed that the expression of LY9 was expressed at higher levels in LUAD tissues than in normal ones (Fig. 5B, p = 2.22e−16), and patients with high hsa-miR-141-3p expression had longer OS (Fig. 5C, p = 0.012). Thereafter, we downloaded lncRNAs predicted to bind hsa-miR-141-3p from the DIANA-tools database LncBase Experimental v.2 and performed correlation and differential expression analyses. We found that RNF216P1 and LINC00943 had a correlation >0.2 and p < 0.05 (Table 3). Then, we performed differential expression and survival analyses of RNF216P1 and LINC00943 with LY9 expression. Results showed that the correlation between LY9 and RNF216P1 was not strong and not statistically significant (Fig. 5I, r = 0.055, p = 0.23). Meanwhile, a stronger and more significant positive correlation was found between LINC00943 and LY9 (Fig. 5E, r = 0.44, p < 2.2e −16). LINC00943 was also expressed at higher levels in LUAD tissues (Fig. 5F, p = 2.8e−07). This higher expression was associated with better OS (Fig. 5G, p = 0.018). Therefore, we hypothesized the existence of a ceRNA regulatory network of LINC00943 - hsa-miR-141-3p - LY9 (Fig. 5L). Finally, we constructed the PPI network of LY9 using the GeneCards database (Fig. 5M). The results indicated that the proteins with co-expression association with LY9 were INPP5D, LAT, INPPL1, PTPN11, SYK, SH2D1A, ZAP70, and SH2D1B. However, these findings need more experiments to validate.

Correlation of LY9 expression levels in LUAD with immune cell infiltration and immune checkpoints

Through differential expression analysis and exploration of prognostic value in the above databases, we found that LY9 was differentially expressed in cancerous and paraneoplastic tissues and significantly correlated with prognostic value in LUAD patients. Therefore, we selected the TIMER 2.0 database to explore the correlation between the LY9 expression level in LUAD, and the infiltration of various immune cells in the TME. Results suggested a significant positive correlation with the level of B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells (Fig. 6E). We also used the CIBERSORT algorithm to calculate the relative expression levels of 23 immune cells in each LUAD patient and to calculate the correlation between these immune cells. Then, we used a differential expression approach (i.e., we divided the patients into high and low groups according to the LY9 expression level), then calculated whether there was a difference in immune cell expression between high and low LY9 expression groups. We used a violin plot to show the results (Fig. 6A). Then, we analyzed the correlation between LY9 expression and these 23 immune cells using the Spearman correlation test. Results are shown in the correlation coefficient bar graph (Fig. 6B). These two methods showed that LY9 expression in LUAD was associated with 13 immune cells (Fig. 6C). They were Plasma cells, T cells CD8, T cells CD4 memory resting, T cells CD4 memory activated, T cells regulatory (Tregs), NK cells activated, Monocytes, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells activated, Mast cells activated and Neutrophils. We also found that LY9 was significantly positively correlated with the B cell markers CD19 and CD79A, the CD8+T cell marker CD8, the CD4+T cell marker CD4, and the neutrophil marker CCR7 (Fig. 6D). This result indicated that LY9 might be a LUAD biomarker and can be used as a target for immunotherapy. Similarly, the TIMER database results showed that the LY9 expression in LUAD was strongly and positively correlated with the expression levels of CD8 and CD4 T cells (Fig. 6E). Additionally, the GEPIA (Figs. 6F–6H) and TIMER (Figs. 6I–6K) databases results suggested that the LY9 expression in LUAD was significantly and positively correlated with the immune checkpoints CD274, PDCD1, and CTLA4. These results demonstrated that LY9 was strongly positively associated with CD8+ and CD4+ T cell infiltration in the tumor microenvironment of LUAD, and was strongly positively associated with the expression of immune checkpoint genes. Therefore, we hypothesized that the expression status of LY9 might affect the therapeutic effect of immune checkpoint inhibitors in LUAD patients.

Validation using GEO dataset, qRT-PCR and IHC experiments

After averaging the Ct value obtained from the three amplifications, the relative quantification was performed by the 2^−ΔΔCt method. The results of PCR experiments showed that the expression of LY9 in LUAD tissues was significantly and positively correlated with CD8 and CD4 expressions (Fig. 7A, CD8: R² = 0.8541, p < 0.001; Fig. 7B, CD4: R² = 0.4325, p = 0.003). Results showed that the expression of LY9 in LUAD was closely and positively associated with CD8 and CD4 expressions. Also, differential expression analysis of the LUAD RNA-seq dataset in the GEO database demonstrated high expression of LY9 in tumor tissues (Fig. 7E, P < 0.001), and therefore, we did not further validate it by experiment. Unfortunately, 26 pairs of LUAD tumor samples were negative for IHC experiments (Fig. 7C and 7D, only two tumor samples results are shown), which means the expression of LY9 at the protein level was relatively low. This may be one of the reasons for the low number of studies on LY9 in solid tumors in recent years. Therefore, the conclusions of this study were based only on the mRNA level of LY9.