Patients and sample collection

Samples were selected from the NYU Thoracic Surgery Archives (NTSA). Established in 2006, the NTSA has prospectively collected serum, plasma, buffy coat, and peripheral blood mononuclear cells, along with lung cancer and matching normal lung specimens under the IRB-approved 8896 protocol. Patients identified on preoperative workup as having a pulmonary nodule suspicious for lung cancer were consented for the collection of blood and snap frozen tissues (tumor and remote lung from the same lobe/segment) in the operating room at the time of their resection. Lung and matching tumor are sterilely cut at the operating room table, transferred to pre-labeled NuncTM vials and immediately snap frozen in liquid nitrogen within 10 min of resection. Samples are de-identified for storage at −80°C until use.

Demographic, clinical, and pathological data are recorded in an encrypted Research Electronic Data Capture (REDCap) database. Patients are seen at 3-month intervals for 2 years, then at 6-month intervals for 1 year, and then annually, with CT scans performed for surveillance in order to document any systemic and loco-regional recurrences, or the development of a second primary tumor. For this study, 48 patients with stage II NSCLC and >3 tumor and normal lung sample aliquots available in the NTSA were selected for microbiome measurement. We excluded patients with <1-month survival after resection, resulting in inclusion of 46 patients (Supplementary Figure 1). For microbiome analysis, we excluded samples with low microbial sequence reads (explained in further detail below), resulting in 39 tumor samples and 41 normal lung samples (34 patients with both tumor and normal lung samples) (Supplementary Figure 1). For peripheral gene expression analysis, we excluded 3 patients missing a buffy coat sample and 1 additional patient with low gene expression counts, leaving 42 patient samples (Supplementary Figure 1).

Outcome definitions

Endpoints were defined according to the consensus agreement in Punt et al. [29]. Disease-free survival (DFS) includes recurrences (loco-regional and systemic), new primaries (same or other cancer), and death from any cause as events. Recurrence-free survival (RFS) includes recurrences (loco-regional and systemic) and death from any cause as events, ignoring new primaries as events. Overall survival (OS) includes death from any cause as events. For all endpoints, person-time is defined as time from surgery to event or loss to follow-up (censored). For two patients missing the date of death, date of recurrence was assigned as the date of death; sensitivity analyses excluding these patients did not substantially change our findings.

Covariates

Data available for these patients were age (years), sex (male, female), race (white, non-white), smoking status (never, former, current), histology (adenocarcinoma, squamous cell carcinoma), stage (IIA, IIB), lymphovascular invasion (yes, no), pleural invasion (yes, no), tumor size (cm), positive nodes (%), and chemotherapy (yes, no). Variables associated with RFS, DFS, or OS in unadjusted or multivariable Cox proportional hazards models (p<0.10) were included in survival models of microbiome-related parameters. These variables were age, sex, race, smoking status, histology, and chemotherapy (Supplementary Table 1).

16S rRNA gene sequencingAssay

DNA was extracted with the DNeasy PowerLyzer PowerSoil DNA Isolation Kit (QIAGEN, Valencia, CA), following the manufacturer’s instructions. PCR amplification was performed on the 16S rRNA gene V4 hypervariable region using the 515F and 806R primers, with a 12-bp unique Golay barcoding [30]. PCR reactions were performed with an initial denaturation of 95 °C for 5 min, followed by 15 cycles of 95 °C for 1 min, 55 °C for 1 min, and 68 °C for 1 min, followed by 15 cycles of 95 °C for 1 min, 60 °C for 1 min, and 68 °C for 1 min, and a final extension for 10 min at 68 °C on a GeneAmp PCR System 9700 (Applied Biosystems, Foster City, CA). The PCR products were purified using a QIAquick Gel Extraction Kit (QIAGEN, Valencia, CA) and quantified using a Qubit 2.0 Fluorometric High Sensitivity dsDNA Assay (Life Technologies, Carlsbad, CA). PCR products were pooled relative to their band intensity. KAPA LTP Library Preparation Kit (KAPA Biosystems, Wilmington, MA) was used on the combined purified PCR products according to the manufacturer’s protocol and the size integrity of the amplicons with Illumina indices was validated with a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) at the Genomics Core at Albert Einstein College of Medicine. High-throughput amplicon sequencing was conducted on a HiSeq (Illumina, San Diego, CA) using 2 × 150 paired-end fragments.

