Genetic and immune characteristics of multiple primary lung cancers and lung metastases.

BACKGROUND: To explore the genetic and immunophenotyping heterogeneities between patients with intrapulmonary metastasis (IPM) or multiple primary lung cancer (MPLC).
METHODS: Whole exome sequencing (WES) and transcriptome sequencing (RNA-seq) were performed on the tissue and blood samples of IPM and MPLC patients to comprehensively analyze the clonal evolution, molecular typing and immunophenotyping.
RESULTS: There was no significant difference in genetic mutation, tumor mutational burden (TMB) value and mutant allele tumor heterogeneity (MATH) value between IPM and MPLC patients. Notably, the loss of heterozygosity (LOH) of human leukocyte antigen (HLA) appeared in all IPM patients, while there was also no significant difference between the two groups. In addition, expression of immune checkpoint-related genes including CTLA-4, BTLA, TIGIT and HAVCR2 in the MPLC group was significantly higher than those in IPM group. At the same time, 86 differentially expressed genes (DEGs) were observed between IPM and MPLC patients with transcriptome sequencing, of which 56 DEGs were upregulated and 30 were downregulated in the IPM group compared with the MPLC group. The cluster analysis revealed that the 86 DEGs could be distinguished in IPM and MPLC samples. Moreover, only the infiltration levels of CD56dim natural killer cells in the IPM group was significantly higher than that in the MPLC group, and the infiltration levels of the remaining 27 immune cell subsets were similar in both groups.
CONCLUSIONS: IPM and MPLC are roughly similar in genetic and immune characteristics indicating that genomics alone may not be able to effectively distinguish between IPM and MPLC, which still needs to be comprehensively evaluated with clinical manifestations, imaging, and pathological characteristics.

INTRODUCTION

Lung cancer is the leading cause of human death worldwide, and about 1.76 million people die of lung cancer every year. Lung adenocarcinoma (LUAD) is one of the main types of lung cancer, accounting for approximately 50% of all lung cancer cases, and its five‐year survival rate is about 50% in the early stages, but only about 20% in the late stages. During the treatment of primary lung cancer, the patient may develop another lung cancer, which is considered as multiple primary lung cancer (MPLC) or intrapulmonary metastasis (IPM). MPLC refers to the simultaneous or consecutive occurrence of two or more primary lung cancers in the lung of same patient. When multiple tumors are histologically similar, it is difficult to distinguish MPLC from IPM. If secondary lung cancer or IPM occurs in a place which previously received radiotherapy, it will be more difficult to distinguish MPLC and IPM due to morphological changes caused by radiotherapy. All of the above has brought new challenges to the current clinical management of MPLC. Notably, the treatment strategies and prognosis of IPM and MPLC are different. MPLC is mainly performed with surgical resection, while IPM is recommended for systemic chemotherapy instead of surgery. Thus, distinguishing IPM and MPLC has clinical significance in improving the accuracy of prognosis assessment as well as therapeutic intervention.

At present, the diagnosis of MPLC is mostly according to the eighth edition American Joint Committee on Cancer (AJCC) staging manual, which is based on the clinical, histopathological and molecular diagnoses of lung tumors. However, there are still many secondary lung cancers misdiagnosed. In recent years, next‐generation sequencing (NGS) has been widely used in the medical field and tumor research, enabling deeper and more comprehensive research on diseases from genome to transcriptome to be conducted. Murphy et al. used NGS to conduct a diagnostic lineage test based on genomic rearrangements from mate‐pair sequencing and found that MPLC had no common genome rearrangements, but IPM did, so that they could distinguish MPLC from IPM. Chen et al. suggested that the combination of histological characteristics and genetic alterations might be an effective method for the diagnosis of MPLC and IPM, and MPLC exhibited unique molecular characteristics which might help distinguish between patients with MPLC or IPM. Therefore, NGS might help distinguish between IPM and MPLC by analyzing a large amount of biological information.,

In this study, we comprehensively analyzed the genetic and immune characteristics of MPLC and IPM using NGS to explore the differences in molecular and immunotyping heterogeneity, aiming to explore whether NGS can be used to effectively distinguish MPLC from IPM, and thereby accurately guiding management of the disease and improving patient prognosis.

METHODS

Patient data

The study was approved by the Ethics Committee of the Affiliated Hospital of Qingdao University, and all patients signed written informed consents.

A total of 11 multifocal LUAD patients in the Affiliated Hospital of Qingdao University were enrolled in this study. The inclusion criteria were as follows: pathologically diagnosed as LUAD; each patient had two or more resectable lesions; the patients had not previously received chemotherapy, radiotherapy or other treatments prior to surgery. Exclusion criteria were as follows: patients with pathological diagnosis other than adenocarcinoma; cases had previous history of other cancers or cancer‐related treatments. Blood samples and tumor tissues were collected from each patient and whole exome sequencing (WES) and RNA sequencing was performed. Tumor tissue samples collected from surgical specimens were fixed in formalin and embedded in paraffin.

