Downregulation of USP18 reduces tumor-infiltrating activated dendritic cells in extranodal diffuse large B cell lymphoma patients

Nine prognostic immune genes were identified in EN DLBCL

The analysis process is shown in Figure 1. To determine the key biomarker related to EN DLBCL, we analyzed RNA-seq profiles and clinical data from 46 DLBCL patients consisting of 25 LN DLBCL and 21 EN DLBCL from TCGA database. The baseline information of the samples is presented in Table 1. As shown in Figure 2A, 2B, DEGs consisting of 583 up- and 585 down-regulated genes between these two groups were illustrated by heatmap and volcano plot. Then, GO and KEGG enrichment analyses were performed to uncover the potential mechanism distinguishing between the development of LN and EN DLBCL. As shown in Figure 2C, 2D, immune-related pathways such as “cytokine-cytokine receptor interaction” and “JAK-STAT signaling pathway” were included in the top ten terms, indicating that immune-related mechanisms were involved in the developmental difference between LN and EN DLBCL.

The DEIGs were obtained by intersecting DEGs and immune-related genes. Ninety-three DEIGs, including 53 up- and 40 down-regulated DEIGs, are displayed in the heatmap in Figure 3A. Then, the DEIGs and prognostic characteristics were submitted to univariate Cox regression analysis to identify prognostic immune genes. As shown in Figure 3B, nine prognostic immune genes, including IFNA21, KIR2DL1, MUCL1, SFTPA2, MX1, USP18, CCL1, IGLV1-36, and GLP1R, were identified.

These nine genes were then integrated into the multivariate regression analysis to build a prognostic predictive model. To avoid overfitting of the predictive model, Lasso regression was performed. The AUC of the ROC curve was 0.892, indicating that all nine genes were essential for the model (Figure 3C). To further assess model fit, we also performed a Schoenfeld residuals test. As shown in Supplementary Figure 1A–1C, the slope of scaled residuals on time was zero, so the proportional hazards assumption in the Cox model conformed to the null hypothesis, indicating high accuracy of the model. According to the model, the risk score of each sample was calculated, and samples were divided into high- and low-risk groups with a median value of 0.786. As shown in Figure 3D, the Kaplan-Meier survival curve showed that patients in the high-risk group had a lower survival rate than those in the low-risk group (P < 0.001), further revealing the good effectiveness of the predictive model.

Then, we generated a risk curve and scatterplot to show the risk score and survival status of each individual with DLBCL. As shown in Figure 3E, patients in the high-risk group showed higher mortality than those in the low-risk group, which also indicated the high efficacy of the model. The expression of prognostic immune genes screened by Lasso regression is displayed by a heatmap in Figure 3F.

USP18 was the key immune gene in EN DLBCL

To further identify the key immune genes, the coexpression of DETFs and nine prognostic immune genes was performed. Five up- and two down-regulated DETFs were identified by intersecting cancer-associated TFs and DEGs, as shown in the heatmap and volcano plot (Figure 4A, 4B). The correlation analysis identified 7 regulatory pairs between DETFs and prognostic immune genes, as shown in Table 2. As shown in Figure 4C–4E, only the expression of CCL1, IFNA21 and USP18 from nine prognostic immune genes was significantly different between EN and LN DLBCL by the Wilcoxon test (P < 0.05). To identify the key immune genes related to EN DLBCL, DETF-related and extranodal involvement (ENI)-related immune genes were intersected. As shown in Figure 4F, two immune genes were found in both groups. Combined with the results of the correlation between DETFs and prognostic immune genes shown in Table 2, the regulatory pair of LIM homeobox 2 (LHX2) and ubiquitin-specific peptidase 18 (USP18) was most significant (R = 0.497, P < 0.001, positive). Consequently, USP18, which is downregulated in EN DLBCL compared with LN DLBCL, was recognized as the key immune gene.

To demonstrate the regulatory mechanism between LHX2 and USP18, chromatin immunoprecipitation followed by high-throughput DNA sequencing (ChIP-Seq) data from the Cistrome database was evaluated. As shown in Supplementary Figure 2, in vitro ChIP-Seq data confirmed the transcriptional regulation patterns between LHX2 and USP18 in multiple cell lines.

