The Personalized Inherited Signature Predisposing to Non-small Cell Lung Cancer in Non-smokers

1. Introduction

Lung cancer (LC) is an important public health problem and a major cause of cancer death worldwide [1]. Despite the strong relationship to tobacco smoking, less than 1 out of 5 heavy smokers develop lung cancer [2], whereas 25% of patients with lung cancer did never smoke [3], suggesting a role for genetic factors in individual susceptibility to lung cancer. So far, no specific environmental or genetic risk factors have been detected in these individuals.

Previous genome-wide association studies (GWAS) have identified up to 45 susceptibility loci for lung cancer in population-based series, mainly including smokers [4,5,6,7,8,9,10,11,12,13,14,15]. On the other hand, several studies found the rare germline T790M mutation in EGFR correlated with familial clustering of lung cancer [16,17,18,19,20,21,22,23]. The EGFR is a major cancer predisposition gene with an estimated 31% risk of lung cancer development in non-smoking carriers [19]. However, the EGFR T790M mutation is weakly oncogenic and became significantly enhanced when associated with a second activating mutation [18] more frequently with female sex and non-smoking status.

Still today, there is a need to better understanding of the genetic factors associated with lung cancer in never-smokers’ patients in non-predisposed families.

In this frame, using the whole-exome sequencing approach, in a previous study analyzing discordant sib-pairs, we identified a combination of germline disruptive mutations in known cancer susceptibility genes in never smoker with lung adenocarcinoma (ADCA) of early onset, that were not present in their healthy control sibs [24]. These results demonstrated that each affected subject had an oligogenic combination of disrupted cancer predisposing genes. This evidence gives experimental support to a model of “private genetic epidemiology” for lung cancer susceptibility that has previously only been hypothesized. The oligogenic nature of the model may therefore explain the non-heritability of the condition.

In the present study, we used a combination of genetic technical tools (whole-exome sequencing analysis and RNA sequencing) coupled with a pedigree model in discordant pairs of non-smokers with lung cancer of early onset. Healthy sibs were used as controls. We identified in affected subjects a unique combination of private “predisposing signatures”, that further confirms and exploits the oligogenic model of the disease.

2. Materials and Methods

2.1. Study design and samples

This study analyzed four non-smoker lung cancer patients with lung cancer (cases), who underwent lung lobectomy in the Thoracic Surgery Unit at the Azienda Ospedaliero-Universitaria Senese (AOUS, Siena, Italy), in comparison with their healthy sibs (controls). The genetic comparison of discordant sibs, that share 50% of the genome, facilitates the identification of variants associated with lung cancer susceptibility. For each patient, formalin-fixed paraffin-embedded (FFPE) samples of tumoral and non-tumoral lung tissues were obtained from the Pathology Unit of the recruiting hospital and analyzed through whole-exome sequencing and RNA sequencing (RNA-seq).

Each patient and sib signed the informed consent declaration for the use of their biological samples and clinical data for research purposes. The study protocol was approved by the Institutional Ethical Committees. Information on histological diagnoses (by the Pathology Unit) was retrieved from the clinical records.

2.2. DNA and RNA extraction

Genomic DNA of cases and controls was isolated from EDTA peripheral blood using QIAamp DNA Blood Kit (Qiagen®, Hilden, Germany) according to the manufacturer’s protocol. DNA was extracted from FFPE lung tumoral and non-tumoral tissue samples using MagCore^® Genomic DNA FFPE One-Step Kit for MagCore^® System (Diatech Pharmacogenetics srl, Jesi, Ancona, Italy). RNA was extracted from FFPE lung tumoral and non-tumoral tissue samples using High Pure FFPE-Tissue RNA Isolation Kit (Roche, Basel, Switzerland) following the manufacturer’s instructions. RNA samples were processed to remove ribosomal RNA (rRNA) using Ribo-Zero rRNA Removal Kit for Human samples (Illumina, Grand Island, USA) following the manufacturer’s instructions. RNA integrity was verified using the Agilent Eukaryote Total RNA Nano Kit (Agilent Technologies, Palo Alto, CA) on Agilent2100 Bioanalyzer (Agilent Technologies). Both DNA and RNA were quantified by spectrophotometry (ND-2000c; NanoDrop Products, Wilmington, DE) and Qubit® Fluorometer with Qubit® dsDNA HS Assay and Qubit® RNA HS Assay Kits (LifeTechnologies), respectively.