Sequence read processing

Sequence reads were processed using QIIME 2 [31]. Briefly, sequence reads were demultiplexed, followed by quality filtering as described in Bokulich et al. [32]. Next, on the forward reads, the Deblur workflow was applied, which uses sequence error profiles to obtain putative error-free sequences, referred to as “sub” operational taxonomic units or amplicon sequence variants (ASVs) [33]. ASVs were assigned taxonomy using a naïve Bayes classifier pre-trained on the Greengenes [34] 13_8 99% OTUs, where the sequences have been trimmed to only include 150 bases from the 16S V4 region, bound by the 515F/806R primer pair. PICRUSt2 was used to predict abundance of functional pathways [35].

Quality control

Samples were processed in two batches, each including a negative PCR control and a positive PCR control (ZymoBIOMICS Microbial Community Standard, Zymo Research, Irvine, CA), and all sequenced together in a single run. Based on examination of the number of sequence reads in the tumor and normal lung samples, negative controls, and positive controls, we excluded samples with low sequence reads similar to the negative controls (<6000 reads/sample) (Supplementary Figure 2). After this exclusion, the median sequence read count per sample after the Deblur workflow was 54,864 (Q1= 33,598, Q3= 74,935). To remove potential contaminant sequences, we used two approaches as outlined in Cao et al. [36]: rare taxa filtering implemented using “PERFect” [37] and contaminant removal based on prevalence in negative controls and frequency among samples implemented using “decontam” [38]. Use of these complementary approaches in conjunction should remove contaminant taxa that are rare or abundant, respectively, and result in data with lower dimensionality, minimal information loss, reduced technical variability, and increased capacity for reproducibility and comparability across studies [36]. We performed contaminant filtering at each taxonomic level (i.e., phylum, class, order, etc.) separately, and then additionally removed contaminants at lower taxonomic levels that were identified as contaminants at higher taxonomic levels (e.g., removing taxonomic classes belonging to a contaminant phylum). After contaminant removal, the median sequence read count per sample was 16,418 (Q1=8,143, Q3=28,334). Contaminant-filtered data was used for all downstream processes, including running of PICRUSt2, calculation of α-diversity and β-diversity, and all statistical analyses described below.

The lowest sequencing depth among the samples, 2468 after contaminant removal, was sufficient to characterize the ranking of diversity among the samples, though not the ranking of number of observed ASVs (i.e., richness) (Supplementary Figure 3). Genus-level composition of the positive controls indicated appropriate detection of expected genera (Supplementary Figure 4). Among 4 samples run in duplicate, only 1 sample achieved sufficient sequencing depth in both duplicates (>2468 reads after contaminant removal); for this sample, compositional reproducibility was excellent, as was reproducibility for the two positive controls, based on similarity in principal coordinate analysis (Supplementary Figure 5). Additionally, we did not observe any compositional batch effect among the samples (Supplementary Figure 5).

Quantitative PCR (qPCR)

A previously identified and published set of Clostridia specific primers (SJ-F/SJ-R) were selected [39]. Though cycling conditions were outlined in the initial publication, a range of annealing temperatures and cycling conditions were assessed to optimize performance. Assessments, similar to eventual qPCRs, were carried out on subsets of sample DNA, a ZymoBIOMICS Microbial Community DNA standard (D6305), an extraneous sample of stool DNA, and finally, reference Clostridium difficile gDNA, strain 4206, acquired from ATCC (BAA-1872D-5). A 2x Applied Biosystems SYBR Green PowerUp master mix (25742) was used and qPCR carried out on the Applied Biosystems Viia7 machine in the Albert Einstein College of Medicine Genomics core. Master mix was prepared according to the manufacturer’s instructions, primer concentrations in recommended range at 500 nM, with only deviation being the addition of MgCl2 to 0.5 uM. All samples were run in triplicate on a single Applied Biosystems 384-well plate (4309849). Lung sample DNA concentrations were measured, then subsequently normalized, and loaded to plate for a total of 15 ng DNA/well. qPCR was set to 40 amplification cycles with an annealing step of 15s at 53 °C and an extension of 60s at 72 °C. All other conditions followed master mix instructions. To generate a standard curve, seven dilutions of C. difficile gDNA were used, ranging from 10ng/well to 10fg/well.