USP18 was positively associated with the MAPK pathway in EN DLBCL

To discover the pathway downstream of USP18, GSVA was conducted, and 27 KEGG signaling pathways between EN and LN DLBCL were identified. The correlations between USP18 and these 27 KEGG pathways were constructed by Pearson correlation analysis, as shown in Figure 5A. To determine the critical signaling pathway, GSEA was also conducted. Three key KEGG pathways, i.e., the arrhythmogenic right ventricular cardiomyopathy (ARVC) pathway, dilated cardiomyopathy pathway, and MAPK pathway, overlapped between GSEA and GSVA (Figure 5B, 5C). Considering the relevance to the disease, we focused on the MAPK pathway in the following analysis. The GSEA of the MAPK pathway is shown in Figure 5D. The correlation between USP18 and the MAPK pathway was fitted by linear regression. As shown in Figure 5E, USP18 was positively correlated with the MAPK pathway (R = 0.294, P < 0.05, positive), suggesting that USP18 might modulate the MAPK pathway in the development of EN DLBCL.

Downregulation of USP18 in EN DLBCL was correlated with the immune gene set of aDCs

To determine the immune responses involved in USP18, ssGSEA was applied. Fifteen immune gene sets were identified in DLBCL patients from 29 immune gene sets that were overexpressed in the tumor microenvironment [11]. Pearson correlation analysis between USP18 and immune gene sets in DLBCL was constructed, as shown by the heatmap in Figure 6A. The top three immune gene sets correlated with USP18 were aDCs (R = 0.694, P < 0.001, positive), type I IFN response (R = 0.673, P < 0.001, positive) and regulatory T cells (Tregs) (R = 0.551, P < 0.001, positive), as shown in Figure 6B–6D. Of these, the relationship between USP18 and aDCs was most significant, indicating that USP18 might affect aDCs in EN DLBCL.

The online database further validated the association between key biomarkers in our analysis

To minimize the bias of the results above, a multidimensional validation was performed. The expression of LHX2 and USP18 and key genes of potential pathways in primary DLBCL, normal nodal tissue, and various cell lines, and their association with prognosis, are summarized in Supplementary Table 1.

First, LHX2 (median rank 695, P < 0.001), IL2RA (median rank 669, P < 0.001), IL21R (median rank 868, P < 0.001) and CHST7 (median rank 782, P < 0.001) were highly expressed in primary DLBCL compared to normal tissue, while IL5RA (median rank 3,659, P = 0.128) showed no difference in any of the four comparisons (Supplementary Figure 3). The GEPIA results showed that the mRNA expression levels of USP18, IL2RA and IL21R were higher in tumor samples than in normal samples (Supplementary Figure 6). At the cellular level, USP18, TLR7, IL21R, GCNT1 and CHST7 were expressed in various tumor cell lines, while the expression of LHX2, IL2RA and IL5RA was low in CCLE (Supplementary Figure 4). The results from The Protein Atlas showed the protein expression of USP18, IL2RA, TLR7, IL21R and GCNT1 in normal lymph node tissue (Supplementary Figure 10).

In addition, an analysis of the genomic and clinical profiles with cBioPortal suggested that LHX2, USP18 and key genes in downstream pathways were prone to mutations, which were associated with poor prognosis (Supplementary Figure 5A–5D). The results also showed that USP18 was coexpressed with LHX2 (R = 0.61, P < 0.001), IL5RA (R = 0.45, P < 0.001), IL21R (R = 0.56, P < 0.001), GCNT1 (R = 0.37, P = 0.024) and CHST7 (R = 0.43, P < 0.001) (Supplementary Figure 5E–5I). Moreover, analysis in the other databases also presented a negative association of key genes and prognosis (Supplementary Figures 7–9). Supplementary Figure 11 shows the PPI network of LHX2, USP18, IL2RA, IL5RA, IL21R, TLR7, GCNT1 and CHST7 generated in String.

USP18 expression and the number of aDCs were low in EN DLBCL tumor tissues

