Two groups of participants were enrolled separately in this study. This included the first group of 13 advanced NSCLC subjects and 6 healthy controls for discovery purpose using transcriptomic microarray analysis. The major findings attained in the discovery group were attempted for published microarray datasets validation and for biological verification in the second independent group of 34 advanced NSCLC subjects and 14 healthy controls along with the analysis of protein expression by immunohistochemistry in resected lung specimens and immunoblotting in cell lysates of BAL (Figure 1). The study was approved by the Ethics Committees of Chang Gung Memorial Hospital, Linkou, and all the participants gave written and signed informed consent.

Bronchoalveolar lavage specimens were collected during routine diagnostic bronchoscopy and BAL was performed from the lung segment in which tumor was located and from right middle lobe or lingular lobe in the normal controls (NC) via instillation of 250 c.c. of saline fluid to harvest the cells. The BAL samples were immediately preserved at 4°C, and the cells were pelleted within 1.5 h via 300 × g centrifugation for 15 min. An aliquot of cells from each sample was separated for the determination of differential cell counts using a Cytospin™ 4 Cytocentrifuge (Thermo Fisher Scientific, Taiwan). The rest of the cells were immediately lysed in buffer RLT and preserved at −80°C for the later RNA extraction.

Total RNA was isolated from each of the BAL samples using RNeasy mini kits (QIAGEN, Valencia, USA), and the RNA purity and integrity were measured using a Nanodrop 1000 (Thermo Fisher Scientific, Taiwan) and Agilent Bioanalyser (Agilent, Santa Clara, CA, USA), respectively. All the RNA samples met the quality criteria of OD₂₆₀/OD₂₈₀ >2.0 and RNA integrity number >9.0. The extracted RNA was labeled with streptavidin-phycoerythrin conjugate and hybridized to an Affymetrix HG-U133 Plus 2.0 microarray (Affymetrix, Santa Clara, CA, USA) as recommended by the manufacturers. The raw data were quality assessed and preprocessed by robust multi-array average normalization using the Bioconductor R package affy. The gene expression levels were further adjusted by removing the unwanted effects from technical batches, sex, age, and smoking using the Bioconductor R package Surrogate Variable Analysis (sva).

Total RNA (1 µg) was utilized for reverse transcription using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Taiwan). Quantitative real-time PCR was performed using the Rotor-Gene 6000 Real-Time PCR System (Corbett-Research, Mortlake, Australia). Each 20 µl reaction contained 5 ng of cDNA, 10 µl of SsoFast EvaGreen Supermix (Bio-Rad, Taiwan), and a final concentration of 300 nM of forward and reverse primer. A set of 14 reported endogenous reference genes were examined for their expression levels in microarray data (Figure S1 in Supplementary Material). Following this approach, HPRT1 was determined as the endogenous reference due to its highly constitutive expression in the samples from both the cancer and control subjects. The primer pairs of target genes included IGJ, IGKC, CEACAM6, MMP7, SPP1, CXCL13, SLC40A1, CPA3, and YES1, which were designed and examined for specificity using Primer-BLAST (Table S1 in Supplementary Material). The amplification efficiency of each primer pair was determined using the standard curve method by serial template dilution to six orders of magnitude. All the reactions were performed in duplicates and the double-delta cycle threshold (ΔΔCt) method was applied for relative quantitation in line with the minimum information for publication of quantitative real-time PCR experiments (MIQE) guidelines (28).

Samples of formalin-fixed paraformaldehyde-embedded sections (4 µm) were dewaxed and rehydrated prior to heat-induced epitope retrieval using Tris-EDTA Buffer at pH 9.0. Primary antibodies of polyclonal rabbit anti-human kappa light chains/IGKC antibody (A0191, DAKO), polyclonal rabbit anti-human Immunoglobulin J chains/IGJ antibody (Atlas Antibodies, HPA044132, SIGMA-ALDRICH), and polyclonal rabbit anti-human CPA3 antibody (Atlas Antibodies, HPA008689, SIGMA-ALDRICH) were used. Visualization was performed using the Envision Mouse/Rabbit-HRP (DAKO) system and chromogen detection using 3-amino-9-ethylcarbazole (AEC).

Aliquots of total protein (30 µg) were loaded with Laemmli sample buffer for electrophoresis using polyacrylamide gel with a 4–20% gradient, and polyvinylidene difluoride membranes were used for protein transfer. The membranes were blocked with protein-free SuperBlock blocking buffer (GeneStar, Taiwan) before overnight incubation with polyclonal rabbit anti-human kappa light chains/IGKC antibody at a dilution factor recommended by the manufacturers. Blots were developed using peroxidase substrate with enhanced chemiluminescence (Pierce™ ECL, Thermo Scientific) and later exposed to a UVP BioSpectrum 810 Imaging System™ (Cambridge, UK).

