Molecular characterization of EBV-associated primary pulmonary lymphoepithelial carcinoma by multiomics analysis

Abstract

Background: Primary pulmonary lymphoepithelial carcinoma (pLEC) is a subtype of non-small cell lung cancer (NSCLC) characterized by Epstein-Barr virus (EBV) infection. However, the molecular pathogenesis of pLEC remains poorly understood.

Methods: In this study, we explored pLEC using whole-exome sequencing (WES) and RNA-whole-transcriptome sequencing (RNA-seq) technologies. Datasets of normal lung tissue, other types of NSCLC, and EBV-positive nasopharyngeal carcinoma (EBV+-NPC) were obtained from public databases. Furthermore, we described the gene signatures, viral integration, cell quantification, cell death and immune infiltration of pLEC.

Results: Compared with other types of NSCLC and EBV+-NPC, pLEC patients exhibited a lower somatic mutation burden and extensive copy number deletions, including 1p36.23, 3p21.1, 7q11.23, and 11q23.3. Integration of EBV associated dysregulation of gene expression, with CNV-altered regions coinciding with EBV integration sites. Specifically, ZBTB16 and ERRFI1 were downregulated by CNV loss, and the FOXD family genes were overexpressed with CNV gain. Decreased expression of the FOXD family might be associated with a favorable prognosis in pLEC patients, and these patients exhibited enhanced cytotoxicity.

Conclusion: Compared with other types of NSCLC and NPC, pLEC has distinct molecular characteristics. EBV integration, the aberrant expression of genes, as well as the loss of CNVs, may play a crucial role in the pathogenesis of pLEC. However, further research is needed to assess the potential role of the FOXD gene family as a biomarker.

Keywords: Biomarker; EBV; Lung cancer; Multiomics; Tumor microenvironment.

Fig. 1

Overview, design and summary statistics. This flowchart describes the inclusion and exclusion of all samples. pLEC Pulmonary lymphoepithelial carcinoma, NPC Nasopharyngeal carcinoma, LUAD Lung adenocarcinoma, LUSC Lung squamous cell carcinoma, WES whole exome sequencing, RNA-seq RNA whole transcriptome sequencing

Fig. 2

Prognosis of pLEC by Kaplan-Meier survival analysis with log-rank test. A Cox-model hazard ratios (x-axis) for Clinical features (y-axis)). Horizontal bars indicate 95% confidence intervals. B The probability of OS in patients with pLEC compared with patients with LUSC, LUAD, and NPC (P = 0.001). C Cophenetic correlation from NMF analysis of pLEC and poorly differentiated LUSC. D The NMF consensus matrix shows that the cluster analysis identified 2 subtypes. P < 0.05 indicates a significant difference. NPC nasopharyngeal carcinoma, LUAD lung adenocarcinoma, LUSC lung squamous cell carcinoma

Fig. 3

Mutation spectrum and CNV landscape of pLEC and comparison with other cancer types. A Mutational landscape of 15 pLEC samples. The dark blue, red, orange, khaki, dark green, light blue, light green, and purple represent missense, nonsynonymous, frameshift deletions, frameshift insertions, splice site mutations, and in-frame insertions, respectively, in-frame deletions, and multiple mutations. The bar graph on the right shows the mutation frequency in the population, and the upper bar graph shows the frequency of non-synonymous mutations in the sample. B High-frequency CNVs of fragments in pLEC. CNV gain (red) is shown on the left and CNV loss (blue) is shown on the right. The figure shows significant amplification or deletion of chromosomes from 1 (top) to 22 (bottom). The green line represents the cutoff value for significance (q = 0.1). C Distribution of non-silent mutation rates in LUSC (red), LUAD (yellow), NPC (cyan) and pLEC (purple). The black lines in the scatterplot represent the median and upper and lower quartiles of sample mutation frequency in each cancer. All cancers were statistically tested by the Wilcoxon rank sum test method. A P value less than 0.05 was defined as a significant difference. D Mutation profiles of the six mutation types for each cancer type. E Mutation spectrum of six mutation types for each cancer type. F Comparison of genomic mutational profiles between pLEC, LUAD, LUSC, and NPC. The bars at the bottom of the graph represent different tumor types. Color blocks represent different types of base changes, and the mutation frequency of the population is shown on the right. G Comparison of CNA in pLEC and other cancer types. H Pathway Alteration Frequencies. The fraction of altered samples for each pathway and tumor subtype. Pathways are ordered by decreasing median frequency of change. Increasing the color intensity reflects a higher percentage. Average mutation counts for each cancer subtype are also provided. Average mutation counts for each cancer subtype are also provided

Fig. 4

Transcriptome comparison of pLEC with other cancer types. A Differential expression analysis of RNA-seq data, where red and blue correspond to up- and down-regulation of pLEC samples relative to NPC samples, respectively. B Evaluation of tumor stemness score. The Wilcoxon rank-sum test was used for statistical testing of the two groups of data. C Sample Regulatory Cell Death Type Score. D Comparative analysis of immune microenvironment between pLEC and NPC, LUSC and LUAD. The upper bars indicate different assessment methods, blue indicates significant downregulation of pLEC and red indicates significant upregulation. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Fig. 5

EBV integrated analysis. A Landscape of EBV integration breakpoints across the genome. Chromosomes are numbered and represented by a series of colors. B Distribution of EBV integration breakpoints across chromosomes. The x-axis is the number of chromosomes, and the y-axis is the number of integration breakpoints. C mRNA cluster analysis of EBV integration sites. The calculation formula of the number (N) of each block is: N = log2 (mRNA expression of a gene in tumor tissue/mRNA expression of a gene in control tissue). Blue represents down-regulated mRNA and red represents up-regulated mRNA

Fig. 6

Integrated analysis of genome and transcriptome and prognostic analysis. A Integrated analysis of transcriptome and genomic CNVs. Each column represents a sample, and each row represents genes amplified by high frequency CNVs. At the top of the chart are annotations of sample clinical information, representing OS time as well as clinical stage. The graph shows the mRNA expression of this gene relative to the control log2-fold change in mRNA expression in cancer samples. FOXD4L2, FOXD4L4, and FOXD4L5 family genes were amplified and upregulated in 100% (14/14) of pLEC patients; mRNA log2-fold changes of FOXD4L2, FOXD4L4, and FOXD4L5 ranged from 0.21 to 1.8, 0.47 to 2, and 0.33 to 1.91, respectively ZBTB16, ERRFl1, GPR62 gene copy number loss and down-regulation. B Kaplan-Meier survival analysis of genes with copy number loss or gain. Statistical significance was estimated by two-sided log-rank test

Fig. 7

Mechanisms underlying an FOXD4L4-associated immune-suppressive tumor microenvironment. A,B FOXD4L4 expression was significantly negatively correlated with infiltrating levels of CD8 + T cells. C The FOXD4L4^high group had significantly greater enrichment scores for activated CD8 T cells and cytotoxic cells. Note: *p < 0.05, **p < 0.01, ***p < 0.001. D-F Boxplot showing cell death scores within FOXD4L4^high and FOXD4L4^low groups