Transcriptomic Microenvironment of Lung Adenocarcinoma

The existence of tumor-associated molecular abnormalities in histologically normal tissue adjacent to lung tumors is well established (1, 2). Gene expression profiling of adenocarcinoma and adjacent nontumor lung tissues from the same patients have revealed important insights about the molecular transitions that occur between normal and tumor lung tissues (3, 4). However, it is unclear whether these molecular alterations are switched on and off or are showing a progressive gradient that correlates with the distance of the nontumor sites relative to the tumor, akin to the airway field of cancerization observed in bronchial airway epithelial cells from patients with lung cancer (5). Whole-genome gene expression profiling of multiple histologically normal lung tissues collected at standardized distances from the tumor has not been characterized yet.

In this study, we sampled nontumor lung parenchyma at standardized distances relative to the tumor and obtained a unique geographic lung tissue sampling collection from patients undergoing lobectomy for stage I lung adenocarcinoma. The transcriptome of the tumor and nontumor lung tissues was investigated. This design was selected to study the molecular transitions that are stable and transient in nontumor lung compared with the tumor. We hypothesized that genes can be classified as upregulated, downregulated, knocked-down, or induced in tumor compared with nontumor lung. We also hypothesized that some genes may show patterns of mRNA expression that correlate with the distance of the nontumor sites relative to the tumor. The outcomes of this study are important to elucidate the transcriptomic microenvironment of lung adenocarcinoma, which has the potential to reveal uncharacterized mechanisms driving lung adenocarcinoma as well as to identify stable and transient biomarkers in tumor or adjacent nontumor lung tissues to personalize management and treatment options.

Figure 1A provides a schematic representation of the sample collection. Lung samples were collected from 12 patients undergoing curative lobectomy for stage I adenocarcinoma including systematic nodal dissection at the Institut Universitaire de Cardiologie et de Pneumologie de Québec – Université Laval (IUCPQ-UL; Quebec City, Canada). Only white patients of self-reported European ancestry were selected. Tumor and adjacent nontumor lung samples were collected at the pathologic department under the supervision of experienced pulmonary pathologists (C. Couture, P. Joubert, S. Pagé, S. Trahan). After surgical removal, lung specimens were immediately processed for pathologic evaluation and grossing. A 1 cm³, tumor sample was harvested as well as the adjacent non-neoplastic pulmonary parenchyma samples (1 cm³) located at 0, 2, 4, and 6 cm from the tumor. Lung samples were snap-frozen in liquid nitrogen and stored at −80°C. The time from surgical removal to storage was less than 30 minutes. Details for tissue processing, staining, and RNA extraction are provided in Supplementary Data. The clinical characteristics of the 12 patients are indicated in Table 1. Patients provided written informed consent and the study was approved by the IUCPQ-UL ethics committee.

For each patient, whole-genome gene expression was evaluated in tumor and nontumor specimens harvested at 0, 2, 4, and 6 cm from the tumor using the Illumina HumanHT-12 v4 Expression BeadChip. The complete dataset was deposited in the NCBI repository and is accessible through GEO Series accession number GSE83213.

The microarray raw data were log2-transformed and quantile normalized using the lumi package in R (6). Standard quality controls (QC) were performed (7). The gene expression profile of one tumor sample clustering with nontumor lung samples was excluded from further analysis (Supplementary Fig. S1). The default detection P value (P < 0.01) in lumi was used to call probes present or absent in each sample. The primary step in all analyses was to filter-out probes absent in all samples. To do so, probes absent in tumors and all nontumor lung specimens (n = 4 or 5 specimens × 12 patients = 57) were excluded. A flow chart of the analysis strategy is illustrated in Fig. 2. 

Identify genes induced and knocked-down in tumors compared with nontumor sites. The absent and present calls for all samples based on the detection P value (P < 0.01) were considered for this analysis. Genes that were called present in at least 10 of 11 tumors that passed QC (more than 90%) and absent in at least 42 of 46 nontumor specimens that passed QC (more than 90%) were identified as gene induced in tumors. In contrast, genes that were called absent in at least 10 of 11 tumors and present in at least 42 of 46 nontumor specimens were identified as genes knocked-down in tumors.

Identify genes up- and downregulated in tumors compared with nontumor sites. Differential gene expression analysis was performed using the significant analysis of microarrays (SAM) method (8) with two class paired analyses and keeping the false discovery rate (FDR) < 5%. A fold-change threshold of 2.0 was used to identify the top genes differentially expressed across tissue sites. All pairwise comparisons between the tumor and nontumor lung specimens at different sites (i.e., 0, 2, 4, and 6 cm) relative to the tumor were performed. Genes that were significantly modulated in at least one comparison were considered up- or downregulated in tumor.