All categorical variables were analyzed using Fisher’s exact test, while Student’s t-test was used for comparisons of means between the two groups. A linear model (Bioconductor R package limma) with adjustments for technical batch, age, sex, and smoking was used for differentially expressed genes (DEGs) analysis, and an original p-value <0.05 was considered statistically significant in the discovery phase of the study. A predictive model was generated by the nearest shrunken centroid learning algorithm using the R package pamr as previously described (29, 30). This approach coupled with the 10-fold cross validation method generated a reduced set of genes/biomarkers determined by their accuracy in classifying the phenotype of interest.

A total 47 advanced NSCLC patients and 20 NC were enrolled, and these study participants were assigned to either the discovery (n = 19) group or the validation (n = 48) group, with the two groups being balanced in terms of clinical and pathologic characteristics (Table 1). The BAL specimens subjected to cytospin preparation were routinely examined using standard morphology criteria, where similar cell distributions of macrophages (85%), lymphocytes (10–15%), and granulocytes (<5%) were identified for all the study subjects.

Principle component analysis (PCA) was used to visualize the processed microarray data, in which removal of the effects from technical batches was confirmed (Figure S2 in Supplementary Material), leaving cancer as the major source of variation for the whole expression matrix (Figure 2B). In our attempt to explore the characteristic transcriptomic profile of BAL cells from tumor-bearing lung segments, we first analyzed the DEGs of BAL cells between the advanced NSCLC patients and the NC in the discovery group, and a total 129 DEGs (log₂ fold change >0.8 and nominal p < 0.05) were identified (Figure 2A). As a next step, we endeavored to understand their biological relevance. To this end, gene ontology enrichment analysis (GOEA) was performed on three publicly available ontology databases (GO, Reactome, and KEGG) in order to identify the links between the DEGs and specific functional categories. We combined network visualization for the enriched GO-terms using Cytoscape Enrichment Map plug-in (31) and found that the topology of gene enrichment was concentrated in a handful of functional domains, including Fc-receptor activation; the positive regulation of lymphocyte activation; humoral immunity and B cell activation; and extracellular matrix, cell adhesion, and proteinase activity (Figure 2C). These functional domains were represented by multiple pathways involving Fc gamma receptor-dependent phagocytosis, circulating immunoglobulin complex, endopeptidase activity, and extracellular matrix interactions (Table 2).

To ascertain whether the DEGs could be credited to specific cell lineages, we subsequently carried out a number of in silico analyses on the published microarray datasets from the Gene Expression Omnibus. The original expression value of each dataset was transformed by standardization and mean-centered to enable comparative analysis (32). Given that the identified DEGs were from BAL cells of tumor-bearing lung segments, it was essential to clarify that the signals were not dominated by potentially contaminated cancer cells. To address this, these DEGs were extracted from microarray data of human resected lung cancer (33) tissue (which consisted of a mixture of malignant, matrix, and infiltrating immune cells, n = 58) and of 9 lung cancer cell lines of the NCI-60 human tumor cell lines panel (34) (which had high purity of malignant cells, n = 26) for PCA and clustering analysis. The results showed that the identified DEGs clearly separate the two datasets along the first principle component (Figure 3A) and also unambiguously assign them to different clusters (Figure 3B), suggesting that the source of the DEGs was most likely non-malignant cells. This finding was reinforced as a set of 14 genes from The Cancer Genome Atlas that specifically characterize malignant cells (35) prevented the discrimination between the resected tumors and the cell lines (Figure S3 in Supplementary Material).

Later, the DEGs were taken to interrogate the normal human tissue gene expression database from Su et al. (36). A tissue is defined as overexpressing a gene when its expression level of the gene is twice the SD above the mean for all tissue types. We found that the identified DEGs were significantly expressed in hematopoietic cell lineages, particularly CD19 B cells, CD33 myeloid cells, BDCA4 dendritic cells, and CD14 monocytes (Figure S4 in Supplementary Material). Hierarchical clustering on the other two datasets of fractionated blood (37) and BAL cells (38) of healthy controls also suggested alveolar macrophages, B lymphocytes, and granulocytes as the major sources of these DEGs (Figure S5 in Supplementary Material).