DNA extraction and whole exome sequencing (WES)

Genomic DNA from blood samples or formalin‐fixed and paraffin‐embedded (FFPE) tissues was extracted using Blood Genomic DNA Extraction Kit (Tiangen Biochemical Technology [Beijing] Co., Ltd.) or BLACK PREP FFPE DNA kit (Analytik Jena AG), respectively, according to the manufacturer's protocol, and stored at −20°C. The fragmentation of genomic DNA was performed with Covaris M220 focused ultrasonicator (Covaris Inc.). Then, DNA library was constructed by KAPA Hyper Prep Kit (Illumina platforms) (KAPA Biosystems) and captured using NimbleGen SeqCap EZ Exome Library (Roche), followed by sequencing using Illumina NovaSeq 6000 platform (Illumina).

WES data processing and analysis

Sequencing reads were aligned to the human reference genome (Hg19, NCBI Build 37.5) using the Burrows‐Wheeler Aligner (version 0.7.17) after removing low quality reads. Duplicate reads were marked by Picard toolkit (version 2.1.0) and realigned using the Genome Analysis ToolKit (version 3.7). Single nucleotide variants (SNV) and short indels in tumor tissue samples were identified by Mutect2 and variants were annotated using ANNOVAR. The identified tumor‐related mutated genes were then classified into 10 signaling pathways and subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses.

Tumor mutational burden (TMB) analysis

Nonsynonymous mutations (SNV and Indel) in the coding region for a given gene were selected for the assessment of tumor mutational burden (TMB) in SCLC. While driver gene mutations and hotspot mutations included in the ExAC/COSMIC database were filtered out, nonsynonymous, frameshift, and stopgain mutation sites within the exonic regions; depth ≥ 40X; reads ≥ 5; and frequency ≥ 5 were chose as TMB candidate sites, and TMB was calculated according to the following formula.

Mutant allele tumor heterogeneity (MATH) analysis

The tumor/normal ratio was determined for each tumor sample and the percentage of tumor cells in each tumor tissue sample was found not to differ significantly. Mutation of tumor samples were corrected using the matched BC samples and tumor/normal ratio from each patient. SNV mutations were selected according to the following criteria: freq ≥ 5; depth ≥ 50X; mutations within the exonic regions were preserved and synonymous mutations were filtered out; mutation sites with the frequencies of more than 10% in BC samples were filtered out; mutation sites with a frequency of 10 times higher than that of less than 10% for the sites in the BC sample were kept. The MATH value of each tumor was calculated from the median absolute deviation (MAD) and the median of its mutant‐allele fractions at tumor‐specific mutated loci: MATH = 100*MAD/median.

Transcriptome sequencing (RNA‐seq)

Total RNA was extracted from the FFPE tumor tissues using RNeasy FFPE Kit (Qiagen Inc) and stored at −80°C. The cDNA library was created using the SMARTer Stranded Total RNA‐Seq Kit v2 (Takara Bio Inc) according to the manufacturer's instructions. After PCR enrichment and purification of adapter‐ligated fragments, RNA sequencing was performed using Illumina NovaSeq 6000 platform.

RNA‐seq data processing and analysis

