Evolutionary selection identifies critical immune-relevant genes in lung cancer subtypes

In an evolving population, proliferation is dependent on fitness so that a numerically dominant population typically possesses the most well adapted phenotype. In contrast, the evolutionary “losers” typically disappear from the population so that their genetic record is lost. Historically, cancer research has focused on observed genetic mutations in the dominant tumor cell populations which presumably increase fitness. Negative selection, i.e., removal of deleterious mutations from a population, is not observable but can provide critical information regarding genes involved in essential cellular processes. Similar to immunoediting, “evolutionary triage” eliminates mutations in tumor cells that increase susceptibility to the host immune response while mutations that shield them from immune attack increase proliferation and are readily observable (e.g., B2M mutations). These dynamics permit an “inverse problem” analysis linking the fitness consequences of a mutation to its prevalence in a tumor cohort. This is evident in “driver mutations” but, equally important, can identify essential genes in which mutations are seen significantly less than expected by chance. Here we utilized this new approach to investigate evolutionary triage in immune-related genes from TCGA lung adenocarcinoma cohorts. Negative selection differs between the two cohorts and is observed in endoplasmic reticulum aminopeptidase genes, ERAP1 and ERAP2 genes, and DNAM-1/TIGIT ligands. Targeting genes or molecular pathways under positive or negative evolutionary selection may permit new treatment options and increase the efficacy of current immunotherapy.

In an evolving population, proliferation is dependent on fitness so that a numerically dominant population typically possesses the most well adapted phenotype. In contrast, the evolutionary "losers" typically disappear from the population so that their genetic record is lost. Historically, cancer research has focused on observed genetic mutations in the dominant tumor cell populations which presumably increase fitness. Negative selection, i.e., removal of deleterious mutations from a population, is not observable but can provide critical information regarding genes involved in essential cellular processes. Similar to immunoediting, "evolutionary triage" eliminates mutations in tumor cells that increase susceptibility to the host immune response while mutations that shield them from immune attack increase proliferation and are readily observable (e.g., B2M mutations). These dynamics permit an "inverse problem" analysis linking the fitness consequences of a mutation to its prevalence in a tumor cohort. This is evident in "driver mutations" but, equally important, can identify essential genes in which mutations are seen significantly less than expected by chance. Here we utilized this new approach to investigate evolutionary triage in immune-related genes from TCGA lung adenocarcinoma cohorts. Negative selection differs between the two cohorts and is observed in endoplasmic reticulum aminopeptidase genes, ERAP1 and ERAP2 genes, and DNAM-1/TIGIT ligands. Targeting genes or molecular pathways under positive or negative evolutionary selection may permit new treatment options and increase the efficacy of current immunotherapy.