As the characteristic transcriptomic profile of BAL cells was noted in the advanced NSCLC patients, we endeavored to understand whether that profile could constitute a common feature that has happened earlier at the early stage of lung cancers. To this end, the expression values of the identified DEGs were extracted from eight published microarray datasets of early stage resected tumors (n = 273) and tumor adjacent normal lung tissues (n = 247) (39–42), as well as normal non-diseased lung tissues (n = 87) (43–46). The original expression value of each dataset was preprocessed by standardization and mean-centered to enable comparative analysis (32). PCA of the expression matrix showed three separated groups that clearly matched the various tissue types (Figure 4A) and unbiased clustering that also differentiated these tissue types from each other (Figure 4B), suggesting what appears histologically normal tumor adjacent lungs are otherwise immunologically different from non-diseased lungs and its paired tumor tissues. As a further step, using Gene Set Variation Analysis (GSVA) as previously described (29), we found that the gene set from the composite of immunoglobulin genes (Figure 1A) showed a gradient increase of expression all the way from the tissues of non-diseased lungs, to tumor adjacent lungs, to the tumors themselves (Tukey HSD test: p < 2.0 × 10⁻⁷ for each pair of comparisons, Figure 4C). In addition, a significantly high expression of the gene set of mast cell proteinase CPA3, TPSAB1, and TPSB2 was also noted in the tumor adjacent lungs (Tukey HSD test: p < 2.0 × 10⁻⁷ for each pair of comparisons, Figure 4D).

To determine the signature that best characterized the transcriptomic profile, the shrunken centroid algorithm coupled with 10-fold cross validation was applied as previously described (29, 30) in the discovery group. A model of the signature reduced to 62 genes (Table S2 in Supplementary Material) was identified with the classification accuracy at 83% (Figure 5A). Overall, 53 (85.5%) of those genes were found to be upregulated in NSCLC, and their global expression across all the subjects in the discovery group was displayed (Figure 5C). Interestingly, these 53 upregulated genes, summarized by the GSVA as a signature of BAL cells from tumor-bearing lung segments, showed moderately high correlations with the gene expressions of multiple checkpoint proteins (Figure 5B), including PD-1/PDCD1 (r = 0.59, p = 0.0078), CTLA4 (r = 0.62, p = 0.0046), TIM-3/HAVCR2 (r = 0.56, p = 0.0130), LAG3 (r = 0.51, p = 0.0260), IDO1 (r = 0.43, p = 0.0660), and TIGIT (r = 0.52, p = 0.0220), where similar correlations, except for that with CTLA-4, were also noted in the published microarray datasets (Figure S6A in Supplementary Material). Furthermore, PCA of this signature presented a dichotomous gene expression pattern of these BAL cells from NSCLC (Figure S6B in Supplementary Material), indicating the existence of inter-patient heterogeneity.

As a next step, validation of a handful of top-ranked and intriguing biomarkers was attempted. This included SPP1, which plays a role in cell–matrix interactions and cytokine activity; CEACAM6, a glycoprotein involved in cell adhesion and migration; MMP7, which encodes a member of the metalloproteinases associated with matrix breakdown; SLC40A1, which encodes a membrane protein functioning in iron metabolism; IGJ, which is responsible for an immunoglobulin J chain polypeptide; IGKC, which is responsible for a constant region in immunoglobulin kappa light chains; CPA3, which encodes carboxypeptidase A3 frequently released by mast cells; YES1, which is responsible for a Src family protein tyrosine kinase; and CXCL13, which encodes a B lymphocyte chemoattractant. Using HPRT1 as an endogenous reference, RT-qPCR confirmed the overexpression of these genes in a manner consistent with the pattern found in the microarray study of the discovery group (Figure S7 in Supplementary Material). As an extension of this validation, we subsequently assessed the reproducibility of the expression pattern in an independent group of 34 NSCLCs and 14 NC, in whom we confirmed that the pattern of gene expression found in the discovery group was recapitulated (Figures 6A–C). These nine genes were later utilized in predictive model training using support vector machine in the discovery group, and the high performance of the resulting model was ascertained by ROC curve analysis (AUC: 0.920, 95% CI: 0.831–0.985, p = 2.98e–07; Figure 6D) when applied in the validation group.

Immunohistochemistry was later carried out for the validation of protein expression in early stage resected NSCLCs focusing on the immunoglobulins and mast cell carboxypeptidase A3, where tissues of tumor adjacent normal lungs and non-diseased lungs from surgical samples of pneumothorax patients being used. Significantly higher staining levels of the proteins IGKC (Wilcoxon signed-rank test, Figures 7D–F,K; p < 0.001), IGJ (Figures 7G–I,L; p < 0.01), and CPA3 (Figures 7A–C,J; p < 0.001) were noted in the tumor adjacent normal lung tissues. The increased staining levels of these proteins were more frequently found in the matrix tissue than in the alveolar air space of the lungs, particularly in the case of IGKC protein (Figure 7E). Immunoblotting also revealed that the lysates of BAL cells of NSCLC presented a significantly higher protein level of IGKC (Wilcoxon signed-rank test, p < 0.01), for both the full-length (25 kD) and the cleaved forms, than the lysates of BAL cells of NC (Figure 8).