Conducted between July 2017 and December 2019, this prospective observational cohort study focused on NSCLC patients undergoing surgical treatment (I–III) at a single tertiary hospital. Surgery consisted of anatomical pulmonary resection and, systematic lymph node dissection for clinical stages II and IIIA and lobe-specific lymph node dissection for stage I. For the LB analyzes, the peripheral blood samples were collected at baseline (CTC1) and one month after surgery (CTC2). Tumor tissue was extracted during surgery, tumoral and peritumoral samples were obtained through a sterile field and placed in RNA Later tubes (AM7020, ThermoFisher Scientific, Waltham, Massachusetts, USA) for preservation.

The pathologic stage was determined by the 8th edition of the international tumor-node-metastasis (TNM) system and, according to the World Health Organization classification, the histologic subtypes of lung cancer were assigned. After surgery, patients were discussed in the hospital's thoracic multidisciplinary oncology committee. Adjuvant chemotherapy was indicated in stages IIb and III following the ESMO guidelines, and it was considered in patients with stage IIA.

The follow-up protocol consisted of surveillance every 6 months for the first 3 years, which included a chest and abdominal CT scan, with optional FDG-PET if required. Then, the follow-up was yearly. Clinical outcomes regarding overall survival (OS) and progression-free survival (PFS) were estimated using the Kaplan–Meier method and both were calculated at 5 years after surgical treatment.

The study, adhering to the Declaration of Helsinki and approved by regional institutional review boards, obtained written informed consent from all participants before enrollment.

Matched tumoral/peritumoral tissue biopsies were collected from patients that underwent surgical treatment. gDNA was extracted using the NucleoSpin® Tissue kit (740952.50, Macherey-Nagel, Düren, Germany), following the manufacturer's instructions. Briefly, 25mg of tissue was used for the extraction. Tissues were homogenized using the TissueLyser LT instrument (Qiagen, Venlo, Netherlands) for 5min at 50oscillations/s. After that, samples were lysed with Proteinase K solution at a final concentration of 2.7mg/ml at 56°C overnight in a shaker. After purification and precipitation, DNA was bound to a silica column, washed twice, and the column was dried out. DNA was eluted in a final volume of 75μl. Concentration was measured using a NanoDrop™ 2000 instrument (ThermoFisher) at a wavelength of 260nm.

Sequencing libraries were constructed using the KAPA HyperCap workflow v3.0 (Roche, Basel, Switzerland) with a customized gene panel, following the manufacturer's protocol. This customized panel was developed to cover 50 genes (exons) involved in four different biological processes: tumor promotion (ROS1, ALK, TP53, APC, RB1, RET, FLT1, PTEN, STK11, TET2, CTNNB1, MLH1, EGFR, DDR2, SLC22A18, AKT1), Chronic Obstructive Pulmonary Disease (COPD) associated genes (ADAMTSL3, ADAM19, PTCH1, ADGRG6, HTR4, FAM13A, CHRNB4, RIN3, TGFB2, IREB2, SERPINA1, CHRNA3, STN1, CHRNA5, AGER), lung development (FGF10, FOXA1, FOXA2, FOXP2, FOXP4, HOXA2, HOXA5, HOXB2, PITX2, WNT2, ALX1, NKX2-1) and inflammation (CD163, S100A4, NLRP3, MUC5AC, EGR1, TLR4, FPR2), spanning approximately 140kb. It was tested on I–IIIB stage NSCLC to elucidate the contribution of non-tumoral gene mutations on cancer disease development, intercepting cancer. Libraries were quantified with Qubit 4.0 Fluorometer using dsDNA HS Assay kit. Individual libraries were diluted and pooled together in equal amounts (by mass, equimolarity) for obtaining equal numbers of sequencing reads (HiSeq X System, Illumina, San Diego, California, USA).