Gene expression

Gene expression was measured in buffy coat collected the day of and prior to tumor resection, using the nanoString® nCounter® PanCancer Immune Profiling Panel at the NYU Thoracic Oncology laboratory. Color-coded barcodes, attached to single target-specific probes corresponding to analytes of interest, hybridize directly to target molecules and are individually counted. The Panel provides multiplex gene expression analysis for 770 genes from 24 different immune cell types, common checkpoint inhibitors, and genes covering the adaptive and innate immune response. Gene expression data was normalized using housekeeping genes with the nSolver software.

Statistical analysisα-Diversity

α-Diversity (within-sample microbiome diversity) was assessed using richness (number of ASVs) and the Shannon diversity index, calculated using the QIIME 2 diversity plugin. The final values for each sample were calculated by averaging the richness and Shannon index from 100 iterations of rarefaction at 2468 sequence reads, the lowest sequence read depth among the samples after the exclusion of samples with low counts and contaminant filtering (described above).

β-Diversity

β-Diversity (between-sample microbiome diversity) was assessed using the Jensen-Shannon Divergence (JSD) [40]. Distances were calculated on the ASV level. Principal coordinate analysis (PCoA) [41] was used for visualization.

Comparison of the tumor and normal lung microbiome

Wilcoxon signed-rank tests were used to examine differences in α-diversity (number of ASVs, Shannon index) between paired tumor and normal lung samples. Permutational multivariate analysis of variance (PERMANOVA) [42] was used to examine whether overall bacterial composition (β-diversity) differed between paired tumor and normal lung samples, using patient ID as strata. PERMANOVA was also used to examine significant predictors of tumor and normal lung microbiome composition. To determine whether paired tumor-normal samples were more alike than unpaired samples, we used the JSD to calculate the average between-pair distance in overall bacterial composition for tumor and normal lung tissue sample pairs. We then calculated the average distances in 5000 permutations of random pairings of tumor and normal lung samples. If the mean pair distance was less than the 5th percentile of the permuted means, we concluded that the paired samples were more alike than random pairings of samples from different patients. For testing of differentially abundant taxa between tumor and normal lung samples, ASVs were agglomerated into phylum, class, order, family, genus, and species levels, in order to perform testing at each taxonomic level. Taxonomic abundance was transformed using the centered log-ratio (clr) transformation [43, 44] after adding a pseudocount, in order to remove compositional constraints of sequencing. Wilcoxon signed-rank tests were used to examine differences in taxonomic abundance between paired tumor and normal lung samples. P-values for taxa were adjusted for the false discovery rate (FDR) [45]; FDR adjustment was done at each taxonomic level (i.e., phylum, genus) separately.

Lung microbiome and survival