Clean data were aligned to the reference genome which is available from the ENSEMBL website (http://www.ensembl.org/index.html) using HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml). The expression level of each gene was then determined by HTSeq (http://www.huber.embl.de/users/anders/HTSeq/doc/overview.html). The quantification of gene expression was performed based on fragments per kilobase of exon model per million mapped reads (FPKM). The differentially expressed genes (DEGs) were identified using DESseq package based on |log2 (Fold Change) | ≥2 and p‐value<0.05. Then, DEGs were compared with gene ontology (GO) and the KEGG databases, and the enrichment analysis of the DEGs was conducted using Metascape (http://metascape.org/). In addition, the infiltration of 28 tumor‐infiltrating lymphocytes (TILs) in each sample was analyzed by ssGSEA method as previously described.

Statistical analysis

The data was analyzed using R 3.6.3 software, and the differences between different groups were analyzed using Wilcoxon Rank‐Sum test or Fisher's exact test. p < 0.05 was considered statistically significant.

RESULTS

WES data analysis of tissue samples

According to the phylogenetic tree analysis, there were four IPM and seven MPLC patients among the 11 LUAD patients (Figure S1). The clinical characteristics of all patients are shown in Table S1. The SNV analysis showed that EGFR mutation was the most frequent (38%) in LUAD patients, and that mutated genes were mainly involved in signaling pathways including p53_, PI3K_ and PTK_RAS_pathways (Figure S2). However, there was no significant difference in mutated genes and their involved pathways between IPM and MPLC patients (p > 0.05) (Figure S3). Moreover, the values of TMB and MATH were not significantly different between IPM and MPLC patients (p > 0.05) (Figure S4).

The germline HLA‐I molecular alleles of patients were classified using HLA‐HD, while there was no difference between IPM and MPLC groups (Figure S5). In addition, the loss of heterozygosity (LOH) of HLA in tumor specimens was further analyzed, and a higher frequency of HLA‐LOH in the IPM group compared with the MPLC group was observed (100% vs. 42.9%, p > 0.05) (Figure FIGURE 1).

FIGURE 1

The loss of heterozygosity of human leukocyte antigen (HLA) in intrapulmonary metastasis (IPM) and multiple primary lung cancer (MPLC) patients shown by phylogenetic tree analysis

Identification and enrichment analysis of differentially expressed genes

A total of 86 DEGs were identified by comparing the gene expression of IPM and MPLC tumor samples using DEGseq software (adjust p < 0.05), and the cluster analysis showed that these DEGs distinguished IPM and MPLC samples well (Figure FIGURE 2a,b). The GO and KEGG enrichment analysis revealed that these DEGs were mainly enriched in the metabolic process or the signaling pathway of cytokine‐cytokine receptor interaction, respectively (Figure FIGURE 2c,d and Figure S6). In addition, the expression levels of immune checkpoint‐related genes in two groups was analyzed, which found that the expression of CTLA‐4, BTLA, TIGIT and HAVCR2 in the MPLC group was significantly higher than those in the IPM group (p < 0.05) (Figure FIGURE 2e).

FIGURE 2

Identification and enrichment analysis of differentially expressed genes. (a) Volcano plot of differentially expressed genes (DEGs) in intrapulmonary metastasis (IPM) and multiple primary lung cancer (MPLC) patients. (b) Cluster heatmap of the DEGs in IPM and MPLC patients. (c) Gene ontology (GO) enrichment analysis was performed with DEGs. (d) The Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis was performed with DEGs. (e) The expression levels of immune checkpoint‐related genes in IPM and MPLC patients

Analysis of tumor‐infiltrating lymphocytes

At the same time, the signature scores of nine immune features of tumor samples, gene‐expression profiling (GEP) and TILs scores were calculated, but no significant differences were found in the IPM and MPLC groups (p > 0.05) (Figure 3a,b). The infiltration of 28 immune cell subgroups in the tumor microenvironment were also detected by the ssGSEA method. As shown in Figure 3c, the infiltration level of CD56dim natural killer cells in the IPM group was significantly higher than that in the MPLC group (p < 0.05), while there was no significant difference in other cell subsets between both groups. Subsequently, the immune status of tumor samples was divided into high degree of immune cell infiltration (TILs‐high, n = 14) and low degree of immune cell infiltration (TILs‐low, n = 10) according to the clustering results of the ssGSEA scores of 28 immune cell subgroups of all samples (Figure S7). Further, according to the immune status of the tumor samples, the patients were divided into the homogeneous group (patients with one immune status) and the heterogeneous group (patients with two immune status). Among 11 LUAD patients, three patients were homogeneous (IPM, n = 1; MPLC, n = 2), and eight were heterogeneous (IPM, n = 3; MPLC, n = 5). Analysis of immune activity found that compared with the heterogeneous group, the homogeneous group had a higher CD8 and Th1 cell infiltration levels, higher CYT and GEP scores as well as higher MHC‐II molecular expression level (Figure S8).

FIGURE 3

Analysis of tumor‐infiltrating lymphocytes between intrapulmonary metastasis (IPM) and multiple primary lung cancer (MPLC) patients. (a) Comparison of the signature scores of nine immune features between IPM and MPLC tumor samples. (b) Comparison of the tumor‐infiltrating lymphocytes (TILs) and gene‐expression profiling (GEP) scores between the IPM and MPLC tumor samples. (c) Comparison of the expression of 28 immune cell subsets between the IPM and MPLC groups. ns p > 0.05, *p < 0.05

DISCUSSION

With the rapid development of imaging and various preoperative diagnostic methods, the detection rate of multiple lung tumors is getting higher, but the distinction between MPLC and IPM is still a difficult problem. Currently, the progress of lung cancer genomics has profoundly changed our understanding of lung cancer at the molecular level. EGFR mutations are common in primary lung tumor, but rarely in the metastases. In this study, higher EGFR mutation was also found in LUAD patients, mainly involved in the signaling pathways including p53_, PI3K_ and PTK_RAS_pathways, while there was no significant difference between the IPM and MPLC groups. Notably, in this study, LOH of HLA appeared in all IPM patients. HLA plays an important role in immune response, and is mainly involved in the processing and presentation of antigens. Studies have shown that when LOH occurs at HLA sites, it may promote immune evasion and lead to immunotherapy resistance, which may be the reason for the poor prognosis of IPM.,

In addition, by analyzing the expression levels of immune checkpoint‐related genes in both groups, it was found that the expression of CTLA‐4, BTLA, TIGIT and HAVCR2 in the MPLC group was significantly higher than those in the IPM group. Abnormal expression and function of immune checkpoint molecules is one of the important causes of cancer. It has been reported that tumor cells can activate immune checkpoints to prevent antigens from being presented to T cells, thereby suppressing the immune function of T cells and leading to immune evasion. Therefore, we speculated that MPLC patients might benefit more from immunotherapy compared with IPM patients.

Studies have found that TILs infiltration plays an important role in the tumor microenvironment,, , which helps tumor cells escape immune surveillance in the late stage. The density, type and proportion of TILs reflect the immune status of the local tumor microenvironment. Thus, we examined the infiltration of 28 immune cell subgroups in the tumor microenvironment using ssGSEA method to explore the difference of TILs between IPM and MPLC patients. The results showed that the infiltration level of CD56dim natural killer cells in the IPM group was significantly higher than that in the MPLC group. Natural killer cells are an important part of the immune system, which can inhibit tumor progression and blood‐borne metastasis. Bauernhofer et al. pointed out that the CD56dim natural killer cells in cancer patients were targeted for apoptosis, leading to low activity of natural killer cells. In addition, MPLC and IPM patients were divided into homogeneous and heterogeneous groups based on the immune status of the tumor samples, respectively, while there was still no significant difference in both groups. However, immune activity analysis revealed that the homogeneous group had higher CD8 and Th1 cell infiltration levels, higher CYT and GEP scores as well as higher MHC‐II molecule expression level compared with the heterogeneous group.

In summary, there were relatively small differences in genetic heterogeneity and immune heterogeneity in IPM and MPLC patients, and it was difficult to effectively distinguish IPM and MPLC patients by NGS sequencing only. Schneider et al. stated that comprehensive genomic and morphological assessment was feasible to identify the MPLC. Saab et al. also pointed out that 65% of MPLC patients could be identified based on clinical manifestations, imaging data, and morphology, but when there were more than three lung tumors with similar growth patterns and lack of major growth patterns, it was difficult to distinguish them based on the above evidence alone. However, 94% of patients can be identified when combining the patient's morphological characteristics and genomics. Thus, genomics are still an auxiliary tool to distinguish between MPLC and IPM.

In conclusion, IPM and MPLC are very similar in genetic and immune characteristics, so although it is difficult to distinguish between MPLC and IPM effectively based on genetic and immune characteristics they still need to be comprehensively evaluated with clinical manifestations, imaging, and pathological characteristics.

CONFLICT OF INTEREST

The authors confirm that there are no conflicts of interest.

Figure S1. The results of phylogenetic tree analysis in 11 LUAD patients.

Click here for additional data file.

Figure S2. Genetic variation analysis of tissue samples. (a) The single nucleotide variation (SNV) landscape of LUAD patients. (b) Map of 10 major signaling pathway mutations.

Click here for additional data file.

Figure S3. The gene mutation characteristics of (a) IPM and (b) MPLC group tumor tissue samples.

Click here for additional data file.

Figure S4. Comparative analysis of (a) TMB and (b) MATH in IPM and MPLC patients. ns p > 0.05.

Click here for additional data file.

Figure S5. Molecular typing analysis of germline HLA‐I molecular alleles. (a) Molecular typing analysis of germline HLA‐I molecular alleles in all patients. (b) Molecular typing analysis of germline HLA‐I molecular alleles in MPLC and IPM patients.

Click here for additional data file.

Figure S6. Identification and enrichment analysis of differentially expressed genes. (a) The GO enrichment analysis was performed with 56 upregulated DEGs. (b) The GO enrichment analysis was performed with 30 downregulated DEGs. (c) The KEGG enrichment analysis was performed with 56 upregulated DEGs. (d) The KEGG enrichment analysis was performed with 30 downregulated DEGs.

Click here for additional data file.

Figure S7 The cluster analysis was performed with the ssGSEA scores of 28 immune cell subgroups of all samples.

Click here for additional data file.

Figure S8. Analysis of Tumor‐infiltrating lymphocytes between homogeneous and heterogeneous groups. (a) Comparison of the expression of 28 immune cell subsets between homogeneous and heterogeneous groups. (b) Comparison of the signature scores of 9 immune features between homogeneous and heterogeneous groups. (c) Comparison of the TILs score, CYT score and GEP score between homogeneous and heterogeneous groups. ns p > 0.05, *p < 0.05, **p < 0.01.

Click here for additional data file.

Table S1. Clinical and histopathological characteristics of LUAD patients.

Click here for additional data file.