Identify genes that are stably and transiently expressed in nontumor lung tissues (n = 3 or 4 specimens × 12 patients = 46). Mean mRNA levels were compared among nontumor tissues collected at 0, 2, 4, and 6 cm using paired t tests. Probes showing no statistically significant difference between any pairwise sites (P > 0.05) were considered stable. All other probes were considered transient.

Identify genes progressively increased and decreased in nontumor lung relative to the distance from the tumor. First, significant probes induced, knocked-down, up- and downregulated in Analysis 1 and 2 were combined to find a set of probes significantly modulated in the tumor compared with the nontumor lung. These probes were overlapped with transient probes identified in Analysis 3. The gene expression patterns of overlapping probes were plotted with respect to the distance relative to the tumor to identify transient genes concordantly differentially expressed in tumor. Second, formal statistical tests were then performed to identify probes that decreased or increased in nontumor tissues by the distance from the tumor. Three models were evaluated including a linear, ordinal, and polynomial models. The linear model assumed a progressive increase or decrease in gene expression from 0 to 6 cm away from the tumor. The ordinal model evaluated whether changes in gene expression followed a gradient by locations classified as categorical ordered labels. Finally, a second-degree polynomial model was used to evaluate whether the observed gene expression changes fitted a curve. In all three models, P values were generated for each probe and adjusted for multiple testing using Benjamini–Hochberg correction.

Because of the uniqueness of our sample collection, external validation was only possible for the comparison tumor versus adjacent nontumor lung tissues (Analysis 2). Two publicly available datasets with whole-genome gene expression data of lung adenocarcinoma and matched nontumor lung specimens were used for this purpose (3, 4). Methods are described in Supplementary Data.

Biological pathway analysis was performed with the Ingenuity Pathway Analysis system (Supplementary Data).

Table 1 shows the clinical characteristics of the 12 patients undergoing curative lobectomy for stage IA lung adenocarcinoma. Lung samples were obtained from 6 men and 6 women with a mean age = 65.2 ± 9.2 years. All were white patients of self-reported European ancestry. Most were former or current smokers with mean pack-years of 27 ± 11, suggesting a relatively high smoking history.

Complete geographic tissue sampling was obtained for the 12 patients (Fig. 1A). Figure 1B shows histologic slides of one representative patient including the tumor specimen and nontumor lung samples collected at 0, 2, 4, and 6 cm. Histologic slides of all patients are provided in Supplementary Fig. S2. Nontumor lung samples were free of tumor cells except for two patients at 0 cm (patients 2 and 7, see Supplementary Fig. S2). These two samples were excluded from further analysis. Each specimen was interrogated using whole-genome gene expression microarray testing a total of 47,323 probes. Probes called absent in all samples (n = 4 or 5 specimens × 12 patients = 57) were excluded from further analysis (n = 17,111) leaving 30,212 probes (21,292 genes) for subsequent analyses.

The goal of this analysis was to identify genes induced and knocked-down in tumors compared with non-tumor sites. All probes for each sample were dichotomized into an absent or present call. The distributions of absent probes for all samples and across all sample sites (i.e., tumor, 0, 2, 4, and 6 cm) are illustrated in Supplementary Fig. S3. On the basis of these distributions, induced probes were defined as being expressed in at least 10 of 11 tumors and nonexpressed in at least 42 of 46 nontumor samples. Using this definition, 31 induced probes (29 genes) were found (Supplementary Table S1). Analogously, knocked-down genes were defined as being absent in at least 10 of 11 tumors and present in at least 42 of the 46 nontumor samples. A total of 39 probes (34 genes) were identified as knocked-down in tumors (Supplementary Table S2). The expression patterns of induced and knocked-down probes with respect to tumor proximity are illustrated in Fig. 3 (see also heatmap in Supplementary Fig. S4).

The objective of the second analysis was to identify genes up- and downregulated in tumor compared with nontumor sites. These analyses were performed with 11 patients as the tumor site failed QC for patient #4. The number of probes differentially expressed (FDR ≤ 0.05 and fold change ≥ 2) between the tumor and the nontumor sites varied from 1,225 to 1,311. Figure 4 illustrates the magnitude of changes across the four comparisons. Absolute fold changes observed for significant genes varied from 2 to 31. Overlapping differentially expressed genes between tumor and nontumor sites were then identified. A total of 365 and 663 probes were found consistently up- and downregulated between tumor and nontumor sites (Supplementary Fig. S5). This corresponds to 312 up- and 566 downregulated genes. Considering all pairwise comparisons between tumor and nontumor sites, 482 genes (563 probes) were upregulated and 843 genes (998 probes) were downregulated in at least one comparison. Up- and downregulated genes as well as fold changes across the four comparisons are listed in Supplementary Tables S3-S4, respectively.

The goal of this analysis was to identify genes with stable expression in nontumor lung tissues, which has important implication for biomarker discovery. Mean expression values of all probes were compared across tissue sites. Comparisons with 0 cm (i.e., 0 vs. 2 cm, 0 vs. 4 cm, and 0 vs. 6 cm) were performed with 10 patients and other pairwise comparisons (i.e., 2 vs. 4 cm, 2 vs. 6 cm, and 4 vs. 6 cm) were performed with 12 patients. Probes were considered stable if they were not significantly different (P > 0.05) between any two sites. This liberal P value threshold was used without correction for multiple testing to have a robust list of genes that did not differ in expression in nonmalignant lung tissues relative to the distance from the tumor. A total of 15,700 unique genes (20,666 probes), were identified using these criteria and classified as stable (Supplementary Table S5). By default, probes not in the stable category were classified as transient (9,546 probes or 7,966 genes, Supplementary Table S6). Please note that transient does not mean statistically significant.

The objective of this analysis was to identify changes in gene expression occurring relative to the distance from the tumor. Probes significantly modulated in tumor compared with nontumor lung tissues were first evaluated. Probes induced (n = 31), knocked-down (n = 39), up- (n = 563), and downregulated (n = 998) identified in Analysis 1 and 2 were combined. These probes were overlapped with transient probes from Analysis 3, excluding stable probes that are clearly not influenced by the distance from the tumor. Figure 5 illustrates 7 induced, 12 knocked-down, 182 up-, and 376 downregulated probes that were considered transient in nontumor lung tissues. In all cases, the main differences in gene expression levels occurred between the tumor and nontumor sites. Changes across nontumor sites were relatively modest and not concordantly modulated with changes observed at the tumor site. More formal statistical tests were then performed to evaluate whether gene expression changes in nontumor sites were significant. These analyses were performed with 10 patients that passed QC for all sites (i.e., 0, 2, 4, and 6 cm). The P values distribution of all probes (n = 30,212) for the linear, ordinal, and polynomial models are illustrated in Supplementary Fig. S6. At a liberal P value of 0.05, 835, 863, and 615 probes were identified with the linear, ordinal, and polynomial models, respectively. However, none of these probes were significant after adjustment for multiple testing (Benjamini–Hochberg correction). Note that these models identified a large number of significant probes when the tumor site was included as part of the analysis (Supplementary Fig. S7). As expected, the best fit was with the polynomial model consistent with large differences observed in gene expression levels between tumor and the nontumor site at 0 cm and the negligible/stochastic changes observed across nontumor tissue sites.

Genes differentially expressed in tumor compared with nontumor tissues identified in this study (Analysis 2) were compared with two publicly available whole-genome gene expression datasets of adenocarcinoma tumors with adjacent nontumor lung tissues. To ensure proper comparison with this study, raw data were obtained from each dataset and analyzed using the same methods for Illumina arrays and equivalent methods for Affymetrix arrays (see Supplementary Data). Overall, 98 up- and 251 downregulated genes overlapped across the three datasets (Supplementary Table S7). A total of 288 of 482 (59.8%) upregulated genes and 579 of 843 (68.7%) downregulated genes in this study were validated in at least one external dataset (Supplementary Fig. S8).

Four gene lists were investigated for biological pathways using IPA including genes induced and knocked-down in tumors (Analysis 1) and genes up- and downregulated in tumors (Analysis 2). The top canonical pathways identified for each gene list are indicated in Supplementary Table S8. Six and 15 pathways were significant for up- and downregulated genes, respectively (adjusted P < 0.05). The top pathway for upregulated genes was the Inhibition of matrix metalloproteases (P = 3.98E−08). The top pathway for downregulated genes was the hepatic fibrosis/hepatic stellate cell activation (P = 5.37E−10). Interestingly, two related pathways were significant for both up- and downregulated genes namely the granulocyte adhesion and diapedesis and the agranulocyte adhesion and diapedesis. A total of 14 and 25 genes in these two pathways were up- and downregulated in tumors, respectively (Supplementary Fig. S9).

This study comprehensively evaluated the transcriptomic microenvironment of stage I pulmonary adenocarcinoma. Briefly, 29 genes were found induced, 34 knocked-down, 482 up-, and 843 downregulated in tumor compared with nontumor lung tissues. These genes are important to understand the molecular transition that occurred from normal lung to adenocarcinoma. The geographic lung tissue sampling of nontumor lung tissues collected at standardized distances relative to the tumor is unique to investigate changes in gene expression with respect to the proximity of the tumor. Interestingly, we demonstrated that the magnitude of gene expression changes between tumor and nontumor sites was similar regardless of the distances with the tumor. At the gene-specific level, we found no gene that were statistically increased or decreased in nontumor lung tissues with the distance from the tumor. Overall, the results are consistent with distinct gene expression profiles for adenocarcinomas and adjacent nontumor lung specimens, but not with significant change in nontumor sites relative to the distance from the tumor.

Genes found differentially expressed in adenocarcinoma compared with normal lung are important to reveal the molecular oncogenesis of the disease. Previous experimental studies have demonstrated the biological relevance of some genes identified in this study. LGR4 (leucine rich repeat containing G protein-coupled receptor 4), also known as GPR48, was among the genes found consistently upregulated in adenocarcinoma. Upregulation of this gene is known to increase invasive activities of lung cancer cells and promote lung metastasis in vivo (9). Other disease drivers identified as upregulated include SPP1 (secreted phosphoprotein 1, also known as osteopontin) and LAMB3 (laminin subunit beta 3) that were shown upregulated at the protein levels in lung cancer tissues compared with adjacent normal lung and to contribute to lung cancer cell invasion and metastasis in vitro and in vivo (10). A number of matrix metalloproteinases were also upregulated in tumor including MMP1, 7, 9, 10, 11, and 12. Mice models demonstrated the role of MMP1 (11), MMP9 (12), and MMP12 (13) in lung cancer development. Similarly, some downregulated genes were shown to experimentally suppress tumor growth or play a protective role against lung tumorigenesis such as THBS1 (thrombospondin 1; ref. 14), SCGB1A1 (secretoglobin family 1A member 1, also known as CC10; ref. 15), IL33 (16), IL1β (17), and BMPR2 (bone morphogenetic protein receptor type 2; ref. 18). Together, these experimental evidences demonstrate that a number of identified genes are potential drivers of disease. Further functional studies are needed to test the other genes and delineate the temporal sequence of molecular events as well as interaction with driver mutations that underlie lung adenocarcinoma development.

The lack of genes progressively increased or decreased in nontumor lung with shorter distances from the tumor goes against our hypothesis and suggests little impact of the tumor on surrounding normal lung tissues. Paracrine effects and even mechanical strains caused by the expansion of the tumor would be expected to modulate gene expression of surrounding tissues. In addition, clinically relevant lung cancer gene expression signatures were identified in nontumor lung tissues (19) and airway epithelial cells (20), suggesting that tumor-associated gene expression alterations can be detected in the entire respiratory tract. Previous studies in human airways of cancer patients have also demonstrated site-dependent gene expression profiles (5, 21). More specifically, 422 gene features were identified as progressively differentially expressed in bronchial airway epithelial cell samples varying in distance relative to the tumor (5) supporting the concept of “Field of injury” or “Field of cancerization”. If the tumor is the cause of the concordant molecular alterations observed in surrounding tissues, we would expect these alterations to occur in adjacent normal lung and not only in airways. Considering the results of this study, the tumor is most likely to be the consequence of alterations caused by other factors, for example, smoking, that may have molecularly altered airways, but may not have reached surrounding adjacent normal lung in stage I lung cancer patients. However, cell type specificity and most particularly the percentage of replicating and nonreplicating lung epithelial cells between airways and whole lung may also explain the different results.

Tumors are heterogeneous tissues containing multiple subclones and varying levels of cancer and normal cells (22, 23), which are characteristics that make the development and validation of robust biomarkers challenging. Accordingly, investigators are progressively moving away from the tumor to develop lung cancer biomarkers (19, 24). For example, Frullanti and colleagues have identified a prognostic lung cancer gene expression signature in nontumor lung tissues (19). Residual nonmalignant lung tissues explanted during surgery represent a unique biospecimen media to develop prognostic factors. Adjacent nontumor lung specimens are more representative of tissues that remain in the patient after surgery and are thus more likely to contain the molecular signature that predicts relapse and death. An issue arising from testing biomarkers in nontumor lung samples is whether the exact location of the specimen assayed relative to the tumor has an impact. In this study, we have provided lists of genes with stable and transient expression levels in nontumor lung with respect to the distance from the tumor. These results represent an important step forward to develop biomarkers in residual nonmalignant lung tissues explanted during surgery. Stable biomarkers up- and downregulated in lung cancer patients may have great prognostic value and may also identify disease drivers.

Alterations driving the development of lung adenocarcinoma are still unclear. In this study, we found 482 genes upregulated and 843 genes downregulated in adenocarcinoma compared with adjacent nontumor lung. A subset of 312 up- and 566 downregulated genes were observed across all comparisons between tumor and nontumor sites. In addition, nearly 60% of upregulated genes and 70% of downregulated genes were validated in at least one external dataset of paired adenocarcinoma and normal lung tissues. This study is thus zooming in on a list of potential key triggering molecular alterations involved in the oncogenesis of the disease. We also observed a predominance of genes down- compared with upregulated in tumor in the three datasets. To summarize, we found 482 up- and 843 downregulated genes in this study, 583 up- and 790 downregulated genes in the dataset by Selamat and colleagues (3), and 223 up- and 430 downregulated genes in the dataset by Landi and colleagues (4; Supplementary Tables S9 and S10; Supplementary Fig. S10). In the three datasets, the number of downregulated genes nearly doubles the number of upregulated genes. These results are not sufficient to demonstrate whether these genes are a cause or consequence of adenocarcinoma development, but are consistent with a global downregulation of anticancer genes and/or pathways. In this study, we noted the overexpression of the oncogene ERBB3 (25) and downregulation of the tumor suppressor gene TGFBR2 (26). However, a number of oncogenes were also found downregulated in tumor including ROS1, KLF4, AFF3, PDGFB, and GATA2. Accordingly, the direction of effect was not necessarily consistent with upregulation of oncogenes and downregulation of tumor suppressor genes.

This study has limitations. First, sample size is not sufficient to stratify analyses by sex, smoking, KRAS mutation status, or other clinically relevant characteristics. Second, we have selected a homogeneous group of lung cancer patients, that is, white patients with stage I adenocarcinoma undergoing curative lobectomy. Accordingly, the results cannot be extrapolated to other histologic subtypes of lung cancer and requires a word of caution before generalizing the results to other ethnic groups and more advanced stages of adenocarcinoma. Third, we have only measured gene expression levels. Other molecular phenotypes including methylation marks, germline, and somatic genetic variants may vary with respect to the spatial proximity from the tumors and will need to be investigated in future studies. Fourth, we have measured gene expression using microarrays. RNA sequencing was a more costly option, but with the potential to provide a more comprehensive characterization of the transcriptomic microenvironment of pulmonary adenocarcinoma including differentially expressed isoforms. Fifth, gene expression profiles were derived from whole tissues containing different proportion of malignant cells in tumor samples (Table 1). Cellular heterogeneity is also expected for nontumor lung samples. Accordingly, our results must be interpreted with caution and further validation in larger sample sets and follow-up experiments using specific cell types will be required. Finally, the transcriptome of nontumor lung was obtained only from patients with lung cancer. A new study design is needed to evaluate how much of the normal-appearing tissue surrounding tumors differs from normal tissue collected from noncancer patients.

In conclusion, this study provides the first whole-genome gene expression profiles of pulmonary adenocarcinoma and adjacent multi-site nontumor lung tissues collected at standardized distances relative to the tumor. Results add to existing data showing the molecular transitions that occur from normal lung to adenocarcinoma and narrow down the number of genes consistently deregulated across nontumor tissue sites and datasets. This study also provides genes with stable mRNA levels in nontumor lung tissues with respect to the distance from the tumor, which are essential for biomarker discovery in this media. Finally, results demonstrate only subtle nonsignificant changes in gene expression across nontumor sites suggesting little effect of the tumor itself on surrounding tissues.

No potential conflicts of interest were disclosed.

Conception and design: Y. Bossé, P. Joubert

Development of methodology: Y. Bossé

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y. Bossé, O. Sazonova, N. Gaudreault, N. Bastien, S. Trahan, C. Couture, P. Joubert

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y. Bossé

Writing, review, and/or revision of the manuscript: Y. Bossé, O. Sazonova, N. Gaudreault, N. Bastien, M. Conti, P. Joubert

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): O. Sazonova, S. Pagé, P. Joubert

Study supervision: Y. Bossé

The authors would like to thank the team at the IUCPQ site of the Respiratory Health Network (RHN) Tissue Bank of the FRQS for their valuable assistance.

Y. Bossé holds a Canada Research Chair in Genomics of Heart and Lung Diseases. This study was supported by grants from the Chaire de Pneumologie de la Fondation JD Bégin de l'Université Laval, the Fondation de l'Institut Universitaire de Cardiologie et de Pneumologie de Québec, the Respiratory Health Network of the FRQS, the Cancer Research Society and Read for the Cure, and the Canadian Institutes of Health Research (MOP – 123369; to Y. Bossé).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.