2.3. Whole Exome Sequencing and data analysis

Whole exome sequencing was performed using the Illumina Nextseq 550 on genomic DNA samples of cases and controls and tumor tissues as previously described [25]. After variant annotation against external datasets, including 1000 genomes (http://www.1000genomes.org/) and dbSNP, in order to identify susceptibility variants, we selected for rare variants (minor allele frequency – MAF ≤ 0.01) with a potentially damaging effect and with a functional link to cancer development. Additional filtering procedures were thus implemented for: i) retrieving exonic rare variants with a potential detrimental impact on protein function, i.e., truncating frameshift insertion and deletion and nonsense variants predicted as deleterious that were present in affected but not in unaffected sibs and vice versa; ii) identifying somatic mutations present in tumor tissues.

2.4. RNA sequencing and data analysis

RNA sequencing was performed using Illumina HiSeq2500 platform (Illumina), in a 2x100bp paired-end (PE) configuration in High Output mode (V4 chemistry), with a total of at least 250 million reads per lane. After quality check, RNA (50 ng) was used to prepare sample libraries. Sequencing library construction included these steps: library construction using Illumina TruSeq RNA Sample Pre Kit (Illumina), library purification using Beckman AMPure XP beads (Beckman Coulter s.r.l., Milan, Italy), insert fragments test using Agilent High Sensitivity DNA Kit on Agilent 2100 Bioanalyzer (Agilent Technologies), quantitative analysis of library (ABI 7500 real time PCR instrument; KAPA SYBR green fast universal 2 9 qPCR master mix. GRN), and cBOT automatic cluster (TruSeq PE Cluster Kit v3-cBotHS). Post-library quality controls were performed using the Agilent RNA 6000 Nano kit (Agilent Technologies) on Agilent2100 Bioanalyzer (Agilent Technologies) and Qubit® Fluorometer (Life Technologies). Libraries were then loaded on HiSeq2500 sequencing platform (Illumina) and sequenced using 2x100bp pair-end High Output Mode (V4 chemistry) per lane. The reads generated on the HiSeq2500 were provided under FASTQ format.

Sequence reads in FASTQ format were processing using the Fastqc software (http:// www.bioinformatics.babraham.ac.uk/projects/fastqc/) for data quality check and removing excess adaptors to get high-quality and clean reads. The high-quality reads were aligned to the GRCh38/hg38 reference human genome (ftp://jgi-psf.org/pub/compgen/phyto zome/v9.0/Ptrichocarpa/assembly/Ptrichocarpa_210.fa.gz) using the TopHat software (version 2.0.9) [26]. Transcript assembling and expression quantification were carried out using Cufflinks (version 2) [27]. Gene expression was expressed as fragments per kilo-base transcript per million mapped reads (FPKM) value [28]. This normalized value was used for visualization on a genome browser (http://genome.ucsc.edu/) [29], as well as to compare read coverage between and throughout different genes. Statistical analysis was performed to compute the mean FPKM level with the associated P-value for lung normal tissues together with the mean FPKM level with the associated P-value for lung cancer tissues. Cuffdiff tool from Cufflinks was used to identify differentially expressed genes [30]. Potential gene fusion events were detected by TopHat-fusion [31] with spanning reads ≥ 10 in cancer and normal tissue. The cancer-specific fusion genes were obtained by excluding the fusion genes that were also identified in distant normal tissue. Gene Ontology (GO) and pathways analyses were performed using the Database for Annotation, Visualization and Integrated Discovery functional annotation tool (DAVID Bioinformatics Resources 6.7, https://david.ncifcrf.gov).

3. Results

3.1. Exome Analysis of Constitutive DNA

We carried out whole-exome sequencing of DNA from blood tissue in four discordant sib pairs (Table 1). Patients had early onset lung cancer in the absence of smoking history, and we used as controls their unaffected sibling. Exome sequencing generated a mean of 28,960,442 reads for sample with a mean read length of 160bp.

Overall, we identified a total of 349,411 genetic variants in eight samples. After removing variants with low coverage and filtering for exonic mutations with MAF≤0.01 or not reported, we obtained 6,086 total variants (Figure 1). We then filtered excluding variants with clinical significance as “benign” or “likely benign” and present in an in-house database of variants identifying 3,235 variants. Among them, we selected deleterious variants applied the Combined Annotation Dependent Depletion (CADD) and the MetaSVM (Support Vector Machine) bioinformatics tools. In this way, we obtained 961 potential deleterious variants of which 370 were non-synonymous and 591 were truncating variants (insertions and deletions (indels) causing exonic frameshifts, and nonsense mutations leading to truncated proteins).

Out of these, 40 variants mapped in genes that have been found previously associated with cancer of lung or other tissues and, interestingly, all of these were present exclusively in the affected sib (Table 2). Validation of these variants was carried out using custom NGS panel for Ion PGM sequencer (Life Technologies) and led to the confirmation of 40 variants probably responsible of lung cancer susceptibility in our cases. All the 40 sequence variants identified were in heterozygous state. No common variants to all four cases were found. However, three variants were common to two out of four patients (ANGPLT4, CARS, ESRRA). The two variants in ANGPLT4 (Angiopoietin Like 4) and ESRRA (Estrogen Related Receptor Alpha) genes were common to case 1 and case 4, while the same variant in CARS (Cysteinyl-TRNA Synthetase) gene was found in case 3 and case 4. Each of the 40 variants were also present in the relative lung tumor tissue in heterozygous state.

3.2. Somatic Mutations in Tumor Lung Tissues

In the 4 tumor lung samples, we identified an average of 2.568 somatic mutations per tumor (ranging from 2.169 to 3.348). When we limited the mutations in the coding regions, the average number of mutations in each tumor was 1.510 (range 1.172-1.865), among which an average of 1092 (range 797-1542) were missense, 374 (range 85-735) were frameshift ins/del and 43 (range 32-56) were nonsense mutations (Figure 2). The number of mutations was not associated with clinical and pathological variables (stage and age at diagnosis). No recurrent mutations, such as mutations in EGFR, KRAS or AKT genes, were present in our tumor samples.

3.3. RNA-seq Analysis of Tumors and Non-Involved Lung Tissues

To detect differences in the expression level, splicing pattern and/or polyadenylation sites that could both help the understanding of the functional role of germline variants in the development of lung cancer and shade light on the pathogenic mechanisms, we performed RNA sequencing (RNA-seq) analysis of RNA extracted from both FFPE tumor and non-involved lung tissues of the four lung cancer patients. RNA-seq generated a mean of 73,865,049 reads per sample and 91.07% of bases sequenced above Q30 quality score (Supplementary Table 1).

To assess the potential functional role of the 40 germline variants identified by WES (whole-exome sequencing) and predicted as deleterious by bioinformatics tools, we combined results from WES experiments with the respective expression profile from RNA-seq. In 16 variants mapping in 16 genes, RNA sequencing data reinforced the pathogenic role of the identified variants showing three different effects (Figure 3). Firstly, 8 genes (ACACA, ANGPTL4, BUB1B, FBN2, MEN1, MMP14, TP73, WWTR1) showed a downregulation in lung cancer tissue indicating a possible “second hit” in tumor suppressor genes responsible of gene inactivation (Figure 3, panel A). Secondly, 3 genes (ACAP2, ENO3, PSCA) showed an upregulation in lung cancer tissue suggesting their role as oncogene (Figure 3, panel B). Lastly, 5 remaining genes (AMIGO3, CARS, DEPTOR, IQGAP2, RNASEL) resulted as downregulated in both normal and cancer tissue indicating a possible transcript instability (Figure 3, panel C). In addition, a nonsense variant that was predicted as deleterious in the putative tumor suppressor URI1 gene, although it was not associated with changes in mRNA levels, was also retained in the panel of candidate genes.

Data obtained from tumor tissues were then compared with those from normal tissues in order to identify differences in the expression level using well-established and accepted analysis tools such as Cufflink-Cuffdiff [30]. Cuffdiff differential expression analysis identified 315 genes significantly down- regulated (fold≤-2 change, P-value ≤0.05) and 347 genes significantly up-regulated (fold change≥2, P-value ≤0.05) in lung tumor tissues compared to normal tissues. We then interrogated these data both for the enrichment of genes involved in peculiar cell functions and for the involvement of specific pathways. GO analysis of upregulated differentially expressed genes revealed a statistically significant enrichment for genes mainly involved in cellular process (p-value=1.29x10^-4), cellular component organization (p-value=1.78x10^-4) and developmental process (p-value=0.002) (Figure 4A). GO terms of downregulated differentially expressed genes were mainly related to acute inflammatory response (p-value=0.015) and cell adhesion (p-value=0.028) (Figure 4B).

Pathways analysis individually performed on the differentially expressed genes between each normal-tumor tissue pair, showed the ECM-receptor interaction pathway as the common involved pathway in all four lung tumor tissues (Table 3).

3.4. Analysis of Tumor Loss of Heterozygosity

We performed whole-exome sequencing of DNA from the four tumor tissues in order to identify potential driver genes. Analysis of exome sequencing in DNA from blood compared to exome sequencing in DNA from tumor tissues allowed to identify the presence of 9 variants in 9 genes that were in heterozygous state in DNA from patient’s blood and in homozygous state in tumor tissue (Figure 5), thus representing a possible second hit responsible of gene inactivation in this tissue.

3.5. Protein-Protein Interaction

To explore possible pathogenic mechanisms in lung cancer, we further investigated the existence of possible interactions among the involved cancer genes that had germline mutations. We found three networks involving mutated genes belonging to different cases (Figure 5). The main network connecting EPHB6 gene (mutated in case 1) with ACACA and ENO3 genes (mutated in case 3), CARS gene (mutated in case 3 and 4) and ACAP2 gene (mutated in case 2). Overall, each patients had deleterious germline variants in genes belonging to this network.

Expanding the analysis of this network with a threshold of 20 interactors, we identified associations between the mutated genes in our patients with many genes involved in cancer development. In particular, the network clustered in two relevant groups, the former involving PRKA interactors, in which are comprised our mutated genes, and the latter involving RAD51 interactors. (Figure 6).

4. Discussion

In this study, we used an integrative approach of next-generation sequencing technologies to dissect genetic susceptibility to lung cancer in nonsmoker patients for the identification of a genetic profile that could be predictive of the individual risk for lung cancer. Whole-exome sequencing technique, coupled with a model of cases and controls deriving from the same kindred, demonstrated that each patient has a combination of an average of 10 (range 7-14) deleterious “private” germline variants in tumor predisposing genes. These mutations were absent in unaffected sibs.

In addition to perform a genome-wide analysis in the search of oligogenic signatures that could differentiate affected from non-affected siblings, we also analyzed, in the same patients RNA-seq data in tumor specimens, comparing tumoral vs non-tumoral tissue. Using this approach, we confirmed the potential functional effect of most of the identified variants with an average of 6 (range 4-8) variants per patients. In particular, we distinguished three class of alterations: 1) variants associated with downregulation of the gene in lung tumoral tissue compared to non-tumoral tissue indicating the presence of a “second hit” in a putative tumor suppressor gene; 2) variants associated with upregulation in tumoral compared with non-tumoral tissue of the gene that could represent a likely oncogene; 3) variants associated with downregulation in both tumoral and non-tumoral tissue indicating possible mechanisms for transcript instability.