To further verify the role of USP18, the expression of USP18 in tumor biopsies of patients with LN and EN DLBCL was detected by IHC staining. The clinical information of DLBCL patients is shown in Table 3. As shown in Figure 7A, the USP18 protein was expressed at lower levels in EN DLBCL tissues than in LN DLBCL tissues. Compared to that in patients with LN DLBCL, the H-score of USP18 in the tumor tissues of EN DLBCL patients was significantly lower (2.125 vs 5.625, P < 0.01) (Figure 7B). The IHC results further confirmed that the downregulation of USP18 expression was associated with EN DLBCL.

To show the distribution of aDCs in DLBCL tissues, we also detected the expression of the aDC marker CD83 on tumor biopsies of DLBCL patients. As shown in Figure 7C, the expression of CD83 was distributed throughout the tissues and was lower in EN DLBDL tissues than in LN DLBCL tissues.

Furthermore, to identify whether aDCs express USP18, we performed immunofluorescence double staining of USP18 with CD83 in DLBCL tissues. As shown in Figure 7D, USP18 was coexpressed with CD83, which indicated that DCs also expressed USP18.

Data preparation and analysis of DEGs

Gene expression profiles and clinical characteristics of primary DLBCL samples were downloaded from TCGA (https://portal.gdc.cancer.gov/). HTseq-count and fragments per kilobase of exon per million reads mapped (FPKM) profiles of DLBCL samples, including 25 LN DLBCL and 21 EN DLBCL samples, were assembled. Immune-related genes were collected from the ImmPort database (https://www.immport.org/) [31]. Cancer-related TFs and ChIP-Seq data were retrieved from the Cistrome Cancer database (http://cistrome.org/) [32]. To identify DEGs between LN and EN DLBCL, the edgeR method was applied [33]; P < 0.05 and log (fold change) > 1 or < -1 were set as the cut-offs. Volcano plots and heatmaps were generated to show DEGs. Finally, GO and KEGG enrichment analyses of DEGs were performed to reveal the potential mechanism of EN DLBCL.

Identification of prognostic immune genes and construction of the predictive model

Volcano plots and heatmaps were created to illustrate the expression of DEIGs, which were extracted from the previous DEG and immune-related gene lists. Then, univariate Cox regression analysis was applied to identify the prognostic immune genes based on DEIGs and clinical information, with cut-offs of P < 0.05 and log (fold change) > 1 or < -1.

To assess the significance of each prognostic immune gene with a β value, which was the regression coefficient of integrated genes in the model, multivariate Cox regression analysis was carried out. The significant factors in the univariate Cox regression analysis were sent to the multivariate Cox regression analysis. The following formula was used to calculate the risk score:

$$\begin{array}{l}\text{Risk Score }=\text{ β}1\times \text{DEIG}1+\text{β}2\times \text{DEIG}2+\\ \text{β}3\times \text{DEIG}3\dots \dots +\text{βn}\times \text{DEIGn}\end{array}$$

In the formula, “n” is the number of prognostic immune genes in the model. “β” is the regression coefficient of each integrated gene. “DEIG_n” is the expression level of each integrated gene. Based on the model, patients were reordered and divided into high- and low-risk groups with the median risk score. To avoid model overfitting, Lasso regression and Schoenfeld residuals tests were performed. The AUC was applied to evaluate the accuracy of the model. Kaplan-Meier survival analysis was performed to compare patient survival between the two risk groups. Next, risk curve, survival state-related scatterplot and heatmap of prognostic immune genes were plotted based on the risk score.

Identification of key immune genes

Volcano plots and heatmaps were created to show the expression of DETFs, which were obtained by intersecting DEGs and cancer-related TFs. Then, to reveal the regulations and associations between DETFs and prognostic immune genes, Pearson correlation analysis was conducted, and only regulatory pairs with a correlation coefficient > 0.300 and P < 0.001 were selected for the next analysis. The intersection of prognostic immune genes in the above regulatory pairs and differentially expressed between EN and LN DLBCL by the Wilcoxon test was performed, as shown in the Venn plot. The immune gene in the regulatory pair with the highest coefficient and differentially expressed between EN and LN DLBCL was recognized as the key immune gene.

Validation of the regulatory mechanism between the key TF and immune gene

The regulatory mechanism between the key TF and immune gene was verified by ChIP-Seq. Two algorithms (JASPAR [34] and ENCODE transcription factor targets) were utilized to illustrate the transcriptional regulation patterns between LHX2 and USP18 to further confirm our hypothesis. LHX2 ChIP-Seq data from an in vitro cell line in the Cistrome database were downloaded to validate the transcriptional regulation patterns of USP18.

Identification of potential downstream KEGG pathways and immune gene sets

To determine the pathways downstream of key immune genes, GSVA was performed to identify differential KEGG pathways between EN and LN DLBCL. GSVA was implemented using the “gsva” package of R and under default settings except for “RNAseq = TRUE”. The GSVA algorithm accepted input from a gene expression matrix (log2-normalized RNA-seq count data) and a specific set of genes. The final output was a data matrix corresponding to each sample with each gene set. Pathways with P < 0.05 were selected and displayed. Pearson correlation analysis was used to uncover the relationship between the key immune genes and ENI-related signaling pathways, as shown by a coexpression heatmap. GSEA was also used to identify ENI-related signaling pathways [35]. Pathways with P < 0.05 was selected. The overlapping KEGG pathways from both GSEA and GSVA, illustrated by a Venn plot, were recognized as potential downstream pathways. The correlation between key immune genes and potential downstream pathways was fitted by linear regression.

ssGSEA was applied to identify immune responses in DLBCL from 29 immune gene sets that were overexpressed in the tumor microenvironment [11, 36]. Pearson correlation analysis was performed to illustrate the relationship between key immune gene and immune gene sets, as shown by the coexpression heatmap. Immune gene sets with P < 0.05 were selected and displayed. The top three correlations between key immune gene and immune gene sets were fitted by linear regression.

Online database validation and construction of regulation network including key TF, immune gene, KEGG pathways, and immune gene sets

For further annotation of identified TF, biomarker, immune gene sets, and signaling pathways, several online databases were used to detect gene and protein expression level. UALCAN [37], UCSC xena [38], Linkedomics [39], Gene Expression Profiling Interactive Analysis (GEPIA) [40], cBioportal [41] and Oncomine [42] were applied to validate the association between gene expression and clinical significance in tissue level in DLBCL. Furthermore, we used Cancer Cell Line Encyclopedia (CCLE) [43] to verify the gene expression in cellular level in DLBCL. Then the human protein altas [44] were applied to show the protein expression level in normal tissue. Finally, String [45] displayed the interaction network among LHX2, USP18 and the downstream pathway.

To show our results more clearly, a network based on the interaction among key TF, immune gene, KEGG pathways and immune gene sets was built by Cytoscape 3.7.1 [46]. Finally, EN DLBCL- related hypothesis built on the bioinformatics was displayed by signaling diagram.

Immunohistochemistry and immunofluorescence double staining

IHC and immunofluorescence staining were conducted according to standard methods on EN and LN DLBCL biopsies. Briefly, 5-μm formalin-fixed and paraffin-embedded (FFPE) sections were deparaffinized and hydrated. The sections were incubated overnight at 4° C in a humidified slide chamber with primary antibodies against USP18 (1:200, ab115618, Abcam) and CD83 (1:50, ab205343, Abcam). Finally, to assess the percentage of positive tumor cells, all the IHC slides were viewed and given a histochemistry score (H-score).

$$H-score={\displaystyle \sum pi(i+1)}$$

i represents the intensity score, and pi is the percentage of cells with that intensity. All immunofluorescence slides were observed with a confocal laser scanning microscope.

Statistical analysis

For descriptive statistics, the continuous variables in normal distribution were expressed as mean ± standard deviation (SD) while the median (range) was used in abnormal distribution. Classified variables were expressed by counts and percentages. Only two-tailed P < 0.05 was considered statistically significant. All statistical analysis was performed using R version 3.5.1 (Institute for Statistics and Mathematics, Vienna, Austria; https://www.r-project.org).

Ethical review committee statement

This study was approved by the Ethics Committee of Shanghai Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine.