Expression of collagen type 1 alpha 1 indicates lymph node metastasis and poor outcomes in squamous cell carcinomas of the lung

Microarray data

GSE40588 was accessed and downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo). The dataset contains 60 tumor-adjacent (histologically normal) noncancerous biopsy samples. 34 (56.7%) were from patients diagnosed with regional lymph node metastasis. 57 of the cases were male vs 3 female with an average age of 58.9 years. 16 (27%) of the cases were designated N1, 18 (30%) were designated N2, and 26 (43%) were designated N0. All samples were analyzed using an Agilent 4x44K gene expression microarray. Primary data are available in the supporting information.

Identification of DEG and HUB genes

Differentially expressed genes (DEGs) were identified using GEO2R (Barrett et al., 2013). Cutoff criteria was set to a P-value <0.05 and logFC (fold change) >0.5. Probe sets which were misidentified were excluded (i.e., gene symbols not in the Genome Database) and genes with two or more probe sets were averaged. The Search Tool for the Retrieval of Interacting Genes was used to construct an interaction network from the DEG’s (Szklarczyk et al., 2015; Doncheva et al., 2019). Cytoscape version 3.4.0 software was used to visualize the output sequence. Molecular Complex Detection (MCODE Bader & Hogue, 2003) was used to arrange the topology, cluster the connected genes, and search through the resulting network. Proteins were included if they had an MCODE score >5, degree cutoff of 2, node score cutoff of 0.2, node density cutoff of 0.1, k-score of 2 and max depth of 100. HUB genes were initially selected if their degree was >10. Functional analysis was performed on both the DEG and HUB genes using the Database for Annotation, Integrated Discovery, and Visualization (DAVID; http://david.ncifcrf.gov) and validated using Cytoscape’s Biological Networks Gene Oncology (BiNGO) tool (Maere, Heymans & Kuiper, 2005). HUB gene and Cytoscape data are available in the supporting information.

Identification of COL1A1 in LSCC

HUB genes were selected which have a degree of connectivity of ≥10. cBioPortal was used to screen the cancer genome atlas pan cancer atlas (Rhodes et al., 2004). This contains 487 unique lung squamous cell carcinoma samples. Following identification of COL1A1 in LSCC, OncoMine was used to evaluate overall expression levels in the sample set. Additionally, UCSC Xena (Goldman et al., 2020) was used to evaluate isoform expression levels. Furthermore, a receiver operating characteristic (ROC) curve was generated by ploting the sensitivity vs specificity of COL1A1 using The Cancer Genome Atlas (TCGA) database.

Immunostaining of clinical samples

Immunostaining of clinical biopsy samples was done using general protocols. Briefly, tissue samples were fixed using a 4% paraformaldehyde solution. Tissues were frozen and 5 µm sections were cut using a cryostat. Sections were placed on a histological slide and dried for 30 min on a slide warmer. Slides were then immediately stained to reduce chances of oxidation and epitope loss. Sections were blocked using a blocking buffer (Fisher Scientific) and then stained with a mouse anti-COL1A1 (Proteintech, China) followed by an anti-mouse secondary antibody (IHC Kit, Maixin Company, China). UltraSensitive SP was used as a detection agent following the manufacturer’s recommended protocols.

Validation in clinical samples

33 unique LSCC tissue biopsies were collected from patients with a median age of 59.09 (±7.30) undergoing lobectomy at the First Hospital of China Medical University. Collections took place between January 2012 and December 2012. The general information of clinicopathological of the patient are presented in Table 1. COL1A1 expression was detected using the aforementioned immunostaining protocol, and levels were evaluated using a general H-score protocol. Briefly, 5 regions of each sample were randomly chosen and visualized at 400x magnification. Intensities were based on a range of 0–3 where 0 is no staining, 1 is weak staining, 2 is moderate staining, and 3 is strong staining. Staining ratios were based on ranges from 0–100. H-score was calculated as the product of intensity and staining ratio (i.e., % of cells stained positive in a given region of interest). The research was approved by the Ethics Committee of the First Hospital of the China Medical University (AF-SOP-07-1. 1-01). Written informed consent was obtained from all of the patients. Raw data for human samples are available as a spreadsheet in the Supplemental Information.

Statistical analysis

Statistical analysis of in vitro work was performed using SPSS (Version 20). Unless otherwise noted, all experiments were independently replicated at least 5 times (n = 5). Data is plotted as the mean and SEM. Results were deemed statistically significant if p < 0.05 using either a Mann–Whitney test. Statistical analysis of bioinformatic work was performed and reported by the integrated bioinformatic tool. All statistical measures are noted either within the figures or within the accompanying captions.

Differential expression of collagen genes is associated with SCC

Accession of genome wide expression data is a powerful tool for evaluating biomarkers in a variety of different diseases. We analyzed GSE40588 (Cao et al., 2017), which profiles tumor-adjacent normal lung tissue from patients with squamous cell lung cancers, using GEO2R. 409 of the genes studied were differentially expressed in late stage LSCC, and 13 of these genes were highly connected HUB genes with >10 connections (Fig. 1B). Out of the identified HUB genes, COL1A1 was the most highly connected, with 31 degrees of unique connections.

We evaluated the HUB genes using Cytoscape’s Biological Networks Gene Ontology (BiNGO). BiNGO assess over-represented or under-represented gene ontology categories and creates an interactive visual map for further evaluation. Over-represented ontologies are presented as darker color labels (i.e., orange) and less represented ontologies are lighter (i.e., yellow). While several of the gene functions were related, the most overrepresented were in extracellular matrix, adhesion, and structural activities (Fig. 2A). These results were expected considering the biological function of the collagen markers analyzed. These results were corroborated using DAVID, which showed a high number of carbohydrate binding and structural molecule activity (Fig. 2B). Further evaluation using cBioPortal showed that each HUB gene studied had a high frequency of alteration (Fig. 2C). While these results suggest that mutations in these genes, particularly overexpression, could be a prognostic marker, analysis of the hub genes on UCSC’s Xena Browser shows many of the genes are also highly expressed in solid normal tissue (Fig. 2D). Several identified genes (e.g., COL4A1, COL18A1, and COL1A2 to name a few) are only marginally overexpressed in primary tumors compared to normal solid tissue. COL1A1 and COL3A1 were the most over-expressed. Based on this data, as well as the high degree of connectivity (Fig. 1B), we chose COL1A1 as the most pertinent biomarker in the LSCC microenvironment.

Overexpression of COL1A1 in the tumor microenvironment is a diagnostic marker for LSCC

We chose to further study COL1A1 as a prognostic and diagnostic marker for SCC. An examination of this gene using the UCSC Xena Browser shows that there are 20 unique COL1A1 mutations in the TCGA, with a missense mutation being the most common (Fig. 3A). In addition, we found that a majority of the genetic alterations correspond with a diploid presentation (Fig. 3B). While mutations in the COL1A1 gene may lead to other non-small cell lung cancers (e.g., large cell carcinomas and adenocarcinomas), it was nearly 2x as likely to mutated in squamous cell carcinomas (Fig. 3C). Furthermore, based on the Oncomine analysis, mean-expression levels are nearly 3.5-fold more in LSCC samples than in normal lung tissue (Fig. 3D, Corresponding p-values are provided in the figures). Comparatively, COL1A1 expression levels were higher in LSCC samples than in normal tissue samples, with a p value of 0.05. This conclusion was corroborated with GEPIA(Tang et al., 2017), which analyzes the TCGA gene expression levels. (Figs. 3E and 3F). Finally, we used computational tools to evaluate the prognostic value of COL1A1. Based on data from UCSC, we found a statistically significant difference (p < 0.05) between cancer and normal tissue. However, there was not enough difference to be able to discriminate between early and late stages of the disease. In addition, high expression of the COL1A1 gene only resulted in a marginal decline in mortality (HR = 1.33). These results suggest that COL1A1 expression may be a significant marker for identifying and validating NSCLC subtypes.

COL1A1 upregulation signals metastasis to lymph nodes

Thus far, we have determined that COL1A1 is a potential diagnostic marker for NSCLC. This was determined using a number of applied bioinformatics strategies. There are several examples of similar workflows being done throughout literature (Cui et al., 2018; Gao & Wang, 2018). However, we chose to further validate these results using clinical biopsy specimens taken from patients with LSCC as well as matching adjacent normal tissue. Results from COL1A1 immunostaining (Fig. 4A) align with our bioinformatic results, suggesting LSCC samples have a much higher level of COL1A1 than did normal tissue (Fig. 4B). While there was no significant difference due to age, treatment, or tumor size (Figs. 4C–4E) there were significant differences in lymphatic metastasis and staging. Tumors with 0 lymphatic metastasis (N0) had a median H-Score of 100, whereas tumors which had metastasized to lymph nodes had a median score of 130. Early tumors (Stage I) also had a significantly lower median H-score, 80, than later stage tumors (Stage II–IV), 140. Compared to the earlier examination of genome sets (Fig. 3G), our data suggests stages II–IV have remarkably much higher COL1A1 expression levels. To corroborate these results, we evaluated the Oncomine database which contained 156 LSCC samples (Fig. 4H). Similar to the data presented in Fig. 4G, the expression levels of COL1A1 increase as the disease progresses. To further verify these results, we examined data in the TCGA database. A ROC curve was produced (Fig. 4I), with an area under the curve (AUC) values of 0.796, indicating high classification performance. These results suggest COL1A1 is a valuable marker which can be used to assist in diagnosing metastasis of LSCC tumors to the lymph nodes.