Concerning oncogenic genes, the comparison of RNA-seq data in the tumoral vs non-tumoral tissue showed the involvement of genes belonging to molecular localization, cell movement and cellular component assembly pathways among the upregulated genes while, among the downregulated genes, the mainly involved pathways were cell localization, transport, negative regulation of cellular processes and response to stress. These results seem to suggest that cancer cells miss proper intracellular molecular localization and are more resistant to stress. Comparing differentially expressed genes in matched tissue pairs, we showed involvement of the ECM-receptor pathway in all of four pairs of siblings. Similar findings have already been reported in prostate cancer [32].

Network analysis concerning the identified cancer susceptibility genes showed three networks that were shared by different patients, suggesting a possible common path for private oligogenic signatures. The main network involving genes mutated in all four patients (ACACA, ACAP2, ENO3, EPHB6, CARS) connected the candidate mutated genes with the group of RAD51 (RAD51 recombinase) and PRKA interactors. RAD51 encodes for a key recombinase that seems to be essential for homologous recombination and DNA repair [33]. It interacts with a large number of proteins involved in double-strand DNA breaks and in cell cycle with important implications in tumorigenesis, such as BRCA1, and TP53 (as reviewed in 34). PRKA encodes for the catalytic subunit of AMP-activated protein kinase (AMPK) that is a cellular energy sensor that maintains energetic homeostasis [35] and has been suggested as a novel target for anticancer therapy since its activation determines a reduction of mRNA translation and protein synthesis [36]. Figure 5 and Figure 6 show the possible network (Figure 5) and “expanded” network (Figure 6) of proposed cooperation among the main mutated genes in the occurrence of “increased susceptibility” to non small cell lung cancer. We have previously discussed the putative role of each gene in the occurrence of lung cancer.