Once FASTQ files were demultiplexed and retrieved from the sequencer, files were ready for data processing. The GATK best practices-based on the Somatic Variant Pipeline (https://www.clinbioinfosspa.es/content/somatic-variant-pipeline) from the Clinical Bioinformatics Area from Fundación Progreso y Salud (Andalusian Govern – Health and Family Counseling) were followed. Briefly, raw data quality was evaluated using the FastQC tool (v0.11.9, Babraham Bioinformatics, Cambridge, UK). High quality reads were aligned against the human reference genome (GRCh38) using Burrow–Wheeler Aligner algorithm (https://github.com/lh3/bwa). After filtering steps, tumor with matched normal mode from Mutect2 pipeline from GATK (https://gatk.broadinstitute.org/hc/en-us) was run for variant calling. Generated filtered VCFs were annotated using VEP (v99, https://www.ensembl.org/info/docs/tools/vep/index.html). Clinical significance was assessed by following the American College of Medical Genetics (ACMG) criteria. R software (v4.2.1) and RStudio (2022.12.0+353 version) were used for file parsing and for plot generation. Maftools17 was used for oncoplot generation and somatic interactions plot.

Peripheral blood (10ml) was collected in EDTA BD Vacutainer® tubes and stored at room temperature for processing within the next 4h. CTCs were isolated by density gradient centrifugation, following a previously described protocol.12 Briefly, blood samples were centrifuged by density gradient with Ficoll Histopaque®-1119 (density=1.119g/ml) (Merck, Darmstadt, Germany) for 45min at 700×g with no break. The buffy coat where CTCs and PBMCs are located was carefully extracted. Isolation of CTCs was performed by positive immunomagnetic selection using the Carcinoma Cell Enrichment kit (130-108-339, Miltenyi Biotec, Bergisch Gladbach, Germany). This selection is based on a multi-cytokeratin antibody (CK3-6H5), that recognizes the cytoplasmic cytokeratins 7 and 8. Cellular samples were spun down onto slides that were stained and characterized by immunofluorescence.

The working area was marked with a hydrophobic pen (Dako). Each sample was hydrated with PBS 1× for 5min and incubated with anti-cytokeratin-FITC antibody, against cytokeratins 7 and 8 (130-060-301, Miltenyi Biotec), diluted 1:100 in PBS 1× – Tween 0.1%+BSA 2% (working solution) for 45min in the dark at room temperature. After this point, darkness was always maintained. After this incubation, samples were washed with PBS 1× for 5min. A blockade was performed with goat serum diluted at 1× in working solution for 1h. To detect the EZH2 protein, primary rabbit anti-EZH2 antibody (#5246S, Cell Signaling Technology) was diluted 1:100 in the working solution and was incubated for 45min. Then, samples are washed in PBS 1× for 5min. Secondary goat anti-rabbit IgG (Alexa Fluor 555, ThermoFisher) was diluted 1:200 in the working solution and was incubated for 45min to reveal EZH2 expression. Nuclei were counterstained using Hoechst33342 (Fluxion) at 1×. Slides were mounted using SlowFade™ Diamond Antifade Mountant (S36972, ThermoFisher) for visualization in Zeiss LSM 710 confocal microscope.

Evaluation of EZH2 in CTCs was analyzed after previous evaluation on A549 control cell line (CCL-185™, ATCC, Manassas, Virginia, USA) for antibody specificity (Supplementary Fig. 1). Because of the heterogeneity observed in patients with CTCs, the percentage expression of EZH2 was calculated as the number of EZH2+ CTCs divided by the total number of CTCs. This calculation aimed to quantify the heterogeneity of CTCs within patients.

CTCs were assessed both as a continuous variable (number of CTCs per 10ml of blood) and as a binary variable (presence/absence and a cut-off of ≥2 CTCs per 10ml). Cut-offs to distinguish low/high CTC numbers and EZH2 expression were determined using the Cut-off Finder web application (fit for mixture model).18 Bivariate analysis examined the relation between clinical features and CTC presence/absence. Clinical variables were evaluated for normality using the Shapiro–Wilk test.

The clinical characteristics were evaluated by Chi-squared (Gender, Procedure, Histology, Respiratory history, Stage, Adjuvancy, Recurrence and Vital status), Fisher's exact test (Tobacco, Resection, Lymph node status, Type of adjuvancy, Site of recurrence, and Cause of death), t-test (Age) and Mann–Whitney U test (PET – SUVmax–, Tumor size (cm), overall survival – OS and progression-free survival – PFS, both in months). Parametric variables are expressed as mean with standard deviation and non-parametric variables as median with interquartile range or down/upper quartile (Q1, Q3).

Wilcoxon Signed-Rank test compared the median of CTC numbers and EZH2 expression between sampling points, to study the dynamics of these variables. Univariate Kaplan–Meier (log-rank test) analyzed the influence of CTC and EZH2 cut-offs on PFS and OS. Fisher's exact test evaluated somatic interactions between gene pairs, while t-tests and Fisher's exact test were further used to compare mutation profiles concerning recurrence and CTC presence, respectively.

Statistical analyses were performed using IBM SPSS Statistics (version 25.0 for MacOS, IBM Corp., Armonk, NY, USA), GraphPad Prism (version 9.0.2 for MacOS, GraphPad software, La Jolla, CA, USA), and R software (v4.2.1) and RStudio (2023.03.0+386 version), with a p-value<0.05 determining statistical significance.

Out of the total of 76 patients, only 27 were eligible for inclusion in the NGS analyses due to meeting the required quality standards for library preparation. Average concentration of extracted gDNA was 316.53ng/μL (range: 43.85–946.85ng/μL). Minimal input ensured for library preparation was 500ng. After library preparation, sequencing yielded a mean coverage of approximately 500× along all the samples. From the 27 patients sequenced, only 4 did not pass variant calling filters and they were removed from subsequent analysis. The total number of variants identified for the 23 remaining patients were 143, with an average of 6 variants per patient (range: 1–25). Most of the variants (61.5%; 88/143) were non-coding and SNPs (91.6%), but only 34.3% (49/143) were known variants. Variants were submitted to ClinVar (Accession numbers: SCV003806285–SCV003806424). Regarding the clinical significance of the variants, 9.8% (14/143) were benign/likely benign, 83.2% (119/143) were variants of uncertain significance (VUS) and 7.0% (10/143) were pathogenic/likely pathogenic (Fig. 1A). Oncoplot in Fig. 1B summarizes data for the mutational profiles of each patient and gene. The number of mutations for each patient and the number of samples (patients) affected by mutations in that given gene are represented on both top and right plots, respectively. The top 5 most highly mutated genes were: TP53, FLT1, MUC5AC, EGFR and NLRP3, presented in 43%, 39%, 39%, 26% and 26% of the samples, respectively. Cancer stage, presence/absence of CTCs, expression of EZH2 and mutational signatures (transitions and transversions) are also provided for each patient (Fig. 1B).

Fig. 1.

Mutational landscape of the sequenced patients. (A) Descriptive flow of all the variants found in the 23 patients (Sankey plot) summarizing all the characteristics. (B) Oncoplot summarizing mutational profile for each patient and gene. Top bar plot represents the number of mutations for each patient filled depending on the consequence of mutation. Right side bar plot illustrates the number of samples affected by mutation in the given gene, filled depending on the consequence of mutation. Clinical stage of each patient is also present below sample names, with corresponding proportions of transitions/transversions.

(1.11MB).

Full description of the pathogenic, likely pathogenic variants and hot VUS found is described in Table 1. Both PCT74 and PCT75 patients recurred and had primary point mutations in EGFR and ROS1, respectively, and both had secondary point mutations in TP53 (according to allele frequencies). Contrarily, patients PCT33 and PCT76, had primary mutations in TP53 and secondary mutations in RB1, and they did not have recurrence, suggesting that secondary mutations in TP53 might have a stronger impact on disease recurrence.

Studying if a pair of genes were mutated simultaneously or if they were mutually exclusive, interactions of the somatic variants were analyzed. The top 5 most significant pairs of genes with co-occurring mutations were: CD163-TLR4, FGF10-FOXP2, ADAMTSL3-FLT1, ADAMTSL3-MUC5AC and MUC5AC-NLRP3 with p-values of 0.0056, 0.0118, 0.0142, 0.0144, 0.0183, respectively. There was only one significant pair of genes with mutually exclusive mutations: RIN3-TP53 (p=0.0457). Interactions among all the genes are represented in Fig. 2. Detailed data of significant pairs of genes highlighted in Fig. 2 can be found in Supplementary Table 2 with additional information for all the significant pairs of genes.

Fig. 2.

Somatic interactions between pairs of genes. Green color is for co-occurrence and pink for mutually exclusive, in boxes. Purple genes are tumor genes, red are inflammation genes, green are COPD-associated genes and black are lung development genes. Number in brackets represents the count of mutations for a given gene. *: p-Value<0.05 for Fisher's exact test.

(0.71MB).

Most variants were VUS (83.2%) and as they might be reclassified in the future as either pathogenic/likely pathogenic (P) or benign/likely benign (B), we studied their impact on recurrence by grouping them with either P or B. Two comparisons were done: P+VUS (pathogenic variants and VUS) and VUS+B (VUS and benign variants). P+VUS showed significant differences (p=0.0404) between no-recurrence and recurrence patients, with a mean number of variants of 4.8 (standard error of the mean – SEM=0.65) in the no-recurrence group versus a mean of 11 (SEM=5.54) variants in the recurrence group. On the other hand, VUS+B showed no significant differences (p=0.0597) between no-recurrence (mean of 5, SEM=0.71) and recurrence patients (mean of 11, SEM=7.09) (Fig. 3).

Fig. 3.

Differences in the number of variants between no recurrence vs. recurrence patients. Number of variants presented in both types of patients, grouped by clinical significance. P+VUS: Pathogenic variants and VUS, VUS+B: VUS and benign variants, P+B: Pathogenic variants and benign variants, Total: all types of variants. *: p-value<0.05, ns is for not significance. Mean with SEM is represented. t test was applied.

(0.06MB).

The possible correlation between the clinical–pathological characteristics of the patients and the presence or absence of CTCs (based on cytokeratin, CK, expression) was investigated before surgery (CTC1, Table 2) and one month after surgery (CTC2, Supplementary Table 3). The presence of CTC before surgery (CTC1) was significantly correlated with age, lymph node status (N), tumor stage, and adjuvant treatment. Regarding tumor stage, most patients with stages II and III were CTC− (p=0.0300). Focused on N status, most patients with lymph node metastasis were CTC− (p=0.0320). Consequently, it was observed that most CTC− patients received adjuvant treatment compared with the CTC+ patients (p=0.0040).

At CTC1, CTCs were detected in 26/76 patients (34.2%) with a median number of 3 in those CTCs+ patients (IQR=3, range=1–8) meanwhile, at CTC2 the presence of CTCs was detected in only 3/76 patients (3.9%) with a median number of 2 CTCs in those CTCs+ patients (IQR=3, range=1–4). In addition, the dynamics of CTCs between the two time points, CTC1 and CTC2, were also analyzed. We observed a significant decrease in the number of CTCs because of the surgery (p<0.0001, Wilcoxon test) (Fig. 4A).

Fig. 4.

Dynamics of circulating tumor cells (CTCs) and expression of EZH2 in CTCs from NSCLC patients. (A) Differences in the number of CTCs between CTC1 (before surgery) and CTC2 (one month after surgery). ****: p value for Wilcoxon test, p=0.0001. (B) Comparison of the percentage of expression of EZH2 at CTC1 and CTC2. Abbreviation n.s. is not significant for Mann–Whitney U test, p=0.2562).

(0.07MB).

The expression of EZH2 on CTCs was measured (see CTC characterization by immunofluorescence methods Section) at CTC1 in the 26 CTC+ patients. We found that 13/26 patients were EZH2+ (50.0%), with a median percentage of expression of 6.25% (IQR: 54.17%, range: 0–100%). At CTC2, EZH2 was detected in 2/3 CTC+ patients (66.67%), with a median percentage of expression of EZH2 of 100% (IQR: 100%, range: 0–100%). Expression of EZH2 showed no significant differences between CTC1 and CTC2 (Fig. 4B).

The prognostic role of circulating tumor cells (CTCs) for PFS and OS was evaluated using either the presence/absence of CTCs or a cut-off of p=0.5430 and p=0.9070, respectively; see Supplementary Fig. 3). Using the CTC cut-off of p=0.1920 and p=0.4380, respectively). OS was measured at 5 years since surgical treatment.