Introduction
Emergence of a clinical cancer population is typically a prolonged Darwinian process in which the malignant cells evolve phenotypic properties that adapt to tissue growth constraints and host immune responses. These evolutionary dynamics occur at the level of phenotypic interactions with environmental selection forces and accumulating genomic changes provide a record of this evolutionary arc. For example, during immunoediting, mutations that protect tumor cells from the immune system enhance proliferation and, therefore, have increased prevalence in the population. In contrast, mutations that increase recognition by the immune system, because they decrease fitness and proliferation, are typically lost from the population. Current molecular biology techniques prioritize the study of "driver" mutations that are observed within some large fraction of the tumor population. These provide important insights into critical pathways and potential therapeutic targets. However, prior theoretical studies have found identifying mutations subjected to negative selection can be equally informative. However, lost mutations cannot be measured directly in a tumor sample and require novel approaches.
In prior studies, the principle of "evolutionary triage (Gatenby et al., 2014)", which links the prevalence of any gene mutation in a population with its contribution to cancer cell fitness (Gatenby and Frieden, 2002;Gatenby and Frieden, 2004). A mutation that increases fitness will also increase proliferation so that it is observed more frequently either within a tumor population or in a cohort of patients with that tumor. A mutation that does not alter fitness will not change proliferation and, therefore, be observed with a frequency that reflects the underlying mutation rate (Kimura, 1968;Gatenby and Frieden, 2002;Gatenby and Frieden, 2004). These represent the well-recognized dynamics of "driver" and "passenger" mutations. Computer simulations also demonstrated some Methodology applied to the TCGA data set to identify genes that are mutated more or less frequently than expected by chance alone. (A) Evolution is a process of neutral (grey circle), positive (green circle), and negative selection (red x-mark) events resulting in the accumulation of mutations that increase cell fitness and removal of mutations that decrease cell fitness. (B). A robust regression model was used to establish the mutation rate based on the fraction of patients with a protein altering mutation in each gene against the size of the gene. Distance from the line for each gene was used to quantify conserved and hypermutated genes. Genes on the regression line have background levels of mutations. (C) Distance from the regression model for each gene. Genes with an SR ≥ 1.0 are hypermutated and genes with an SR ≤ -1.0 are conserved genes (under mutated). All others have an expected number of mutations (Background mutations).
genes must continue to function normally for optimal cancer cell fitness. In these genes, non-synonymous mutations cannot improve function and, so, will frequently decrease fitness and disappear from the population. These "essential" genes will usually have an observed prevalence of mutations that is less than expected by chance alone (i.e., less frequently than a passenger gene of comparable size) ( Figure 1A).
The dynamics of negative selection are well recognized, but prior investigations across multiple cancer types found few genes are broadly conserved (Martincorena et al., 2018;Zapata et al., 2018). However, computer simulations demonstrated evolutionary conservation in cancer populations varies depending on local tissue environment selection forces and prior molecular changes in their evolutionary arc (Gatenby et al., 2014;Cunningham et al., 2015). Thus, genes conserved among cancers emerging from different tissue types and/or through different initiating driver genes are predicted to be uncommon. Recent studies of TCGA data in EGFR-mut, KRAS-mut, and wild-type (no known driver mutations) (WT) lung adenocarcinomas, confirmed this prediction as few genes were found to be conserved even among subtypes that shared a common tissue of origin (Freischel et al., 2021).
Here we investigate evolutionary triage in genes associated with tumor-immune interactions. Epithelial cells are active participants in the host immune response (Schleimer et al., 2007). Antigen presentation, immune-modulating protein expression (e.g., checkpoint ligands), as well as cytokine and chemokine secretion by epithelial cells can promote or reduce host immune functions (Schleimer et al., 2007;Larsen et al., 2020). Additionally, epithelial cells express surface and cytoplasmic pattern recognition receptors (PRR), which recognize molecular patterns derived from foreign pathogens and markers of stressed cells triggering intrinsic immune responses within epithelial cells, including release of inflammatory mediators and initiation of pyroptosis (Amarante-Mendes et al., 2018).
We hypothesized that, within these diverse and complex interactions, we could determine the most important cellular functions and molecular pathways by identifying genes under strong positive and negative evolutionary selection. Furthermore, by investigating cancers arising from the same organ but with different initiating oncogenic mutations, we could investigate divergence in immune-evasion strategies arising from different initiating genetic events.
Here we explore this hypothesis by investigating the mutational frequency of genes with known or expected functions in host-immune interactions in two cohorts of lung adenocarcinoma patients in the TCGA database, mutant KRAS (n = 163) and no known driver (WT) (n = 313). Evolutionary selection of multiple genes with well-documented roles in LUAD supports our hypothesis and encourages further investigation of other evolutionarily identified genes and molecular pathways that have not been extensively investigated.
We identified one sample each from patients belonging to the lung adenocarcinoma cohort, and classified patients based on known driver or recurrent mutations in KRAS (G12, G13, Q61, A146), BRAF (V600, N581, G464, G466, G469, G596, D594), and EGFR (L858, S768, L861, G719, T790, indels in exons 18-21). Samples were excluded if they matched criteria for more than one of these genes. Samples that did not meet any of the mutation criteria were classified as WT (no driver mutations in EGFR, KRAS, or BRAF.) Cohorts with EGFR-mutant LUAD (n = 58) and BRAF-mutant LUAD (n = 31) had fewer samples, so we elected to focus primarily on the mKRAS (n = 163) and WT (n = 313) cohorts. Tumor and normal sequence alignment files (BAM) were downloaded from the Genome Data Commons, and gene-level depth of coverage was calculated by calculating bases covered by sequencing from the above files across each of the RefSeq coding genes (with 25 base pair flanking regions). A base was considered sufficiently covered if the depth of coverage was≥14 in tumor sample and≥8 in normal samples (as has been previously described: https://www. synapse.org/#!Synapse:syn1695394). The fraction of each gene's protein coding bases (using the longest RefSeq transcript) covered by sufficient sequence data was calculated for each sample using the Negative Storage Model [PMC7157186].

Identification of over and under-mutated genes
To identify over-and under-mutated genes, we calculated the fraction of patients in each cohort with any protein-altering (nonsynonymous, truncating, canonical splice-site) mutations in each gene. To control for undercounting of mutations due to poor coverage, we divided the fraction patients mutated by the median fraction of coding bases with sufficient coverage. Genes with median coverage <50%, ANNOVAR annotation errors, or gene expression <2.0 log2 counts were excluded, as these genes either have artificially low mutation rates or are not likely to be essential due to low expression. Gene size was controlled by calculating a linear model between corrected mutation fraction and gene size (protein coding bases from largest RefSeq transcript.) ( Figure 1B). Corrected mutation rate and gene size were both square-root transformed. The standardized residual of each gene to the regression line was determined: the most positive genes are the classical over-mutated or driver genes, and in each cohort the most highly mutated gene had Frontiers in Genetics frontiersin.org 03 the highest standardized residual (e.g., KRAS was the most overmutated gene in the mKRAS cohort.) Genes with an SR ≥ 1.0 are overmutated and genes with and SR ≤ -1.0 are conserved genes ( Figure 1C). To minimize false positives, our analysis primarily focuses on genes with standardized residual (SR) values for conserved genes SR ≤ -2.0 to achieve p = 0.05 and SR ≤ -1.65 to achieve p = 0.1.

Gene conservation differs between the wild-type and mutant KRAS cohorts
More genes were highly conserved (SR ≤ -2.0) in mKRAS than WT cancers (160 and 248 respectively) (Supplementary Table S2). Among the immune genes examined (n = 576) 19 were conserved (SR ≤ -1.65) in WT cohort and 11 in mKRAS (Table 1). Many of the conserved genes were unique to each cohort with only one overlapping pair of genes. Complement genes C4A (WT SR = -2.863, mKRAS SR = -1.711) and C4B (WT SR = -2.864, mKRAS SR = -1.711), which encode the acidic and basic forms of complement factor 4, an activator of the classical complement pathway, were conserved in both cohorts. No genes from the list of chemokines and chemokine receptors were highly conserved in either cohort (Supplementary Figures S1A,B).

Differential conservation in antigen processing genes
Endogenous proteins are degraded into antigen precursors in the cytoplasm by the (immuno)proteasome. Transporter associated with antigen processing (TAP) transports the protein fragments into the endoplasmic reticulum (ER) for further processing by aminopeptidases and loading onto MHC-I proteins. Dysregulation of this process alters the peptide pool presented and loss of integral component results in diminished MHC-I expression. Transporter 1, ATP Binding Cassette Subfamily B Member, TAP1 was conserved in the WT cohort (SR = -2.566) but not in the mKRAS group (SR = -0.561) ( Figures 3A,B). TAP1 co-localizes with TAP2 and mediates the flow of peptides from the cytosol into the endoplasmic reticulum prior to cleavage by ERAP proteins. TAP2 gene expression and conservation score could not be determined from our dataset.
The complexity of evolutionary selection is notable in the endoplasmic reticulum aminopeptidase genes, ERAP1 and ERAP2, which shape the peptidome by cleaving peptide precursors after transport into the ER ( Figure 3C). ERAP1 is among the most conserved genes in WT cancers (SR = -2.797) and ERAP2 among the most conserved in mKRAS (SR = -2.085) ( Figures 3A,B). Interestingly, no ERAP1 mutations were found in any mEGFR tumors or comparable TCGA data from melanoma patients (data not shown), suggesting its function is essential in multiple cancer types.

Polio virus receptor gene conservation in wildtype tumors
Tumor cells express immune modulating proteins that activate or inactivate (checkpoint proteins) immune cells. None of the immune modulating proteins were significantly conserved in the mKRAS group. Polio Virus Like Receptor 2 (PVRL2, Nectin2, CD112) and PVR (CD155), conserved in WT (SR = -2.054 and -1.782, respectively), were also conserved in mKRAS cohort although did not reach significance (SR = -1.458 and -1.234, respectively) ( Figures 4A,B). The conserved tumoral PVR proteins represent known tumor evasions strategies and can activate or inhibit (Wu et al., 2021) NK cells and T cells ( Figure 4C).
MHC class I polypeptide-related sequence A (MICA) a stress-induced activating ligand for Natural Killer Cell Lectin Like Receptor, KLRK1 (NKG2D) was conserved in the wild-type cohort (SR = -1.699). Ligation of NKG2D activates NK cells and gamma-delta-T cells and co-stimulates CD8 T cells (Holdenrieder et al., 2006). However, MICA can be cleaved from the cell surface and inactivate distant NK cells, NKT cells, gamma-delta T cells, and CD8 T cells (Uhlenbrock et al., 2014). Numerous proteases, including MMPs, ADAM10, and ADAM17, are involved in MICA cleavage. Additional NKG2D ligands, RAET1E (ULBP4) and ULBP2, were also conserved in the wild-type cohort but did not reach the SR cutoff of -1.65 (WT SR = -1.37 and -1.32, respectively).

Conservation of cytokine signaling genes
Malignant epithelial cells secrete and detect cytokines in the tumor microenvironment that regulate the immune response and support pro-growth signals. The potential for false positives is increased in genes with a conservation score above SR = -1.65. However, we note that both subunits of the type I IFN receptor, interferon-alpha receptor 1 (IFNAR1) and IFNAR2, are conserved in mKRAS samples (SR = -1.490 and SR = -1.418, respectively) ( Figures 5D,E). Stimulation of IFNAR1/2 activates Janus kinase 1 (JAK1), (neutral in mKRAS) and tyrosine kinase 2 (TYK2) (mKRAS SR = -1.070). Type II interferon-gamma receptors IFNGR1 and INFGR2 were also conserved in mKRAS but did not meet the -1.65 threshold (SR = -1.371 and SR = -1.069, respectively) ( Figures 5D,E). IFN-γ signaling in epithelial cells induces the expression of multiple genes containing γ-interferon activated sequences. MHC-I, checkpoint inhibitor PD-L1, and immunosuppressive enzyme, IDO are upregulated on cancer cells by IFN-γ treatment (21). All interferon receptor genes were neutral in the wild-type cohort.
Two of the nine interferon regulatory factor (IRF) genes were highly conserved and expressed in the wild-type cohort. IRF5, conserved in WT (SR = -2.002) and neutral in mKRAS, regulates type-I interferon genes and IFN stimulated genes (ISG) ( Figures  5C,D). Similarly, activated IRF-3, conserved in WT (SR = -1.864) translocates to the nucleus following phosphorylation by TKB1 (WT SR = -1.010) and IKBKE (WT neutral) kinases and induces expression of type I interferons and ISGs. IRF3 is activated by a variety of pattern recognition receptors which are neutral in both cohorts ( Figures 5A,B,F) (Zhu et al., 2010). Interleukin 1 receptor-1 (IL1R1) is the most conserved among the cytokine signalling genes in WT tumors (Supplementary Figure S1C). Following ligation by either IL1A or IL1B, IL1R1 associates with interleukin-1 rector accessory protein (IL1RAP) (WT SR = -1.116) and activates many downstream proinflammatory pathways including NF-κB and MAPK (44).

Over-mutated genes are shared between WT and mKRAS cohorts
Mutations can inhibit or enhance protein function. Positive selection occurs when the mutation increases fitness. WT and Frontiers in Genetics frontiersin.org 06 Frontiers in Genetics frontiersin.org 07 mKRAS tumors averaged 348 and 300 total mutations per tumor compared mEGFR LUADs with 86 (data not shown). Frequently mutated genes, defined as genes mutated in greater than 10% of patients in each cohort, numbered 208 and 157 in the WT and mKRAS cohorts, respectively. Genes with the highest mutation scores in our cohorts reflect known drivers and frequently mutated genes including KRAS, TP53, KEAP1, and STK11 (Supplementary Table S3) (Gong et al., 2020), demonstrating

FIGURE 3
Gene conservation in antigen presentation genes. Gene conservation score, measured as the standard residual from the regression line (SR) graphed against genes expression levels (Log2) for wild-type (no identifiable driver mutations) (WT) (A) and mutant KRAS (mKRAS) (B). Conserved genes (SR ≤ -1.0) and over-mutated genes (SR ≥ 1.0) are colored and labelled with text. ERAP1 and ERAP2 sculpt the peptide pool prior to presentation by MHC-I proteins. Dysfunctional ERAP alters the peptides presented (C). Standardized residual (SR). Created with Biorender.com.
Frontiers in Genetics frontiersin.org the validity of this approach. Of the immune related genes we examined (n = 576), 29 were over-mutated in WT (SR ≥ 1.65) and 26 in mKRAS with 16 genes common between the subtypes (Table 2).

FIGURE 5
Gene conservation in pro-inflammatory genes. Gene conservation score, measured as the standard residual from the regression line (SR) graphed against genes expression levels (Log2). Innate immune signaling genes in WT (A) and mKRAS (B). Interferon signaling genes in WT (C) and mKRAS (D). Conserved genes (SR ≤ -1.0) and over-mutated genes (SR ≥ 1.0) are colored and labelled with text. Type-I and Type-II interferon receptors are conserved in mKRAS (E). Interferon regulatory factor 3 (IRF3) is a central component in multiple pattern recognition receptor pathways (F). Standardized residual (SR). Created with Biorender.com.
Frontiers in Genetics frontiersin.org

Discussion
The evolutionary arc of each tumor is likely unique with variations resulting from molecular properties of the initiating cells, heterogeneity in host responses, and stochastic variations in accumulated genetic events (Toki et al., 2018;Wang et al., 2021). However, prior theoretical studies have suggested commonalities should be found among tumors originating from the same tissue and/or initiation driver genes. Here, we demonstrate mKRAS and WT lung cancers generally exhibit convergent evolutionary selection among broad mechanisms of immune interactions and some gene families, but the specific genes selected, particularly those conserved, are frequently subtype specific. Such convergence of broad (phenotypic) properties through divergent molecular pathways is observed in nature. For example, cavefish develop a characteristic phenotype (no eyes or pigment) but through diverse genetic pathways depending on the genomic state of the founder population and stochastic perturbations (Gatenby et al., 2011).
Magalhães (Magalhães, 2021), in a PubMed analysis, found virtually every human gene has been associated with cancer. Within this vast data set, how can clinically important genes be selected? Here we propose a "Darwinian road map" can identify critical genes under evolutionary selection in clinical cancer cohorts. We examine gene conservation in immune related genes. Several of the genes identified in our studies have been noted previously. Preclinical and clinical work targeting many of these pathways are ongoing. The divergence of evolutionary selection in mKRAS and WT suggest pre-clinical experiments and clinical investigations with targeted drugs should be carefully designed with an understanding of the Darwinian importance of specific genes in specific tumor subtypes. For example, among immune modulators, common targets for immunotherapy (e.g., CD274, CTLA4) show neutral patterns of evolutionary selection. PVR, PVRL2, and MICA, which likely regulate NK cell function, are conserved in wild-type cohort. Nine anti-TIGIT (PVR and PVRL2 receptor) antibody therapies are currently being evaluated in 43 trials of advanced solid tumors, including NSCLC (Ge et al., 2021).
Each cohort conserves only one of two ERAP genes, which modify peptides for presentation by MHC class I. Clearly, the ERAP function is critical for optimal fitness in each cohort, and we note ERAP1/2 expression was detected in all 9,125 tumor specimens in The Cancer Genome Atlas database (TCGA), and deep deletions are rare (0.6-0.8%) (Compagnone et al., 2019. Targeting ERAP has been suggested by others and ERAP modulators are under development to modulate the antitumor immune response (Stratikos, 2014). However, our data raise an intriguing question: why do WT tumors highly conserve ERAP1 while mKRAS tumors equally highly conserve ERAP2? IRF3, conserved in wild-type cohorts, regulates the expression of type I IFN genes and IFN stimulated genes (ISG) through binding interferon-stimulated response elements (ISRE) in their promoters. IRF3 is activated by several innate immune receptors ( Figure 5F), including the cGAS/STING pathway, a novel immunotherapy target with ongoing clinical trials in cancer (Jiang et al., 2020).
Finally, in prior theoretical studies, computer simulations predicted therapies that disrupt conserved genes can have similar therapeutic efficacy to targeting driver genes (Gatenby et al., 2014). Thus, conserved genes in this study may be valuable clinical targets. Furthermore, simulations predicted combination therapy targeting a driver gene and a driverspecific conserved gene was highly lethal and often produced an evolutionary state in which no adaptive strategy was available, resulting in population extinction (Gatenby et al., 2014).

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: mc3.v0.2.8.PUBLIC.maf.gz: https:// gdc.cancer.gov/about-data/publications/mc3-2017. Genome Data Commons: https://portal.gdc.cancer.gov/repository. affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.