We used Cox proportional hazards models to determine whether α-diversity (number of ASVs, Shannon index) in the tumor or normal lung was associated with RFS, DFS, or OS, adjusting for age, sex, race, smoking status, histology, and chemotherapy. The community-level test of association between the microbiota and survival times (MiRKAT-S) [46] was used to test the association of overall bacterial composition (β-diversity), as measured by the JSD, in tumor or normal lung with RFS, DFS, and OS, adjusting for aforementioned covariates. Taxa and functional pathways associated with RFS, DFS, or OS were assessed independently in the tumor and normal lung samples, using repeated cross-validated elastic-net penalized Cox proportional hazards regression, as previously described [47]. For taxonomic analysis, ASVs were agglomerated into phylum, class, order, family, genus, and species levels, in order to perform testing at each taxonomic level. Taxonomic and functional pathway abundance was transformed using the centered log-ratio (clr) transformation [43, 44]. Agglomerated taxa missing taxonomic classification at the respective taxonomic level were removed (e.g., genera missing genus level classification), resulting in the inclusion of 6 phyla, 10 classes, 19 orders, 27 families, 41 genera, 13 species, and 787 ASVs. There were 359 functional pathways. We conducted 500 × 10-fold cross-validated elastic-net penalized Cox regression using the “cv.glmnet” function in the “glmnet” R package [48], with an α value of 0.5 to allow groups of correlated predictors to be selected together. Non-penalized covariates (age, sex, race, smoking status, histology, and chemotherapy) were included in each model. We summed the number of times each taxon or functional pathway was selected out of the 500 repetitions. For all tested taxa and functional pathways, we also fit standard Cox proportional hazards models for RFS, DFS, and OS, adjusting for the covariates listed above. P-values for these models were adjusted for the FDR [45]; FDR adjustment was done at each taxonomic level (i.e., phylum, genus) separately. We focused further on taxa and functional pathways selected ≥25% of the 500 times (125 times or more) and with FDR-adjusted q < 0.20. We used Spearman’s correlation to examine associations between the relative abundance of taxa and functional pathways.

Peripheral gene expression and survival

Gene expression data was log2 transformed for analysis. As shown above for microbiome, we used 500 × 10-fold cross-validated elastic-net penalized Cox regression to identify genes for which expression was related to RFS, DFS, or OS. We focused further on genes selected ≥25% of the 500 times (125 times or more) and with FDR-adjusted q < 0.20. The “topGO” package in R was used to determine gene ontology (GO) enrichments for survival-related genes. We used Spearman’s correlation to examine associations between survival-related genes and survival-related taxa and functional pathways.

Survival risk prediction

We considered whether the identified lung microbiome and/or peripheral gene expression biomarkers may improve risk prediction for RFS, DFS, or OS over standard covariates (age, sex, race, smoking status, histology, and chemotherapy) alone. Annual time-dependent area under the receiver operating characteristic (ROC) curve (AUC), starting at 1 year after surgical resection, were estimated using the R package “timeROC”, and 95% confidence intervals were constructed from the distribution of AUCs from 1000 bootstraps of the data. Permutation tests were used to test whether the biomarkers provide significant additional survival risk discrimination over standard covariates. Briefly, we generated null distributions of difference in time-dependent AUCs, by generating 5000 permutations of the microbiome and/or gene data while not disrupting the link of the survival time/censoring with the standard covariates. If the test statistic (difference in AUCs between microbiome and/or gene model and standard covariate model) was greater than the 95th percentile of the null distribution at a given time, we deemed that the biomarkers added significant (p < 0.05) risk discrimination over the standard covariates at that time.

Validation cohorts

To validate observed findings of the current study, we included 16S microbiome sequencing data from two former studies: our previous pilot study of tumor and normal lung of NSCLC patients, referred to here as Peters et al. [21], and a previously described study of lower airway brushings of NSCLC patients, referred to here as Tsay et al. [17]. Sequencing data from Peters et al. [21] and Tsay et al. [17] was processed in QIIME2 as described above (see the “Sequence read processing” section). Contaminant filtering was performed at each taxonomic level using “PERFect” [37]. For Peters et al. [21], 16 tumor and 17 normal lung tissue samples were available, and covariates included in analysis were limited to age (years) and sex (male, female) due to low sample size. For Tsay et al. (2021), a total of 183 samples from 73 patients were available. Samples were categorized as “involved” (i.e., same side as tumor) or “non-involved”, and patients were categorized as having “local” (I-IIIA) or “advanced” (IIIB-IV) stage disease. Covariates included in the analysis were age (years), sex (male, female), race (white, non-white), smoking status (never, former, current), histology (adenocarcinoma, squamous cell carcinoma), chemotherapy (yes, no), and surgery (yes, no). To validate associations of selected taxa with survival, we fit Cox proportional hazards models for OS with clr-transformed taxa as predictors, adjusting for covariates (and accounting for within-patient clustering of samples in the case of Tsay et al.).