Our study has a number of strengths and limitations. One of the main strengths is that we used an integrative next-generation sequencing approach combine whole-exome sequencing technique (germline and tumor DNA) with transcriptome sequencing (tumor and matched normal tissue). In addition, the original selection strategy of discordant sibs has been used. Among the limitation, variants in noncoding regions, copy number variations, genome rearrangements and common susceptibility variants have been missed due to the study design. Secondly, the number of samples analyzed is relatively small. However, the private nature of the oligogenic susceptibility reduces the potential contribution of additional discordant sib-pairs.

Anyway, the “novelty” that present results strongly suggest is the following. Opposite to what occurs in inherited multitumoral syndromes, such as Familial Adenomatous Polyposis (FAP), where germline mutations of a single gene (i.e. APC gene) determines colon cancer in 100% of affected patients, in the present model, that is likely to occur in most frequent sporadic cancers, variable oligogenic combinations of germline mutations, that change from an individual to another, all together are responsible for the occurrence of a “susceptible” or “resistant” phenotype towards a given cancer.

Overall, our findings showed that private oligogenic signatures could be part of an individual way to cancer. We suggest that every patient has his/her peculiar signature of germline mutations, governing personal cancer susceptibility. This signature may play a major role in the development and growth of lung cancer, namely in nonsmoker patients and may therefore explain the non-heritability of the condition. In fact, lung cancer in non-smokers rarely shows familial aggregation but rather seems sporadic or occurs in phratries. These two possibilities resulted perfectly explained by our private oligogenic model of inheritance (Figure 7). The proposed model may have important implications in the evaluation of individual risk for lung cancer and may lead to a “germline genetic signature”, which, in the modern era of personalized medicine, could be of benefit for early detection of cancer.

5. Conclusions

In conclusion further study is mandatory to confirm present findings in larger studies. Anyway, our “focused analysis” in a small number of patients could contribute to a deeper insight into the complexity of the various subtypes of lung cancer and the variable interplay among gene programs involving biological processes, which could seem apparently distinct and/or distant from each other.