Interestingly, when evaluating EZH2 expression in the CTCs using a cut-off of <70% or ≥70% for the previously calculated percentage of expression (EZH2low/EZH2high, respectively), significant differences in terms of OS were observed at CTC1 (p=0.0180, Log-Rank, Mantel–Cox) with EZH2high patients having an 87% greater risk of death. In terms of PFS no significant differences were found for EZH2 characterization (p=0.9130) (Fig. 5C and D).

Fig. 5.

Prognostic role of CTCs and EZH2 expression on CTCs. Kaplan–Meier at CTC1 of the effect of the number of CTCs based on a cut-off (high≥2/low<2) in the relapse-free survival (A) and overall survival (B). Kaplan–Meier of the influence of the EZH2 expression cut-off in the progression-free survival (C) and overall survival (D). Log-Rank (Mantel–Cox) test was applied in all the panels.

(0.37MB).

At CTC2 neither presence/absence of CTCs nor CTC2low/CTC2high showed significant differences in terms of PFS or OS (Supplementary Fig. 4). At this point, Kaplan–Meier curves could not be calculated for EZH2low/EZH2high in CTCs+ patients due to the very small sample size.

Presence or absence of CTCs at baseline was analyzed based on the tissue molecular profiles. While we did not find statistical differences between the presence or absence of CTCs and VUS in the studied genes, we did observe a noteworthy trend. There were two out of three patients without CTCs in the follow-up but with VUS in ADAM19 that locally recurred (p=0.1739). The only patient with CTCs (both at baseline and the follow-up) and VUS in ADAM19 did not recur at time of these analyses (Supplementary Table 4).