Immunological surveillance against cancer across mammals

Domain: earth_science
Published: 2025-11-24
Editorial outcome: major_revision

Abstract
Contrary to expectations based on their higher cell numbers, larger and longer-lived species do not face dramatically increased risk of cancer. This strongly suggests that evolution has fashioned natural cancer resistance mechanisms, yet our knowledge remains limited on what these mechanisms might be. The cancer immunological surveillance hypothesis, proposed by Burnet and Thomas in the 1950s, highlights immunity as a key factor determining species-specific cancer resistance. Here we address the original, evolutionary interpretation of this hypothesis by investigating the relationship between cancer mortality risk and markers of efficient antigen presentation. Our results show that the expansion of the MHC class I gene complex, as well as increased selection for diversity at these genes is associated with sharply decreasing cancer mortality risk across mammals. This suggests that the efficient presentation of diverse peptides in somatic cells is important for cancer suppression across mammals, providing pioneering evidence that supports the cancer immunosurveillance hypothesis across species.

Introduction
Cancer is associated with mutations in genes that control the growth, division, or functioning of cells. These alterations are generally reflected by the phenotype of cells, including a set of tumor-specific or tumor-associated antigens exhibited at the cell surface that differ (in structure, distribution, or density) from those of healthy cells, making them immunogenic1. Such immunogenicity implies the recognition of tumor antigens by the adaptive immune system, potentially triggering the elimination of cells carrying mutations, often before tumor formation could take place. Critical to the functioning of this process is the major histocompatibility complex (MHC, in humans named human leukocyte antigen complex, or HLA)2. The MHC gene complex is present in all amniotes and encompasses some of the most variable functional genes described across these lineages1,3,4, with the number of duplicated genes being unusually high and their alleles exhibiting hyper-polymorphism5. Genes of the MHC code for membrane proteins whose primary role is antigen presentation in the immune response2.

MHC class I (MHCI) molecules have a particularly important role in cancer defense6. These proteins are present on the surface of all nucleated cells in vertebrates4, and serve the function of presenting intracellular protein fragments to the immune system. In healthy cells, this implies the presentation of fragments of proteins normally present in cells, eliciting no eliminative response from the immune system, through mechanisms of self-tolerance7. In case of an intracellular infection (e.g., viruses, intracellular bacteria) or genetic alterations (e.g., cancer driver mutations), peptides unusual to healthy cells are transported and exhibited on the surface of cells by MHCI molecules. Abnormal antigens are then recognized by CD8 T lymphocytes and (unless T cell responses are restrained by immune checkpoints8) cells carrying these are destroyed by CD8 effector T cells1. In case of cancer, this process is called T-cell mediated cancer immunosurveillance, and it is a mechanism widely recognized as an important pathway that nips many cancers in the bud6. The cancer immunosurveillance hypothesis is not new, having its roots in the early 20th century work of Paul Ehrlich9, later being suggested by Lewis Thomas10, and further developed by Sir MacFarlane Burnet11. The two of them explicitly framed immunological surveillance as an evolutionary hypothesis, suggesting that the need to recognize and eliminate tumor cells was the primary selective force leading to the emergence of adaptive immunity, a major evolutionary event that happened about 500 million years ago12. The role of adaptive immunity in cancer defense is now well established, with extensive experimental work demonstrating much higher rates of carcinogen-induced, as well as spontaneous tumors in mice lacking adaptive immunity as compared to wild-type mice13.

The efficiency of T cell-mediated cancer immunosurveillance is contingent on multiple factors, and the MHCI genotype could be one of these14. The generally observed high variability in MHC alleles within diverse animal taxa has been suggested to indicate the evolutionary advantage of MHC heterozygotes over homozygotes4,15,16. This is also supported by observations that codons maintaining the structural integrity of the MHCI molecules are highly conserved, while codons of the peptide-binding regions are extremely variable15,17. Heterozygote advantage might emerge due to superior protection against multiple infectious agents or mutations4,18. For instance, heterozygosity in MHC I genes at the level of individuals is associated with a larger repertoire of protein fragments that are able to bind to the peptide-binding groove of MHCI molecules3,18,19, potentially resulting in more efficient elimination of cells with abnormal intracellular proteins, such as those emerging due to infections or genetic damage16,19,20. In the case of cancer defenses, the importance of MHC class-I proteins is striking6,7. For instance, in human cancer patients, it was recently demonstrated that heterozygosity at HLA-I loci (and even nucleotide diversity of HLA-I alleles across fully heterozygous individuals) is a good predictor of the efficacy of checkpoint blockade immunotherapy19,20. These studies demonstrate that HLA-I genetic diversity plays a central role in efficient immunological surveillance of abnormal cells by increasing the repertoire of peptides presented on MHCI (the immunopeptidome). Moreover, abundant evidence suggests that suppression of MHC I expression and the resulting impairment in antigen presentation is a mechanism widely used by tumor cells to evade T cell-mediated cancer immunosurveillance7,21. The significance of MHCI in cancer prevention is evident not only at an individual level but also within populations of certain animals. For example, the suppression of classical MHC I molecules on the surfaces of cancer cells, and the presence of non-classical MHC molecules (which inherently lack variations between individuals) are believed to collectively contribute to the heightened vulnerability of Tasmanian devils to clonally transmissible cancer known as devil facial tumor disease22–24, initiated by cancer cells spreading as allografts among individuals. Due to its prominent role in cancer defense, MHC-I-mediated immune surveillance is currently also the primary target mechanism for cancer vaccines25,26. Moreover, recent evidence in humans has suggested that MHCI genotype predicts the mutational landscape of tumors in human cancer patients14. When looked across species, clear differences in MHCI genotype are evident regarding both loci number and allele diversity27,28, mostly believed to be shaped by pathogen-mediated balancing selection. Nevertheless, we have no understanding on how cross-species variability in MHC genetics is linked with cancer risk.

Cancer susceptibility varies across animals, with some species highlighted for their exceptionally low risk of cancer (e.g., naked mole rats29). Moreover, long-lived, large bodied species are often highlighted for their surprisingly low risk of cancer compared to theoretical expectations based on their large cell numbers (due to their large size and continuous cell turnover throughout their long lives)30–32. This observation, known as “Peto’s paradox”, brings into evidence that evolution solved the problem of carcinogenesis in large bodied and/or long-lived species33. A recent study has highlighted that a key mechanism through which this might happen is a slower rate at which the mutational burden of somatic tissues augments, at least in long-lived taxa, pushing cancer to mostly emerge in late age-classes34. Such deceleration in mutational accumulation during the evolution of long-lived species could occur through (1) a lowered rate at which mutations emerge (e.g., better protection of DNA), (2) an improved DNA repair and/or neutralization of damaged cells (e.g., by apoptosis/cellular senescence)35, (3) more efficient immunological surveillance and elimination of cells with non-silent mutations, or (4) a combination of these factors. In this study, we aimed to explore the third component, by testing whether the efficiency of cancer immunosurveillance (reflected by the complexity, defined here by higher gene copy number or allelic diversity of the MHC I gene complex) shapes cross-species variation in cancer mortality risk (CMR) or inherent cancer predisposition (reflected by body size or longevity) across mammals. We predict (1) decreasing cancer prevalence with increasing complexity or diversity of the MHCI gene complex across species; and (2) increasing complexity or diversity of the MHCI gene complex with increasing longevity and/or body mass in order to compensate for their increased inherent risk of cancer.

Here we built on our previous work30 in which we established a methodology to estimate CMR of hundreds of zoo-held mammalian species, using two metrics i.e., CMR and incidence of cancer mortality (ICM). CMR represents a simple measure of the prevalence of cancer mortality across deceased animals, while ICM is the incidence of cancer mortality at 90% adult longevity, thereby being controlled for age and censoring (see ref. 30 for details). In order to establish the genetic diversity at MHCI loci, we obtained sequences of classical MHCI genes from GenBank, particularly focusing on exons 2 and 3 (i.e., regions coding the extracellular α1 and α2 domains, including the peptide-binding region of the MHCI molecules). Given the heterogeneity in sampling intensity across species, traditional measures of MHC I diversity (e.g., allele diversity) could not be used here. Therefore, we estimated selection for MHC diversity as the ratio of non-synonymous (dN) to synonymous (dS) nucleotide substitution rates, a widely used measure of evolutionary pressure on the coded proteins36. Importantly, only about 25% of codons from the considered exons 2 and 3 code for the peptide-binding groove of the MHC I molecules37,38, responsible for direct interactions with antigens. These functional peptide-binding MHC residues are the most likely targets of positive diversifying selection4. Therefore, to focus on sites most likely involved in antigen recognition, we further estimated selection for MHCI diversity at a subset of codons showing the strongest evidence of positive selection. Specifically, we identified the 20 most positively selected sites (PSS) across the MHCI sequences from all considered species and used the estimates of the dN/dS ratio across this set of sites as a proxy of selection for functionally relevant (in terms of antigen binding properties) variation of MHCI molecules (hereafter referred to as selection for MHCI diversity at PSS). Moreover, in order to measure the complexity of the adaptive immune system we estimated the expansion of the MHC class I gene complex over the evolutionary timeline of mammals. We focused on a single gene family (i.e., groups of genes that share a common ancestor) that includes key genes of the MHC-I antigen-presenting system (AZGP1, MHC-A, MHC-B, MHC-C, MHC-E, and MHC-G) and we estimated the number of protein-coding genes belonging to this gene family in each mammalian species using Orthofinder39.

Results
Selection for MHC I diversity and cancer mortality risk
We obtained MHC I sequence data to estimate selection for MHC I diversity in taxa where cancer mortality could also be assessed. Data were available for 28 species for cancer mortality rate (CMR) and 27 species for intrinsic cancer mortality (ICM) (Fig. S1). Using this information, we then explored the link between selection or MHCI diversity and measures of cancer mortality. Species with zero CMR/ICM estimates were excluded from multivariate models (resulting in final sample sizes of 25 and 24, respectively), because zero cancer mortality estimates often reflect poor sampling (see ref. 30), while the excess of zeros skews the distribution of disease risk in a way that violates model assumptions. Moreover, the low number of species with no detected cancer cases (n = 3) did not allow us to run binomial models on cancer incidence. We thus analyzed non-zero CMRs as a function of selection for MHC I diversity in phylogenetic logistic regressions. We detected negative correlations between selection for MHC I diversity and both CMR (β (SE) = −3.62 (2.13), p = 0.1029), as well as ICM (β (SE) = −5.89 (2.41), p = 0.0232, Table S1 and Fig. 1a, b). The coefficient of determination indicated that selection for MHC I diversity alone explains 24% and 35% of variation in CMR and ICM, respectively. Selection for MHC I diversity at PSS showed considerably stronger associations with CMR (β (SE) = −0.73 (0.21), p = 0.0021) and ICM (β (SE) = −0.64 (0.29), p = 0.0366) (Table S2 and Fig. S3). Selection for MHC I diversity at PSS explained 33% and 28% of cross-species variation in CMR and ICM respectively.

Expansion of the MHCI gene family and cancer mortality risk
Allelic diversity within a gene can be shaped by copy number variation, where distinct alleles at individual loci enhance population-level diversity and enables single organisms to harbor diverse allelic repertoires. This is especially important for MHCI, exhibiting striking inter-specific variation in locus number27,28,40 from e.g., only three classical loci in humans to up to over a dozen loci in some rodents41. Given this variability, we explored whether the expansion of the MHC I gene complex shows a similar negative correlation with cancer mortality across species.

To quantify MHCI expansion, we estimated the number of protein coding genes in the MHCI gene family using Orthofinder39. Data were obtained for 127 and 119 mammalian species with available cancer mortality data (CMR and ICM respectively) and sequenced genomes. MHCI gene family size ranged from only one protein coding gene in the Southern three-banded armadillo (Tolypeutes matacus) to 44 in the White-faced saki (Pithecia pithecia), with a mean of 14.60 genes across taxa. We then analyzed non-zero cancer mortality (CMR: n = 106; ICM: n = 100) as a function of MHCI gene family size in phylogenetic logistic regressions. Notably, two rodents, the naked mole rat (Heterocephalus glaber) and the Patagonian mara (Dolichotis patagonum) emerged as high leverage outliers in our statistical analyses (Fig. 2). Both species exhibit an unusually low MHC I gene family size while also having exceptionally low cancer mortality. Including these species in our models altered model fits considerably, suggesting they may represent biologically distinct cases rather than generalizable trends. Indeed, the immune system of naked mole rats is known to be highly distinctive compared to other mammals (e.g., lack of natural killer cells, immune cell populations dominated by splenic immune cells: 56% compared to 8% in mice42, low frequency of T cells42,43). The immune system of the Patagonian mara remains largely unstudied, its unique placement in our analysis raised concerns about its comparability to other taxa in our database. Due to these factors, we excluded both species from our statistical models to ensure robust and generalizable conclusions. Nonetheless, these outliers are visualized in figures and their placement is discussed from an evolutionary perspective. Following these exclusions, our results indicated a steeply declining cancer mortality with the expansion of the MHCI gene family across mammals (CMR (β (SE) = −0.51 (0.14), p = 0.0004; ICM (β (SE) = −0.47 (0.17), p = 0.0068, Table S3 and Fig. 2). The coefficient of determination indicated that MHCI gene family size alone explains 37% and 15% of the cross-species variance in CMR and ICM respectively.

Confounding effects
Previous works have emphasized that cross-species variation in the number of MHC loci has co-evolved with body mass and longevity, at least in birds44. While CMR appears to be unrelated to body mass and adult life expectancy30 (but see ref. 45), these life-history traits might have the potential to confound associations detected between the complexity or diversity of MHCI and CMR. To explore this, we rerun the above models, while including adult body mass and adult life expectancy as covariates in phylogenetic logistic regressions exploring the link between cancer mortality measures and selection for MHCI diversity or MHC I gene family size. Model results containing these confounding effects provided highly congruent results to single-predictor models both in case of selection for MHC I diversity (Table S4), as well as in case of MHCI gene family size (Table S5). These models thus confirmed the negative association of the MHCI complexity and diversity with CMR across the mammalian phylogeny.

Our previous analyses on the cross-species variation in cancer risk across mammals has revealed an important role of the species’ diet30. We have shown that primarily consuming mammalian prey items significantly increases the risk of cancer mortality across mammalian species. Given the important role of mammal-based diet, we tested whether controlling for this effect in our analyses changes our conclusions regarding the role of MHCI in shaping cross-species variation in cancer mortality. The models indicated that the effect of selection for MHC I diversity (Table S6), as well as of MHCI gene family size (Table S7) on cancer mortality metrics remain robust and are not influenced by the consumption of mammalian prey across species.

The phylogenetic tree used in our analyses represents only one of the multiple hypotheses reflecting the potential evolutionary trajectory of the taxa included in these analyses. Therefore, to explore how sensitive our results are to the phylogenetic uncertainty, we run each model using 100 equally plausible phylogenies. Results on the association between CMR and selection for MHC I diversity or MHCI gene family size were highly consistent between alternative phylogenetic hypotheses (MHC diversity: Fig. 1c, d, MHC diversity at PSS: Fig. S3c, d, MHCI gene family size Fig. 2c, d).

Evolution of MHCI in function of life-history traits
In order to explore the extent to which selection on MHCI is influenced by the evolutionary history of the studied species, we explored the phylogenetic signal of the estimated MHCI traits. The phylogenetic inertia of selection for MHC I diversity was low when measured across all sites (n = 28, λ = 0.11, P = 0.6842, significance marks the difference from λ = 0), but showed moderate phylogenetic signal when measured across PSS (n = 28, λ = 0.55, P = 0.0066). Similarly, MHCI gene family size showed a very strong phylogenetic signal (n = 127, λ = 0.91, p < 0.0001, Fig. S4). We then explored whether selection on MHCI is more intense in large or long-lived taxa and if MHCI could be a factor contributing to the better cancer defense of such species. To do this, we gathered information on average adult body mass, adult life expectancy30 (i.e., average time lived by individuals following their sexual maturity) and 80% adult longevity (post-sexual-maturity age when 80% of the adult population is dead, representing a more reliable longevity metric than maximum lifespan46), and we tested the effects of these traits on MHCI characteristics. To test the correlation between selection for MHC I diversity or MHC gene family size and life-history traits, we built single-predictor phylogenetic regressions. We used MHCI characteristics as dependent variables and life-history traits as explanatory variables in subsequent models. Selection for MHC I diversity was unrelated to adult body mass, but tended to increase with adult life expectancy (R2 = 0.16), and increased significantly with increasing 80% adult longevity (R2 = 0.33, Fig. 3a–c and Table S8). On the contrary, selection for MHC I diversity at PSS showed no association with any of the three life history traits (Fig. 3d–f and Table S8). Similarly, MHCI gene family size was unrelated to body mass, adult life expectancy and 80% adult longevity across the studied mammals (Fig. 4a–c and Table S8).

Methods
Cancer mortality risk
We estimated CMR using data provided by the Zoological Information Management System (ZIMS), managed by Species360 (a non-profit organization custodian of zoo-aquaria data)61 following the methodology of Vincze et al.30. In short, we obtained a new data extraction of individual-level date of birth, date of death and cause of death (categorical data, including neoplasia as an option) from the husbandry module of ZIMS for animals that lived in zoos after the 1st of January 2010 (data extracted in June 2023). We excluded records of species that might have been subject to domestication (see ref. 30 for the list of excluded taxa), individuals below the age of sexual maturity (obtained from literature, see Source Data for species-specific resources) and all species where less than 20 individuals have relevant death information (RDI, information regarding the cause of death). Using the resulting database, we estimated CMR, as simple ratio between individuals who died of cancer divided by the total number of individuals with RDI in each species. Second, we estimated the cumulative ICM, a metric of CMR at 90% adult longevity, controlling for potential biases due to disregarding left-truncation (i.e., cancer before individuals enter the study) and right-censoring (individuals alive, thus with unknown fate at data extraction)30 using survival modeling. The database used to estimate CMR/ICM data for this study includes 80 species from 13 mammalian orders and contains information on 88,464 individuals, of which 34,128 were dead at the time of data extraction and 8503 had available RDI30.

Selection for MHC I diversity
In order to estimate selection for MHC diversity, we relied on publicly available resources from the NCBI GenBank Nucleotide Database. As sequence data in GenBank may originate from different sources (e.g., from different populations with varying demographic history) and may have been collected under different sampling schemes (the information on the number of genotyped individuals often remains unknown), it cannot be effectively used to quantify and compare standard measures of MHC diversity (e.g., allelic richness) across species. However, GenBank data provides a sample of species-specific MHC sequences, which can be used to detect molecular signature of historical selection at these genes and to reliably quantify variation in selection patterns among species.

In order to infer selection patterns at the mammalian MHC class I, we first downloaded GenBank sequences of MHC class I exons 2 and 3 (i.e., the two exons coding for the peptide-binding groove in MHCI) for species with available estimates of CMR (the last accession on June 2023). In three taxa with exceptionally high numbers of sequences (e.g., >3700 available sequences for Macaca fuscicularis), we downloaded 100 random sequences per exon. All sequences were aligned using Geneious software v. 10.0.5 (Biomatters Ltd., Auckland, New Zealand) and trimmed to our targeted exons. Identical allelic variants were removed from the alignments. The average length of alignments was 256.3 ± 2.9 [SE] bp and 261.9 ± 3.9 [SE] for exons 2 and 3, respectively, which covered 94–97% of full exon length.

To exclude the possible bias caused by the presence of non-classical MHC I sequences, we identified putative non-classical sequences based on the conservation of amino acid residues responsible for anchoring (docking) antigen peptide ends in the MHC peptide-binding groove. The anchor residues are highly conserved in classical MHC molecules across divergent vertebrate lineages62,63, but may experience radical amino-acid changes in non-classical genes. Here, we identified putative non-classical sequences based on the presence of at least one or two amino acid alterations at the anchor residues within exon 2 (n = 3 residues) and exon 3 (n = 7 residues), respectively. In total, we identified up to 12 putatively non-classical sequences per species (on average 4.4 ± 1.1 [SE] % of all sequences per species), and we excluded these from all further analyses. Species with fewer than four classical sequences per exon available were not considered in further analyses. After removal of putative non-classical sequences, consensus sequences from virtually all (96%) alignments better blast-matched alleles from classical (A-C) rather than non-classical (E-G) human MHCI loci. However, we were still not able to distinguish which classical loci our sequences were derived from. Such information could be only retrieved for a couple of model species, for which the detailed architecture of the MHC region is well resolved, and this was generally not feasible for non-model species. Here, we thus assume that the retrieved sequences were derived from a spectrum of different loci, possibly both paralogs and orthologs. Nevertheless, we believe this is unlikely to bias our estimates, as (1) we used the same methodology in every species; and (2) generally, we retrieved exons from GenBank (not from genomic resources) that were already annotated, avoiding potential issues with manual annotations.

Since recombination signal may bias estimates of nucleotide substitution rates (no unique tree topology can describe evolutionary history of the entire sequence64), we first aimed to identify recombination breakpoints and recombinant sequences in our data. For small alignments (<50 seq), we used the Genetic Algorithm for Recombination Detection (GARD, Datamonkey web server65) and partitioned all available sequences into recombinant fragments for further analyses. Since GARD performed poorly for larger alignments (i.e., 50–100 seq), recombinant sequences in these alignments were detected using the RDP4 software66. Briefly, recombination events were identified based on the consensus from at least two different recombination detection algorithms implemented in RDP4. On average, 15.8 ± 3.4 [SE] % of sequences were identified as recombinant and were removed, resulting in the final sample sizes lower than 100 sequences per exon in all the species.

The average final size of alignments (classical sequences with no recombinants) was 22.4 ± 3.0 [SE] sequences per exon per species, but ranged from 4 to 91 sequences. To account for this uneven sampling across species (i.e., variation in the number of available sequences), we applied a standardized subsampling procedure for species with >10 available sequences per exon. Specifically, we estimated selection for MHC I diversity using ten randomly selected (subsampled) sequences per species. Subsampling was conducted using fasta.sample function in the FastaUtils R package67 and iterated ten times per exon per species. In these cases, the estimates of nucleotide substitution rates (see below) were averaged across all iterations to obtain a representative measure of selection for each species.

Following the elimination of potential recombination-related bias and subsampling, we estimated the relative rate of non-synonymous (dN) and synonymous (dS) base pair substitutions using Fast Unconstrained Bayesian AppRoximation68 (FUBAR, Datamonkey) as a primary estimate of selection for MHC I diversity. Although the dN/dS ratio is traditionally used as a measure of positive selection across lineages (different alleles are advantageous and become fixed in independent lineages)36, pathogen-mediated balancing (diversifying) selection acting within lineages tend to leave the same hallmark signature at the MHC genes69. Hence, species-specific dN/dS ratios have commonly been used to estimate the strength of historical selection for MHC diversity (e.g., refs. 70–73, where dN/dS > 1 indicates diversifying selection (faster accumulation and maintenance of non-synonymous than synonymous substitutions), dN/dS < 1 indicates purifying selection (slower accumulation of non-synonymous than synonymous substitutions), while dN/dS = 1 indicates no signature of selection (non-synonymous and synonymous substitutions accumulate at similar rates).

The dN/dS estimates were first computed for all available codons within each alignment. Nonetheless, only about 25% of codons from the exons 2 and 3 code for the peptide-binding groove of the MHC I molecules37,38, which are truly responsible for direct interactions with antigens (so-called peptide binding region). Therefore, in order to confirm that the potentially detected associations are driven by functional (rather than structural) variation at the MHC, we re-extracted dN/dS ratios using only codons showing the strongest evidence of positive selection. Specifically, based on the values of Bayes factors, we identified the 20 most PSS across the MHCI sequences from all considered species (10 PSS per exon) and used the estimates of the dN/dS ratio across this set of sites as a proxy of selection for functionally relevant variation (selection for MHC I diversity at PSS). Selection for MHC I diversity at PSS was calculated similarly to the overall measures of selection for MHC I diversity (across all sites), using subsampling of 10 sequences, iterated 10 times, and averaged. While both dN and dS were initially calculated separately for exons 2 and 3, we averaged the values across the two exons to obtain a single species-specific estimate of dN/dS ratio at the MHC class I. Selection for MHC I diversity could be calculated for a total of 28 mammalian taxa with available cancer risk estimates.

Our analyses revealed relatively low dN/dS values, especially when all MHCI codons were considered (mean = 0.41, n = 28). This might be attributed to the fact that only a few codons of the peptide binding region are directly responsible for direct interactions with antigens and, thus, are under positive diversifying selection, while codons of the structural or cross-membrane part of MHC molecules remain generally invariant62. Also, pathogen or cancer-driven selection may selectively target only a fraction of codons of the peptide-binding region, while some other residues may be evolutionarily conserved due to their specific functions in peptide binding, for example, being responsible for binding peptide ends74. In fact, the mean number of PSS was only 3.62 and 3.39 at exon 2 and exon 3, respectively across the 28 species studied here (see values of PSS and NSS, respectively in Source Data). In contrast, we have evidence for an average of 5.70 and 8.00 negatively selected codons at the exon 2 and 3, respectively. Taking all this into account, there is much stronger signature of negative than positive selection across the entire exon 2/3 sequence of the mammalian MHC-I, which results in low overall dN/dS ratios. At the same time, our data shows strong evidence of positive diversifying selection at some codons and the mean dN/dS of the 20 most PSS per species is 2.99 across the taxa explored here (see Table Source Data), providing support for strong diversifying selection acting at a fraction of MHC-I codons.

MHCI gene family size
Protein-coding sequences from wild mammals with known CMR/ICM estimates (n = 127) were obtained from genomes annotated with TOGA75 (Tool to infer Orthologs from Genome Alignments) via the Senckenberg data repository (https://genome.senckenberg.de/download/TOGA/). TOGA is a genome annotation software which uses similarities of intronic and intergenic regions to identify orthology components by comparison to a well-annotated reference genome (here, human, GRCh38/hg38). We used TOGA-annotated genomes because of their enhanced annotation completeness, especially of conserved mammalian genes75 and annotation consistency, improving comparability among the genomes. If multiple assemblies were available for a single species, only the one with the highest contig N50 score was retained. We also downloaded the human proteome from the Ensembl Genome Browser76 to identify genes of interest in the target species among the orthogroups. Exclusively, the longest transcript variant of each gene per genome was used for downstream analysis. To assign the protein sequences of each species to orthogroups, we used OrthoFinder v2.5.539 with default settings. The sequence similarity search tool DIAMOND77 was used to perform a comprehensive comparison of the protein sequences. Subsequently, homologous genes were clustered into orthogroups defined as sets of genes originating from a single gene in the common ancestor of the analyzed species39. Orthofinder assigned the protein sequences into 22,327 orthogroups, including one orthogroup (OG0000011) containing the human HLA class I genes (AZGP1, HLA-A, HLA-B, HLA-C, HLA-E, and HLA-G). We used the latter orthogroup as an estimate of MHCI gene family size and we used this variable as predictor of CMR/ICM in phylogenetic logistic regressions.

Phylogeny and statistical analyses
To account for the statistical non-independence of species due to common descent, we obtained a sample of 1000, equally plausible phylogenetic trees, from the posterior distribution published by Upham et al.78, covering 5911 species. We then obtained a rooted consensus tree using the “sumtrees” Python library79. To explore the effect of phylogenetic uncertainty, all models were run using 100 equally likely phylogenetic hypotheses and the results were checked for consistency across these. All statistical analyses were carried out in R statistical and programming environment, R v.4.3.180, using the R function gls from the package “nlme” and Brownian motion evolutionary models. Partial coefficients of determinations were calculated using the function “R2.pred” from the R package rr281, based on models presented in Tables S1 and S3. All statistical tests performed are two-tailed. No correction for multiple testing was performed. Data and R code needed to reproduce the analysis are publicly available (see “code availability” section below). Animal silhouettes used to visually represent mammalian orders were downloaded from PhyloPic (http://www.phylopic.org).

To test the association between CMR and MHCI, we built PGLS models. We used cancer mortality risk (CMR and ICM) as dependent variables and selection for MHC diversity (overall or at PSS), or MHCI gene family size as explanatory variables. Given that the precision of estimates of CMR is contingent on the underlying sample sizes, we controlled for this biasing effect by including the log number of individuals with RDI per species as model weight in all models of CMR.

Confounding effects
To test the correlation between life-history traits and MHC diversity, we retrieved information on body mass from ZIMS and from literature resources30 and we estimated species-specific adult life expectancy and 80% adult longevity using the individual husbandry database obtained to estimate CMR/ICM. We calculated adult life expectancy as a simple average of post-sexual-maturity lifespan of dead individuals. For 80% adult longevity we calculated age-specific survival rates using the Kaplan–Meier procedure (survfit function in the R package survival82). Individuals older than their age at sexual maturity at the beginning of data registration were left-truncated at their age at this date; individuals reaching sexual maturity after this date were left-truncated at their age at sexual maturity. Individuals still alive at the time of data extraction were considered right-censored (samples per species varied from 94 to 8460 individuals), while known fate individuals were assigned as dead (n = 34,128), irrespective of whether RDI was available for the individual or not.

We also performed multiple sensitivity analyses, to explore the robustness of our results. (1) To test the potentially confounding effects of body mass and life expectancy on the relationship between CMR and selection for MHC I diversity, we introduced these two variables as covariates in PGLS models with CMR as dependent and MHC diversity as explanatory variables. (2) To check the robustness of our results to the phylogenetic hypothesis, we repeated single-predictor phylogenetic logistic regressions using a random sample of 100 trees provided by Upham et al.79. (3) Given that we previously have demonstrated that the consumption of mammalian prey items significantly increases the risk of cancer across mammals, we aimed to test whether this effect might confound our results. To do so, we obtained information on the diet of each species83 and we re-categorized mammal-based diet at two levels: never/rarely consumed or representing the primary/secondary food item of the species (see also ref. 1). Mammal consumption was introduced as a categorical explanatory variable to phylogenetic logistic regressions between cancer mortality and selection of MHCI diversity or MHCI gene family size.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, the authors found a negative correlation between cancer risks and “selection for MHC I diversity” across\nmammal species (but not at the PBR region, which is surprising). This result supports the cancer immunosurveillance\nhypothesis. The manuscript is well written and the results are highly significant to the field. However, I am quite surprised by\nthe measure of MHC diversity (or “selection for MHC diversity”), and the low values detected. This needs to be explained in\nmore details in the manuscript.\nThe dn/ds is not a common measure of MHC diversity. The rational for using this measure is relatively well explained in the\nmethod sections, but I would suggest to add some details explaining why you use dn/ds instead of gene copy number or\nallelic diversity in the main text also.\nLine 208 : It seems surprising to pool classical MHC I and non-classical MHCI genes together in the analysis. I would\nsuggest to analyse classical MHC I genes only.\nLine 222: Why did you focus on viruses ?\nLine 333: “Nonetheless, most species do not exhibit extremely high MHC I diversities”. Is this statement based on your\nresults or on the literature? If it is based on your measures: is dn/ds a measure of diversity or selection for diversity ? If it is\nbased on the literature, please add a reference.\nLine 343: This result is very surprising. Many studies show that dn/ds ratio is higher at the PBR than outside the PBR region,\nwhich kind of confirms the identification of the PBR region. In your data, are dn/ds ratios higher for the PBR than outside the\nPBR region ?\nLine 355-356: Please add a reference.\nLine 406: Good quality genomes of many mammal species are published at NCBI. A measure of gene copy number can\nprobably be retrieved from many of these species (see He et al. 2024 for an example in 50 ruminant species). Please\nexplain why you did not use this method to get information about gene copy numbers.\nLine 414: “n = 3 species”. Do you mean you have 3 species with good quality data ?\nIn your data, selection for MHC diversity ranges between 0.25 – 0.7 indicating negative selection. There is no indication of\ndiversifying, positive selection. So is MHC diversity selected ? Please explain.\nI am also quite surprised by the dn/ds values. They seem very low compared to the values found in the literature (usually > 1\nindicating diversifying, positive selection). A comparison with the values found in the literature is needed.\nMHC class II presentation has also a central role in tumor evolution. A similar study with MHC class II genes may also be\nequally interesting.\nFig S2 is not in the pdf file. So we cannot see the dn/ds values for the PBR region. Were these values higher than the dn/ds\nratios for the whole region ?\nThe structure of the manuscript is difficult to follow. A “Results” heading may help.\nIn the main text, some reference numbers are indicated in superscript, while other are not.", "reviewer_2": "This manuscript presents a comparative study in mammals, examining the evolutionary links between MHC class I diversity\nand cancer risk. This study significantly contributes to the emerging knowledge on cancer ecology and evolution, providing\nimportant insights into the evolution of cancer resistance. The design and methods are sound, including testing relevant\nconfounding variables or potential mediators of the observed links and examining robustness to phylogenetic uncertainty.\nThe authors interpret their results with a high level of knowledge and care, and provide a well-grounded discussion without\noverextending their conclusions. The manuscript is written in a clear, accessible and concise way, making the text engaging\nand easy to follow. Overall, this is a high-quality manuscript that will certainly appeal to a wide audience. I would endorse\nthe manuscript publication with only a few minor comments below.\nThe fact that the PGLS models used to model the association between cancer risk and MHC diversity are logistic regression\nmodels is only mentioned in the section of the Methods describing sensitivity analyses and in captions of supplementary\ntables. I recommend to clearly state this in the main section about statistical analysis (‘Phylogeny and statistical analyses’).\nIf I understand correctly, what you call the zero-inflated models (e.g., line 530), are in fact two separate models. If that is the\ncase, your approach should not be called zero-inflated models in my opinion because it differs from the actual zero-inflated\nmodels in at least two aspects: 1) the (structural) zeros and the non-zero component are not estimated simultaneously, and\n2) your approach does not model sampling zeros. If you really fitted zero-inflated models, please specify the package used\nbecause zero-inflated models are not implemented in ‘nlme’ package as far as I am aware.", "reviewer_3": "Report on Immunological surveillance against cancer across mammals\nIn this manuscript, the authors combine a previously published dataset on cancer incidence and mortality in mammals with\ndata on MHC genetic variation from the NCBI nucleotide datase in order to investigate whether cancer-related mortality is\nassociated with MHC class I diversity. The authors thus aim to test a long-standing hypothesis that the increasing risk for\ntumors in larger and longer-lived species was an important factor in the evolution of adaptive immunity.\nWhile the hypothesis is very intriguing and the proposed mechanisms are plausible, such cross-species comparisons are\ninherently difficult, because of the difficulty of unbiased and systematic sampling.\nThe previously published dataset on cancer risk among 191 species is very comprehensive and uniquely suitable for this\ntest. Unfortunately, the data on MHC class I diversity from public databases is much less comprehensive, leading to a final\ndataset with only about two dozen species. MHC-I diversity is here measured as the ratio of non-synonymous over\nsynonymous mutations per species, which is an accepted measure for positive selection.\nOverall, the authors find a negative correlation between cancer risk and positive selection in MHC-I sequences and interpret\nthis as proof for the original cancer immunosurveillance hypothesis.\nWhile the authors employ several sanity checks to increase the confidence in their results, several uncertainties remain. I am\ndetailing several conceptual and analytical issues below. Overall, I am not fully convinced that the reported association is\nrobust and meaningful.\nThe most immediate critical problem is the comparability of the MHC-I sequence data among species. The authors admit\nthat the number of sequences obtained per species in the NCBI database is not a good proxy for MHC-I diversity of a\nspecies, as it is inherently biased by the number of individuals sampled. They argue that the chosen dN/dS ratio is more\nrobust. However, as can be tested from the data in the supplementary table, the authors‘ measure of positive selection is\npositively correlated with the number of sequences obtained per species, which is not surprising.\nThe main correlation and thus the main result is thus dependent on the species-specific sampling of MHC sequencies in the\npublic database, and this is not at all standardized or controlled. It is thus unclear what the observed correlation is based on,\nand whether it reflects a biological phenomenon at all.\nConceptually, it is also not clear to me how tumor-based mortality would select for a high allelic diversity in the\npopulation/species. As the authors outline in their introduction, cancer immunosurveillance requires the detection of as many\nneoantigens as possible, thus selecting for MHC heterozygosity. But heterozygosity can be achieved by just two alleles, and\nit doesn’t require the large number of alleles observed in most species. If anything, cancer immunosurveillance should select\nfor more MHC gene copies (with different alleles), as this would allow presenting more different tumor antigens. However,\nthe authors do not see a correlation between MHC gene number and cancer risk.\nGenerally, the authors fail to discuss their findings in the light of what is know about MHC diversity in other, non-mammalian\nspecies. Many songbirds show a vast allelic diversity in their MHC genes, but are not known for either high cancer rates nor\nlarge bodies or particularly long life. Similarly, many fish species show high MHC diversity, compared to mammals, although\nthey liver shorter and are mostly smaller. How do the authors reconcile their findings with these observations?\nIt is not always clear what the authors see as cause and what as consequence. For instance, the sentence in L123 sounds\nas if MHC diversity is expected to drive the risk for cancer, but the cancer immunosurveillance hypothesis would expect the\nopposite.\nIn their original paper on cancer risk across 191 mammalian species, the authors identify several life-history traits that\ncorrelate with cancer risk and mortality. However, they are here ignored. Instead, these traits should here be used as\ncovariates, in order to control for their effect.\nGenerally, there is very little detail on how the MHC data was obtained from the NCBi database, nor can I find information on\nsource about the MHC gene number per species. This information needs to be made available.\nPredicted viral richness is not a good proxy for pathogen-mediated selection. Prevalence and transmissibility would be more\nsuitable. Using this data here is thus of little value and can not be taken as proof that MHC diversity is not driven by\npathogen-mediated selection."}, "author_response": "RESPONSE TO REVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nIn this manuscript, the authors found a negative correlation between cancer risks and “selection\nfor MHC I diversity” across mammal species (but not at the PBR region, which is surprising).\nThis result supports the cancer immunosurveillance hypothesis. The manuscript is well written\nand the results are highly significant to the field. However, I am quite surprised by the measure\nof MHC diversity (or “selection for MHC diversity”), and the low values detected. This needs to\nbe explained in more details in the manuscript.\nWe believe that the low overall dN/dS values may be attributed to the fact that only a few PBR\ncodons are directly responsible for direct interactions with antigens (peptides) and, thus are under\npositive diversifying selection. It has been long known that some PBR residues remain generally\ninvariant due to their specific functions in peptide binding, for example, being responsible for\nbinding peptide ends (see Kaufman et al. 1992). At the same time, pathogen (or cancer) driven\nselection may selectively target only a fraction of PBR codons. Apparently, at the mammalian\nMHC-I this fraction is relatively small. We have now added the number of positively and negatively\nselected sites for both exon 2 and exon 3 to our dataset that we aim to publish along with our\nmanuscript. We also discuss the relatively low dN/dS values detected when all codons are\nconsidered. We show that the mean number of positively selected sites (PSS) is only 3.61 and 3.39\nat exon 2 and exon 3, respectively (see values of PSS and NSS, respectively in Table S9). In\ncontrast, we have evidence for an average of 5.70 and 8.00 negatively selected codons at the exon\n2 and 3, respectively. Taking all this into account, there is much stronger signature of negative\nthan positive selection across the entire exon 2/3 sequence of the mammalian MHC-I, which\nresults in low overall dN/dS ratios. At the same time, our data shows a strong evidence of positive\ndiversifying selection at some codons and the mean dN/dS of the 20 most positively selected sites\nper species is 2.27 across the taxa explored here (see Table S9), providing support for strong\ndiversifying selection acting at a fraction of MHC-I codons. Please see line 519-530.\nKaufman, J., Andersen, R., Avila, D., Engberg, J., Lambris, J., Salomonsen, J., ... &\nSkjødt, K. (1992). Different features of the MHC class I heterodimer have evolved at\ndifferent rates. Chicken BF and beta 2-microglobulin sequences reveal invariant surface\nresidues. Journal of immunology (Baltimore, Md.: 1950), 148(5), 1532-1546.\nThe dn/ds is not a common measure of MHC diversity. The rational for using this measure is\nrelatively well explained in the method sections, but I would suggest to add some details\nexplaining why you use dn/ds instead of gene copy number or allelic diversity in the main text\nalso.\nThank you for these remarks, we have added a short argument for this choice in the main\ntext too. See line 156-158. Please note that we have now also estimated the number of\nprotein coding genes within the MHC I gene family as an alternative measure of the\ndiversity and complexity of the MHC I complex across species.\nLine 208 : It seems surprising to pool classical MHC I and non-classical MHCI genes together in\nthe analysis. I would suggest to analyse classical MHC I genes only.\nThank you for this comment. We perfectly agree with this and we have now rerun all\nanalyses only relying on classical MHC-I loci. Our results and conclusions remained\nqualitatively unchanged. See lines 445-457.\nLine 222: Why did you focus on viruses?\nViruses are obligate intracellular organisms and therefore relevant to MHC-I\nimmunosurvillance. Yet, we recognize that viruses are not the only intracellular\npathogens (e.g. intracellular bacteria, fungi, protozoa, haemosporodians also exist) and\nthe lack of effect of viral pathogen pressure on the link between MHC I and cancer\nmortality is not definitive proof of the lack of overall effect of intracellular pathogen\npressure. For this reason and following the criticisms of Reviewer 3 we decided to\nremove this analysis from the manuscript and only discuss this potential caveat in the\ndiscussion.\nLine 333: “Nonetheless, most species do not exhibit extremely high MHC I diversities”. Is this\nstatement based on your results or on the literature? If it is based on your measures: is dn/ds a\nmeasure of diversity or selection for diversity ? If it is based on the literature, please add a\nreference.\nThank you for spotting the lack of information here. This statement is indeed based on\nthe literature (e.g. Migalska et al. 2019). High MHC diversity has been highlighted\nrepeatedly as maladaptive, especially since it renders more T-cell receptors self-\nreactive/it increases the probability of T cell self-reactivity, leading through negative\nthymic selection to a drastic reduction of the T-cell receptor repertoire that would be\ncapable of dealing with pathogens. We now provide a reference for our statement. See\nline 334.\nMigalska, Alvaro, Sebastian, Jacek, Radwan (2019). 1. Major histocompatibility complex\nclass I diversity limits the repertoire of T cell receptors. Proceedings of the National\nAcademy of Sciences of the United States of America, doi: 10.1073/PNAS.1807864116\nLine 343: This result is very surprising. Many studies show that dn/ds ratio is higher at the PBR\nthan outside the PBR region, which kind of confirms the identification of the PBR region. In your\ndata, are dn/ds ratios higher for the PBR than outside the PBR region ?\nYes, the MHC-I diversities are considerably higher in the PBR than overall (please see Fig\nS3, Table S9, but also our response to Reviewer 1’s first comment). As the Reviewer\npoints out, this confirms that the codons we selected are indeed coding the peptide-\nbinding region that should in theory be under the highest selection. Yet, we cannot be\nconfident that the selected codons are all in the PBR, or if all PBR codons have been\ncorrectly identified and selected. Overall, we did find the association between cancer risk\nand dN/dS at PBR weaker, which is contrary to all our expectations. While this result was\ndiscouraging, the stability of the results using dN/dS across the entire MHC-I region, as\nwell as the newly presented analyses on the association between MHC-I gene family size\nand cancer risk made us mostly doubt the proper identification of the PBR regions\nacross this wide taxonomic scale.\nLine 355-356: Please add a reference.\nGiven that we now estimated the number of protein-coding genes within the MHC I gene\nfamily, which provides a similar measure to the number of MHC I loci, we decided to\nremove the analyses, including the number of MHC I loci numbers that we previously\nobtained from primary literature resources. The primary reason is that the latter could\nonly be obtained for a very limited number of species and thus the analyses provided no\nconclusive evidence on the presence or absence of a role of increased loci numbers in\ncancer defence. We thus removed the argument mentioned by the Reviewer from the\ntext.\nLine 406: Good quality genomes of many mammal species are published at NCBI. A measure\nof gene copy number can probably be retrieved from many of these species (see He et al. 2024\nfor an example in 50 ruminant species). Please explain why you did not use this method to get\ninformation about gene copy numbers.\nThank you for this comment. We have now replicated our analyses using high-quality\ngenomes that we downloaded from NCBI. We used Orhtofinder (v2.2.5, Emms & Kelly\n2019), a computational tool designed for inferring orthogroups (i.e., groups of genes that\nshare a common ancestor) across different species. Using this tool we could quantify the\nnumber of protein-coding genes in the MHC-I gene family across all species where the\ngenome was available. We identified a single orthogroup, in which most MHC-I genes\nwere clustered (including MHC-I A, B, C, E, and G), thus providing a precise estimate of\nthe number of protein-coding genes in this gene family across species. Using the size of\nthis orthogroup as a predictor, we could demonstrate that the expansion of the MHC-I\ngene family is coupled with a significantly reduced cancer mortality risk across\nmammals. We have now incorporated this analysis in our manuscript, hoping to\nconvince readers who might have similar scepticism about our previous analyses.\nPlease see the detailed description of the Orthofinder protocol, model fitting and results\nin the revised ms.\nFigure Redacted\nDavid, M., Emms., Steven, L., Kelly. (2019). 1. OrthoFinder: phylogenetic orthology\ninference for comparative genomics. Genome Biology, doi: 10.1186/S13059-019-1832-Y\nLine 414: “n = 3 species”. Do you mean you have 3 species with good quality data ?\nIn your data, selection for MHC diversity ranges between 0.25 – 0.7 indicating negative\nselection. There is no indication of diversifying, positive selection. So is MHC diversity\nselected ? Please explain.\nThank you for spotting this. We have now rephrased this sentence to “taxa with\nexceptionally high numbers of sequences“.\nI am also quite surprised by the dn/ds values. They seem very low compared to the values\nfound in the literature (usually > 1 indicating diversifying, positive selection). A comparison with\nthe values found in the literature is needed.\nAs we explained above, it seems that only a fraction of all MHC I codons are under\npositive selection in mammals, leading to relatively low dN/dS values when all codons\nare considered. Please see our discussion on this issue on Lines 515-530.\nMHC class II presentation has also a central role in tumor evolution. A similar study with MHC\nclass II genes may also be equally interesting.\nThank you for this remark and for the interesting idea. While we want to maintain the\nfocus of our current work on MHC-I, as we originally planned, we might indeed repeat the\nsame line of analyses on other immune genes, including MHC II in the future.\nFig S2 is not in the pdf file. So we cannot see the dn/ds values for the PBR region. Were these\nvalues higher than the dn/ds ratios for the whole region ?\nWe’re truly sorry for this mistake. We have now added figure S2 to the manuscript file.\nBesides, we have also added a figure (S3) to the manuscript showing the association\nbetween overall dN/dS and dN/dS at putative PBR sites. As the Reviewer can see, the\ndN/dS values at PBR are considerably and consistently higher in all species than overall\ndN/dS.\nThe structure of the manuscript is difficult to follow. A “Results” heading may help.\nThank you for this remark. We have now added a “Results” heading and tried to clarify\nthe structure.\nIn the main text, some reference numbers are indicated in superscript, while other are not.\nThank you, we have now corrected this formatting issue.\n========================================================================\n======\nReviewer #2 (Remarks to the Author):\nThis manuscript presents a comparative study in mammals, examining the evolutionary links\nbetween MHC class I diversity and cancer risk. This study significantly contributes to the\nemerging knowledge on cancer ecology and evolution, providing important insights into the\nevolution of cancer resistance. The design and methods are sound, including testing relevant\nconfounding variables or potential mediators of the observed links and examining robustness to\nphylogenetic uncertainty. The authors interpret their results with a high level of knowledge and\ncare, and provide a well-grounded discussion without overextending their conclusions. The\nmanuscript is written in a clear, accessible and concise way, making the text engaging and easy\nto follow. Overall, this is a high-quality manuscript that will certainly appeal to a wide audience. I\nwould endorse the manuscript publication with only a few minor comments below.\nWe are grateful to Reviewer2 for the positive feedback on our work and for the remarks\nthat helped us improve our article.\nThe fact that the PGLS models used to model the association between cancer risk and MHC\ndiversity are logistic regression models is only mentioned in the section of the Methods\ndescribing sensitivity analyses and in captions of supplementary tables. I recommend to clearly\nstate this in the main section about statistical analysis (‘Phylogeny and statistical analyses’).\nThank you for spotting this. We have now stated in the main section that we used logistic\nregressions. See lines 185-186.\nIf I understand correctly, what you call the zero-inflated models (e.g., line 530), are in fact two\nseparate models. If that is the case, your approach should not be called zero-inflated models in\nmy opinion because it differs from the actual zero-inflated models in at least two aspects: 1) the\n(structural) zeros and the non-zero component are not estimated simultaneously, and 2) your\napproach does not model sampling zeros. If you really fitted zero-inflated models, please\nspecify the package used because zero-inflated models are not implemented in ‘nlme’ package\nas far as I am aware.\nThank you for this remark, we fully agree that the text was imprecise. Please note\nhowever that these models were now removed from the manuscript following the\nrecommendation of Reviewer 1.\n========================================================================\nReviewer #3 (Remarks to the Author):\nReport on Immunological surveillance against cancer across mammals\nIn this manuscript, the authors combine a previously published dataset on cancer incidence and\nmortality in mammals with data on MHC genetic variation from the NCBI nucleotide datase in\norder to investigate whether cancer-related mortality is associated with MHC class I diversity.\nThe authors thus aim to test a long-standing hypothesis that the increasing risk for tumors in\nlarger and longer-lived species was an important factor in the evolution of adaptive immunity.\nWhile the hypothesis is very intriguing and the proposed mechanisms are plausible, such cross-\nspecies comparisons are inherently difficult, because of the difficulty of unbiased and systematic\nsampling.\nThe previously published dataset on cancer risk among 191 species is very comprehensive and\nuniquely suitable for this test. Unfortunately, the data on MHC class I diversity from public\ndatabases is much less comprehensive, leading to a final dataset with only about two dozen\nspecies. MHC-I diversity is here measured as the ratio of non-synonymous over synonymous\nmutations per species, which is an accepted measure for positive selection.\nOverall, the authors find a negative correlation between cancer risk and positive selection in\nMHC-I sequences and interpret this as proof for the original cancer immunosurveillance\nhypothesis.\nWhile the authors employ several sanity checks to increase the confidence in their results,\nseveral uncertainties remain. I am detailing several conceptual and analytical issues below.\nOverall, I am not fully convinced that the reported association is robust and meaningful.\nWe are very grateful for Reviewer 3 for the time taken to carefully read our work and for\nthe critical reception of our work. Given the low sample size we could previously obtain\non dN/dS ratios to characterize MHC I, we have now therefore taken extra steps to\nincrease the statistical power of our analyses and to validate our results. Most\nimportantly, we performed a re-extraction and re-estimation of the cancer mortality\nestimates using more individual animal records from more species, in order to use the\nmost extensive database possible on cancer mortality and to increase the overlap of the\nMHC I and cancer databases. This has increased considerably the statistical power of\nour analyses . Second, we performed additional analyses based on genomic information,\nwhere we quantified the number of protein-coding genes in the MHC-I complex of\nspecies with accessible cancer mortality estimates. Analysis using this alternative\nmeasure of MHC I complexity we could show that the expansion of the MHC-I gene family\nis strongly, negatively linked to cancer mortality risk across 58 mammalian species and\ntherefore strongly corroborates our previous results. Please see our more specific\nresponses below. We hope these new analyses will convince the reviewer that the\nreported associations are indeed robust and meaningful.\nThe most immediate critical problem is the comparability of the MHC-I sequence data among\nspecies. The authors admit that the number of sequences obtained per species in the NCBI\ndatabase is not a good proxy for MHC-I diversity of a species, as it is inherently biased by the\nnumber of individuals sampled. They argue that the chosen dN/dS ratio is more robust.\nHowever, as can be tested from the data in the supplementary table, the authors‘ measure of\npositive selection is positively correlated with the number of sequences obtained per species,\nwhich is not surprising.\nThe main correlation and thus the main result is thus dependent on the species-specific\nsampling of MHC sequencies in the public database, and this is not at all standardized or\ncontrolled. It is thus unclear what the observed correlation is based on, and whether it reflects a\nbiological phenomenon at all.\nWe appreciate the Reviewer’s comment on the potential biases introduced by the number\nof sequences available in GenBank. However, we would like to clarify why a positive\ncorrelation between dN/dS ratios and the number of sequences retrieved per species is\nexpected due to underlying biological reasons rather than a mere artifact of sampling\nbias.\nTo illustrate this, let’s consider two hypothetical populations of 1000 individuals each,\nwith equal sampling effort (e.g., 10% of individuals genotyped). Suppose these\npopulations differ in MHC polymorphism, with total allelic richness of 200 and 50 alleles,\nrespectively. If the number of MHC alleles retrieved is proportional to the sampling effort,\nwe expect to recover 20 and 5 alleles from each population, respectively. As these\nsequences are deposited in GenBank, the difference in sequence numbers between the\ntwo populations is driven primarily by differences in MHC polymorphism rather than just\nsampling effort. Since MHC polymorphism is a direct consequence of the strength of\ndiversifying selection (dN/dS), it follows that dN/dS ratios will be positively correlated\nwith the number of sequences deposited in GenBank.\nOf course, sampling effort does contribute to the number of sequences available,\nmeaning that the number of sequences in GenBank reflects both the underlying\npolymorphism level (driven by selection) and the extent of sampling. However, including\nthe number of sequences as a linear covariate in the analysis would risk removing\ngenuine biological variation in MHC diversity that we aim to capture. We hope this\nexplanation clarifies our reasoning and why we believe our approach remains valid in\ncapturing meaningful biological patterns. Moreover, the newly implemented genomic\nanalyses further supports our claim on the importance of immunological surveillance in\ncancer defense, without the problem of sampling effort, strengthening our conclusions\nbased on dN/dS.\nConceptually, it is also not clear to me how tumor-based mortality would select for a high allelic\ndiversity in the population/species. As the authors outline in their introduction, cancer\nimmunosurveillance requires the detection of as many neoantigens as possible, thus selecting\nfor MHC heterozygosity. But heterozygosity can be achieved by just two alleles, and it doesn’t\nrequire the large number of alleles observed in most species. If anything, cancer\nimmunosurveillance should select for more MHC gene copies (with different alleles), as this\nwould allow presenting more different tumor antigens. However, the authors do not see a\ncorrelation between MHC gene number and cancer risk.\nA high MHC diversity (polymorphism) at the population level is considered to be useful\nagainst new pathogens (pandemics), because it ensures that at least some individuals\nwill be able to mount a protective immune response by presenting antigens to CD8 T\ncells, and thereby survive. But cancer is not a transmittable disease, so it is indeed a\nvalid question: how can a large number of MHC I alleles in a species (nearly 30,000\ncurrently described in humans, Robinson et al. 2024) protect better against cancer at the\nlevel of the individual? MHC heterozygosity is a plausible explanation. In fact, even if in\nsome cases heterozygosity can be achieved by just two alleles, on a population level this\nis insufficient. If there were only two MHC I alleles (with each a 50% frequency in the\npopulation), half of the individuals would be homozygotes. In contrast, a high number of\ndifferent MHC I alleles (and haplotypes) reduces the frequency of homozygotes\nsignificantly in the population (unless there is inbreeding). Moreover, the new genomic\nanalyses we present does indeed highlight that expansion of the MHC-I gene family is\nassociated with reduced risk of cancer, providing further support for the cancer\nimmunosurveillance hypothesis (see Fig 2, Table S3, S5).\nRobinson, J., Barker, D. J., & Marsh, S. G. (2024). 25 years of the IPD IMGT/HLA D\natabase. Hla, 103(6), e15549.\n‐\nGenerally, the authors fail to discuss their findings in the light of what is know about MHC\ndiversity in other, non-mammalian species. Many songbirds show a vast allelic diversity in their\nMHC genes, but are not known for either high cancer rates nor large bodies or particularly long\nlife. Similarly, many fish species show high MHC diversity, compared to mammals, although\nthey liver shorter and are mostly smaller. How do the authors reconcile their findings with these\nobservations?\nThis is a very interesting question and idea. Indeed it would be very interesting to\ncompare classes of vertebrates in regards to the general characteristics of their MHC I\nand cancer risk. Yet we were unable to find any work that would provide a\ncomprehensive comparison of MHC I complexity among vertebrate classes, or even just\nbirds and mammals.\nRegarding the Reviewer’s comments, some passerines indeed appear to have\nexceptionally high numbers of classical MHC I loci. We also know that birds generally\nface an extremely low risk of cancer compared to mammals (see e.g. Compton et al.\n2025, which is also supported by our unpublished data on ca. 400 species of birds),\nwhich would be in line with the conclusions of our analyses in mammals. Moreover, birds\ndo have very long lifespans compared to size-matched mammals (see e.g. Costantini\n2008). While this is interesting, we do need to keep in mind that not all birds have highly\nexpanded MHC I complexes, and most basal birds, like Galliformes harbour only one or\ntwo classical MHC I loci (Minias et al. 2016). Moreover, in classes of vertebrates outside\nbirds and mammals, the risk of cancer is extremely little known, so discussions\nregarding its link to the structure of MHC I within the taxonomic class would be very\ndifficult. Overall, we feel that putting our results in a broader context is very appealing,\nbut despite our best efforts, we feel that without former analyses and data to back up our\nclaim, the best we can do is to speculate, that we prefer to avoid. We remain open to\nsuggestions regarding this discussion given the above limitations.\nCompton, C., & Compton, C. (2020). A Short History of Cancer: How Did We Get Here?.\nCancer: The Enemy from Within: A Comprehensive Textbook of Cancer’s Causes,\nComplexities and Consequences, 83-107.\nCostantini, D. (2008). Oxidative stress in ecology and evolution: lessons from avian\nstudies. Ecology letters, 11(11), 1238-1251.\nMinias, P., Bateson, Z. W., Whittingham, L. A., Johnson, J. A., Oyler-McCance, S., & Dunn,\nP. O. (2016). Contrasting evolutionary histories of MHC class I and class II loci in\ngrouse—effects of selection and gene conversion. Heredity, 116(5), 466-476.\nIt is not always clear what the authors see as cause and what as consequence. For instance,\nthe sentence in L123 sounds as if MHC diversity is expected to drive the risk for cancer, but the\ncancer immunosurveillance hypothesis would expect the opposite.\nIndeed, it’s rather difficult to establish causality here. It is likely that elevated risk of\ncancer increases selective pressure on MHC I, so high MHC I diversity can be the\nconsequence of high cancer risk. On the other hand, efficient MHC I based\nimmunosurveillance can decrease cancer risk, so in this sense low cancer risk is a\nconsequence of high MHC I complexity/diversity. We tried to clarify this sentence now.\nIn their original paper on cancer risk across 191 mammalian species, the authors identify\nseveral life-history traits that correlate with cancer risk and mortality. However, they are here\nignored. Instead, these traits should here be used as covariates, in order to control for their\neffect.\nThank you for this remark. Indeed, in our previous article, we have shown that diet,\nespecially the consumption of mammalian prey, is associated with an elevated risk of\ncancer across mammals. This might be highly relevant in the case of the MHC gene\ncomplexes, as it has been hypothesized that high levels of observed polymorphism at\nMHC loci are maintained by natural selection acting to preserve alleles for pathogen\nresistance, which might be highly relevant in carnivorous animals. To explore the\npotential confounding effect of diet we have tried to introduce the consumption of\nmammals (categorical variable coded as 0, if mammals don’t or rarely occur in the diet of\nthe species or 1, if mammals constitute the primary or secondary food items of the\nspecies, just like we did on our previous article; please see data source and\nmethodology in Vincze et al. 2022). We have included mammalian diet as a predictor of\ncancer risk besides dN/dS ratios as well as the size of the MHC-I gene family. In both\ncases, mammal-consuming diet seemed to be unrelated to the risk of cancer (no\nsignificant effect was detected) and the strength of the association between selection for\nMHC-I diversity or MHC-I gene family sizes with cancer risk remained qualitatively\nunchanged. We now present two of the simpler models in the appendix and discuss this\nresult shortly in the paper (see lines 245-250). Moreover, to the best of our knowledge, no\nother life-history trait was shown to explain a significant portion of the cross-species\nvariation in cancer risk, therefore we found no other variables that might be important to\ncontrol for here. Yet we remain open to further suggestions.\nGenerally, there is very little detail on how the MHC data was obtained from the NCBi database,\nnor can I find information on source about the MHC gene number per species. This information\nneeds to be made available.\nThank you very much for spotting the lack of information here. We removed the number\nof classical MHC I loci data from the analyses because its availability was limited to 20\nspecies and the number of protein-coding genes within the MHC I gene provides a highly\ncongruent measure to this with a much larger sample size. Nonetheless, following your\ncriticism, we now provide full sequence alignments in fasta files that we used to\ncalculate dN/dS ratios, and we provide the accession number of all genomes analyzed to\nobtain the number of protein-coding genes in the MHC I gene family, to increase the\nreproducibility of our study.\nPredicted viral richness is not a good proxy for pathogen-mediated selection. Prevalence and\ntransmissibility would be more suitable. Using this data here is thus of little value and can not be\ntaken as proof that MHC diversity is not driven by pathogen-mediated selection.\nThank you for this remark. Indeed, the lack of effect of estimated viral richness on the\nlink between MHC I and cancer is not proof of the lack of effect of any pathogen-\nmediated selection, especially given the wide diversity of potential intra-cellular\npathogens. Yet, obtaining data on the prevalence, and/or diversity of all these potential\nparasites for all taxa that we studied would be rather difficult. For this reason, and to\navoid misleading our readers, we decided to remove these analyses from the manuscript\nand discuss the potential effect of pathogen-mediated selection on the link between MHC\nI and cancer in the discussion. Please see lines 355-368."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed my criticism. I have only a few further comments.\nLine 108: A paper that may be of interest for you: “Association of HLA diversity with the risk of 25 cancers in the UK Biobank; Wang, Qiao-Ling et al.; eBioMedicine, Volume 92, 104588”\nLine 134: Please add a discussion about the recent studies that did not detect evidence for Peto’s paradox.\n- “Compton, Z. T., Mellon, W., Harris, V. K., Rupp, S., Mallo, D., Kapsetaki, S. E., ... & Boddy, A. M. (2025). Cancer prevalence across vertebrates. Cancer discovery, 15(1), 227-244.”\n- “Butler, G., Baker, J., Amend, S. R., Pienta, K. J., & Venditti, C. (2025). No evidence for Peto’s paradox in terrestrial vertebrates. Proceedings of the National Academy of Sciences, 122(9), e2422861122.”\n- Bulls S. E., et al. , Unraveling the relationship between cancer and life history traits in vertebrates. bioRxiv (2024). 10.1101/2022.07.12.499088\nLine 143: please define “complexity”\nLine 177: What are “appropriate” sequences ?\nLine 372: “naked mole rat (NMR)”. In the manuscript, naked mole rat is sometimes abbreviated, sometimes not. I would suggest not to use the abbreviation.\nFigures : I suggest to add more details in the Y-axes “selection for MHC-I diversity” and “log MHC-I family size” (e.g., dN/dS and number of genes).", "reviewer_2": "The authors have done a good job addressing the comments of all reviewers. In particular, the inclusion of the MHC I copy number variation in response to a comment of one of the reviewers has further strengthened the conclusions. In my opinion, the manuscript is ready for publication.", "reviewer_3": "I have read the revised manuscript and the response by the authors, and appreciate the effort the authors have made to respond to the concerns of the other reviewers and myself. Some of the spurious and potentially misleading analyses have been dropped, and new data with more balanced analyses have been added. I therefore do think that the manuscript and especially the data and analyses have improved significantly.\nHowever, there are still a number of points, both in the response to reviewers, and in the new analyses, that I think need to be clarified, before I can unanimously recommend the publication of this manuscript.\nI am detailing these points below:\n1) Reviewer #1 rightly raised the concern that the dN/dS ratios observed here and interpreted as signatures of 'Selection for MHC diversity' are ranging from 0.25 to 0.7, thus not indicating selection for diversity (dN/dS > 1) at all. Instead, they only indicate relaxed purifying selection, approaching neutrality, which is a very different thing.\nThe authors then respond to this by arguing that selection is known to mostly target the peptide binding region (PBR), which is correct. However, when they calculate dN/dS rations for the PBR sites only, the values are still largely below 1, and only 7 of the species investigated there show values above 1 (and probably not all significantly so). So what the authors use as their main parameter for 'MHC diversity' in reality indicates something else: reduced selection. Furthermore, the actual association between dN/dS and cancer risk disappears when using only the PBR sites (Table S2), which is concerning, as here the association should become stronger, if it was really based on the functional variation of the HLA genes.\n2) In response to my concern about the unbalanced sampling/sequence numbers per species in the NCBI database as a source for the calculation of dN/dS ratios, the authors argue that the number of different MHC sequences in NCBI should in fact reflect the actual MHC polymorphism of the species, providing a mock example that assumes equal sampling between species.\nI spare us an equally likely mock example where unequal sampling between species with equal MHC polymorphism would lead to different MHC sequences numbers. Instead, I want to pick up a point that the authors are making in their response: They argue that the number MHC sequences deposited in NCBI is a valid measure of MHC polymorphism per species ('proven' by their mock example mentioned above). In that case, the number of MHC sequences deposited in NCBI should be a more direct predictor of cancer risk, and thus should be used in the main analyses instead of the more indirect dN/dS ratio. The authors are correct in pointing out that the dN/dS ratio correlates positively with the number of sequences used for its calculation, so the effect should be the same, if not stronger. Interestingly, one paragraph further in their response, the authors then downgrade their argument (and mock example) from above and concede that the number of sequences might indeed be affected by sampling bias, which is why they don't want to include it in their models. However, if the number of sequences is biased by unbalanced sampling, so must be the dN/dS ratio, at it is predicted by the number of sequences used for its calculation.\nEither way, I remain unconvinced that the unsystematic data of MHC sequences in NCBI is a reliable data source for MHC polymorphism of a species. In consequence, I do not have sufficient confidence in the presented results on the sequence numbers or dN/dS ratios to warrant such interpretations as done here.\n3) The authors have added a new analysis on MHC gene copy number across species, which might be able to alleviate some of my concerns, and which I therefore appreciate. However, the MHC region is inherently difficult to assemble, and thus many genome assemblies are lacking in quality and resolution of their MHC regions, correlated with the general effort that went into the given genome assembly. I would feel more comfortable with the estimate of MHC gene family expansion, if the author had also done similar analyses on other known multi-gene families and shown that they do not all correlate with each other (which would indicate a general problem with heterogeneity in the genome assembly quality).\n4) In response to a more conceptual criticism about how cancer could select for many different MHC alleles, the authors argue that heterozygote advantage alone would be sufficient to explain such polymorphism. However, the contrary has been proven in simulation studies (see below for examples), so, I don't see this concern resolved.\nDe Boer, R. J., Borghans, J. A., van Boven, M., Kesmir, C., & Weissing, F. J. (2004). Heterozygote advantage fails to explain the high degree of polymorphism of the MHC. Immunogenetics, 55(11), 725-731. https://doi.org/10.1007/s00251-003-0629-y\nEjsmond, M. J., & Radwan, J. (2015). Red Queen Processes Drive Positive Selection on Major Histocompatibility Complex (MHC) Genes. PLoS Computational Biology, 11(11), e1004627. https://doi.org/10.1371/journal.pcbi.1004627\nRelated to this point, I maintain that throughout the manuscript, the causality between MHC polymorphism and cancer risk/mortality remains unclear (to me): The cancer immune surveillance hypothesis, as I understand it, postulates that an increasing risk for cancer (due to some life-history-related factor that is higher in vertebrates) drove the evolution of the adaptive immune system. The authors here then postulate that this also applies to the evolution of MHC polymorphism, but I would argue that this is a different thing than the evolution of adaptive immunity per se). The cancer immune surveillance hypothesis would be satisfied with having a few MHC gene loci in a genome, and wouldn't require them to be polymorphic, since cancers don't co-evolve.\nSo I do think that the authors need to clarify this conceptual background at set their findings into context. In the end, do not think that their data shows that cancer risk drives the evolution of MHC polymorphism, but it might have contributed to the evolution of adaptive immunity, including the evolution of (neo)antigen-presenting proteins."}, "author_response": "REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have adequately addressed my criticism. I have only a few further comments.\nWe are truly grateful to Reviewer 1 for the meaningful and helpful suggestions to improve\nour work.\nLine 108: A paper that may be of interest for you: “Association of HLA diversity with the risk of 25\ncancers in the UK Biobank; Wang, Qiao-Ling et al.; eBioMedicine, Volume 92, 104588”\nThank you very much for sharing this article. It’s highly relevant, and we now cite this in our\narticle.\nLine 134: Please add a discussion about the recent studies that did not detect evidence for Peto’s\nparadox.\n- “Compton, Z. T., Mellon, W., Harris, V. K., Rupp, S., Mallo, D., Kapsetaki, S. E., ... & Boddy, A.\nM. (2025). Cancer prevalence across vertebrates. Cancer discovery, 15(1), 227-244.”\n- “Butler, G., Baker, J., Amend, S. R., Pienta, K. J., & Venditti, C. (2025). No evidence for Peto’s\nparadox in terrestrial vertebrates. Proceedings of the National Academy of Sciences, 122(9),\ne2422861122.”\n- Bulls S. E., et al., Unravelling the relationship between cancer and life history traits in vertebrates.\nbioRxiv (2024). 10.1101/2022.07.12.499088\nWe truly appreciate the reviewer’s comments, especially for being so up-to-date with current\npublications in comparative oncology. Yet, please note that we purposefully do not mention\nthese articles for the following reasons:\nIn Compton et al. (2025), we did not resolve Peto’s paradox. But indeed, we failed to elaborate\nthat the weak correlation detected between cancer risk and body mass or lifespan does not\nresolve the paradox. Peto predicted a 100-fold increase in cancer risk for every 100-fold\nincrease in body size. What we have shown in Compton et al. and what Butler et al. have\nshown is a 9-fold or 5.7-fold increase (respectively) for a 442,618-fold increase in body size\nbetween mice and elephants. There is a major discrepancy between what Peto predicted and\nthe shallow associations that the statistical models were able to reveal, meaning that the\nparadox is very much valid. Moreover, Due to substantive disagreements regarding\nmethodology and interpretation, we chose not to cite this study. We expand on this\ndisagreement in our recent commentary in PNAS (Compton et al. 2025).\nOn the other hand, Bulls et al. is a preprint, and we are generally careful about citing work\nthat has not undergone peer review and might change considerably when that happens.\nReferences\nCompton, Vincze, Mellon, Tollis, Abegglen, Schiffman, Giraudeau, Boddy, & Maley,\nParadoxical indeed, Proc. Natl. Acad. Sci. U.S.A. 122 (19) e2504512122,\nhttps://doi.org/10.1073/pnas.2504512122 (2025).\nLine 143: please define “complexity”\nThank you for spotting the lack of clarity here. We now specify that we talk about, e.g. higher\ngene copy number or allelic diversity. See lines 143.\nLine 177: What are “appropriate” sequences ?\nThank you for spotting the lack of clarity here. We now rephrased this part to: “We obtained\nMHC I sequence data to estimate selection for MHC I diversity in taxa where cancer\nmortality could also be assessed. Data were available for 28 species for cancer mortality rate\n(CMR) and 27 species for intrinsic cancer mortality (ICM)”.\nLine 372: “naked mole rat (NMR)”. In the manuscript, naked mole rat is sometimes abbreviated,\nsometimes not. I would suggest not to use the abbreviation.\nWe now spell out the species name throughout the manuscript.\nFigures : I suggest to add more details in the Y-axes “selection for MHC-I diversity” and “log\nMHC-I family size” (e.g., dN/dS and number of genes).\nThank you for this suggestion. We now provide more information on the figure axes as\nrecommended.\nReviewer #2 (Remarks to the Author):\nThe authors have done a good job addressing the comments of all reviewers. In particular, the\nauthors have done a good job addressing the comments of all reviewers. In particular, the\ninclusion of the MHC I copy number variation in response to a comment of one of the reviewers has\nfurther strengthened the conclusions. In my opinion, the manuscript is ready for publication.\nWe are truly grateful to Reviewer 2 for the help in improving the work and for the positive\nfeedback on our revision.\nReviewer #3 (Remarks to the Author):\nI have read the revised manuscript and the response by the authors, and appreciate the effort the\nauthors have made to respond to the concerns of the other reviewers and myself. Some of the\nspurious and potentially misleading analyses have been dropped, and new data with more balanced\nanalyses have been added. I therefore do think that the manuscript and especially the data and\nanalyses have improved significantly.\nHowever, there are still a number of points, both in the response to reviewers, and in the new\nanalyses, that I think need to be clarified, before I can unanimously recommend the publication of\nthis manuscript.\nI am detailing these points below:\nWe are grateful to Reviewer 3 for the time and energy dedicated to reviewing and advising\nour analyses. We do – again – agree with the points raised by the Reviewer, and we made\nsignificant changes to respond to the criticism outlined here. Please see the modifications\nexplained in detail below. Overall, we are happy to say that even after changing how we\nmeasure selection for MHC-I diversity, and despite the relatively modest sample size we have\nhere, our results are surprisingly robust.\n1) Reviewer #1 rightly raised the concern that the dN/dS ratios observed here and interpreted as\nsignatures of 'Selection for MHC diversity' are ranging from 0.25 to 0.7, thus not indicating\nselection for diversity (dN/dS > 1) at all. Instead, they only indicate relaxed purifying selection,\napproaching neutrality, which is a very different thing.\nThe authors then respond to this by arguing that selection is known to mostly target the peptide\nbinding region (PBR), which is correct. However, when they calculate dN/dS rations for the PBR\nsites only, the values are still largely below 1, and only 7 of the species investigated there show\nvalues above 1 (and probably not all significantly so). So what the authors use as their main\nparameter for 'MHC diversity' in reality indicates something else: reduced selection. Furthermore,\nthe actual association between dN/dS and cancer risk disappears when using only the PBR sites\n(Table S2), which is concerning, as here the association should become stronger, if it was really\nbased on the functional variation of the HLA genes.\n2) In response to my concern about the unbalanced sampling/sequence numbers per species in the\nNCBI database as a source for the calculation of dN/dS ratios, the authors argue that the number of\ndifferent MHC sequences in NCBI should in fact reflect the actual MHC polymorphism of the\nspecies, providing a mock example that assumes equal sampling between species.\nI spare us an equally likely mock example where unequal sampling between species with equal\nMHC polymorphism would lead to different MHC sequences numbers. Instead, I want to pick up a\npoint that the authors are making in their response: They argue that the number MHC sequences\ndeposited in NCBI is a valid measure of MHC polymorphism per species ('proven' by their mock\nexample mentioned above). In that case, the number of MHC sequences deposited in NCBI should\nbe a more direct predictor of cancer risk, and thus should be used in the main analyses instead of the\nmore indirect dN/dS ratio. The authors are correct in pointing out that the dN/dS ratio correlates\npositively with the number of sequences used for its calculation, so the effect should be the same,\nif not stronger. Interestingly, one paragraph further in their response, the authors then downgrade\ntheir argument (and mock example) from above and concede that the number of sequences might\nindeed be affected by sampling bias, which is why they don't want to include it in their models.\nHowever, if the number of sequences is biased by unbalanced sampling, so must be the dN/dS ratio,\nat it is predicted by the number of sequences used for its calculation.\nEither way, I remain unconvinced that the unsystematic data of MHC sequences in NCBI is a\nreliable data source for MHC polymorphism of a species. In consequence, I do not have sufficient\nconfidence in the presented results on the sequence numbers or dN/dS ratios to warrant such\ninterpretations as done here.\nWe are grateful to the Reviewer for this detailed and constructive comment. We fully agree\nthat this issue requires further investigation, and we have undertaken additional analyses to\naddress your concerns more thoroughly.\n1) Interpretation of dN/dS ratios and the use of PBR-restricted analyses\nWe have carefully explored why dN/dS values restricted to the human-defined peptide-\nbinding region (PBR) do not show stronger associations with cancer risk, as one would expect\nif functional diversity in these regions underpins immunosurveillance efficiency.\nFirst, we investigated the distribution of positively selected sites (PSS) across species in\nrelation to the human-defined MHC-I PBR. We found that only 47.8% of all detected PSS (85\nout of 178 across all species) overlapped with the human PBR. This supports our initial\nsuspicion that using the human PBR definition may fail to adequately capture the full\nspectrum of positive selection in non-model and phylogenetically divergent mammalian taxa.\nWe had previously suggested this as a potential explanation for the weaker associations\nobserved in PBR-restricted analyses.\nImportantly, over half of the signature of positive selection is located outside the human-\ndefined PBR, suggesting that relying on this human-based annotation may overlook\nfunctionally relevant diversity in other species. Additionally, while the human MHC-I PBR\nincludes approximately 40 residues across both exons, we detect a maximum of 21 PSS across\nboth exons in our dataset. This indicates that only a subset of the PBR is under positive\nselection in most species, which likely contributes to the generally low dN/dS ratios observed\neven when restricting analyses to the PBR. At the same time, it is worth stressing that we\ndetected the presence of positively-selected sites in all but two species, providing a clear\nindication that the overwhelming majority of our study taxa experienced some degree of\npositive diversifying selection at the MHC-I.\nIn light of these findings, we have decided to remove the analyses based solely on the human-\ndefined PBR from the manuscript. Instead, we now introduce a more robust and species-\ninformed measure: the dN/dS ratio calculated over the 20 most positively selected sites across\nall taxa. As positive selection is expected to primarily target MHC residues functional in terms\nof antigen binding properties, this approach is likely to capture the most functionally relevant\nvariation, independent of any human-based site classification.\n2) Unbalanced sampling and use of NCBI sequences for dN/dS estimation\nWe also acknowledge and agree that uneven sampling across species in public databases like\nNCBI’s might introduce bias in dN/dS estimates. To address this concern, we have\nimplemented a new analysis using a standardized subsampling procedure. Specifically, we\nrecalculated dN/dS values by randomly subsampling 10 sequences per exon per species (10\niterations per exon), both for all MHC-I sites and for the 20 most positively selected sites (as\nexplained above).\nWhile the sampling effort remains imperfectly balanced (in some species, resampling was not\nfeasible due to a smaller total sample size, i.e., fewer than 10 sequences), this subsampling\napproach substantially reduced variation in sequence number across species and improved\nthe comparability and reliability of our estimates. These updated (resampled) estimates of\nselection for MHC-I diversity were now used to run all models presented in the manuscript.\nValues of selection estimated this way vary from 0.26 to 0.61 across all sites, but from 1.44 to\n3.98 across the 20 most positively selected sites, providing further evidence for the presence of\ndiversifying selection at the MHC-I of our study taxa.\nOverall, replacing our original measures of selection for MHC-I diversity with those derived\nfrom subsampled sequence datasets yielded results that were highly consistent with our\nprevious findings. Specifically, selection for MHC-I diversity continued to show a robust\nnegative correlation with cancer mortality across mammals. These new analyses clearly\nindicate that our key results and conclusions were not driven by variation in sampling effort\nacross species. We hope these revisions address the reviewer’s concerns and further clarify\nthe rationale and robustness of our analytical approach.\n3) The authors have added a new analysis on MHC gene copy number across species, which might\nbe able to alleviate some of my concerns, and which I therefore appreciate. However, the MHC\nregion is inherently difficult to assemble, and thus many genome assemblies are lacking in quality\nand resolution of their MHC regions, correlated with the general effort that went into the given\ngenome assembly. I would feel more comfortable with the estimate of MHC gene family expansion,\nif the author had also done similar analyses on other known multi-gene families and shown that they\ndo not all correlate with each other (which would indicate a general problem with heterogeneity in\nthe genome assembly quality).\nWe fully understand the Reviewer’s concern, but we remain confident in the robustness of\nthis analysis. In fact, one of the authors (Louise Maille) conducted an independent, genome-\nwide investigation aimed at identifying multigene families associated with cancer mortality\nrisk across mammals. This work involved rigorous filtering for genome quality and employed\nphylogenetically controlled single-predictor models. Out of 884 multigene families tested and\nwith stringent correction for multiple comparisons, the MHC gene families emerged among\nthe only nine that consistently and significantly predicted cancer mortality risk across\nmammals. These findings provide independent and genome-wide support for the relevance of\nMHC genes in modulating cancer risk across species. This work is currently being finalised to\nbe submitted for publication, so unfortunately, the data and the results cannot be shared as\nsupplementary material here.\nBesides the reference to the above work, we carefully considered further how it would be best\nto demonstrate that MHC-I gene family expansion is not simply a byproduct of genome\nassembly quality or a general trend across all multi-gene families. Ultimately, we conducted a\ntargeted analysis using all gene families identified in a genome-wide Orthofinder run across\n106 mammalian species (where we had data on cancer mortality risk). To avoid confounding\neffects from immune-related gene clusters, we focused on gene families that did not contain\ngenes previously classified as immune-related in the ImmPort database (Bhattacharya et al.\n2018) and (2) exhibited sufficient cross-species variation, defined as at least 30% of species\ndiffering from the cross-species mean in copy number, and at least 30% of the species had at\nleast one gene in the gene family. This filtering yielded 686 gene families for comparison. We\nthen calculated Spearman rank correlations between MHC-I gene family size and the copy\nnumber of each of these 686 gene families, which allowed us to assess the strength and\ndirection of any associations while accounting for non-Gaussian distributions. The resulting\nhistogram of correlation coefficients (see below) shows a distribution centred near zero (mean\n≈ 0.08), with both positive and negative values. While a slight bias toward positive correlations\nexists (likely reflecting general genome-wide expansion trends, such as genome-wide\nduplication events, co-expansion of functionally linked families, etc), the presence of\nnumerous negative correlations clearly indicates that not all gene families co-expand with\nMHC-I. This supports our conclusion that MHC-I gene family expansion is not merely a\nreflection of overall genome assembly characteristics, but rather exhibits specific and\nbiologically meaningful variation.\nReferences\nBhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, Hu Z, Zalocusky KA,\nShankar RD, Shen-Orr SS, Thomson E, Wiser J, Butte AJ. ImmPort, toward repurposing of\nopen access immunological assay data for translational and clinical research. Sci Data. 2018\nFeb 27;5:180015. doi: 10.1038/sdata.2018.15.\n4) In response to a more conceptual criticism about how cancer could select for many different\nMHC alleles, the authors argue that heterozygote advantage alone would be sufficient to explain\nsuch polymorphism. However, the contrary has been proven in simulation studies (see below for\nexamples), so, I don't see this concern resolved.\nDe Boer, R. J., Borghans, J. A., van Boven, M., Kesmir, C., & Weissing, F. J. (2004). Heterozygote\nadvantage fails to explain the high degree of polymorphism of the MHC. Immunogenetics, 55(11),\n725-731. https://doi.org/10.1007/\nEjsmond, M. J., & Radwan, J. (2015). Red Queen Processes Drive Positive Selection on Major\nHistocompatibility Complex (MHC) Genes. PLoS Computational Biology, 11(11), e1004627.\nhttps://doi.org/10.1371/\nRelated to this point, I maintain that throughout the manuscript, the causality between MHC\npolymorphism and cancer risk/mortality remains unclear (to me): The cancer immune surveillance\nhypothesis, as I understand it, postulates that an increasing risk for cancer (due to some life-history-\nrelated factor that is higher in vertebrates) drove the evolution of the adaptive immune system. The\nauthors here then postulate that this also applies to the evolution of MHC polymorphism, but I\nwould argue that this is a different thing than the evolution of adaptive immunity per se). The\ncancer immune surveillance hypothesis would be satisfied with having a few MHC gene loci in a\ngenome, and wouldn't require them to be polymorphic, since cancers don't co-evolve.\nSo I do think that the authors need to clarify this conceptual background at set their findings into\ncontext. In the end, do not think that their data shows that cancer risk drives the evolution of MHC\npolymorphism, but it might have contributed to the evolution of adaptive immunity, including the\nevolution of (neo)antigen-presenting proteins.\nWe appreciate the reviewer's thoughtful critique and the opportunity to clarify the conceptual\nframework of our study. We have now carefully read the cited articles and we improved the\ntext to reflect that:\n1. Our manuscript does not propose heterozygote advantage in cancer defense as the sole\nevolutionary mechanism driving MHC diversity. We now clearly state (lines 370–397)\nthat MHC diversity is shaped by multiple and interacting evolutionary processes\n(Radwan et al. 2020), including selection pressures from parasites, pathogens, mate\nchoice and potentially cancer. We also explicitly acknowledge that our study does not\naim to resolve the long-standing debate on the relative contributions of different\nmechanisms maintaining MHC polymorphism, but rather highlights cancer as an\nunder-explored and potentially contributing factor.\n2. We also tried to clarify that our goal in this article was not to revisit the evolutionary\norigins of MHC polymorphism, but to examine whether the species-specific MHC I\ncomplexity predicts cancer risk across taxa. While related, this question is distinct\nfrom asking whether MHC polymorphism predicts cancer predisposition.\n3. While we acknowledge that models of heterozygote advantage often fail to explain the\nextreme MHC polymorphism observed in natural populations, this might not be\nindependent of the inherent simplifications in such models. Specifically, many models\ndo not fully incorporate the complexities of host immune ecology, including T-cell\nreceptor repertoire selection, the biochemical constraints on evolving broadly binding\n(“generalist”) MHC molecules, and the potential immunopathological costs associated\nwith such generalism. For example, generalist MHC alleles exhibit reduced expression\nlevels and lower peptide-binding affinity or specificity, potentially as a mechanism to\nmitigate self-reactivity (Chappell et al., 2015). Moreover, although recent\nmathematical and empirical studies have improved our understanding, heterozygote\nadvantage remains one of the hypotheses for the maintenance of MHC diversity.\nIndeed, it’s true that most cancers do not co-evolve with the hosts as pathogens do, but\n(1) transmissible cancers do very much act as pathogens, they can be transferred\nbetween hosts and (2) many cancers (approximately 20% in humans) is of pathogenic\norigin. At least these specific cases might reinforce the selection mechanisms already\nknown in case of pathogens. We now express this in the manuscript.\n4. Second, indeed, the cancer immunosurveillance hypothesis proposes that the adaptive\nimmune system evolved in response to the need to detect and eliminate transformed\ncells. And while adaptive immunity has multiple components, MHC is an essential part\nof it. On the contrary, MHC polymorphism is indeed not necessary for adaptive\nimmunity to exist, yet MHC polymorphism can substantially enhance the efficacy of\nadaptive immune mechanisms, particularly in a context of somatic evolution, where\nthe healthy and transformed cells are very similar (unlike in the case of pathogens).\nThe beneficial role of high MHC I diversity in cancer is also well supported by multiple\nrecent publications indicating lower incidence of diverse tumor types in individuals\nwith more diverse MHC (Wang et al. 2023), as well as a higher efficiency of\nimmunotherapies in patients with higher allele diversity (Chowell et al. 2019). We now\ncite these supporting works to solidify our argument on the importance of MHC-I in\ncancer defense.\nReferences\nChappell, PA, El Kahina Meziane, Michael Harrison, Łukasz Magiera, Clemens Hermann,\nLaura Mears, Antoni G Wrobel, Charlotte Durant, Lise Lotte Nielsen, Søren Buus, Nicola\nTernette, William Mwangi, Colin Butter, Venugopal Nair, Trudy Ahyee, Richard Duggleby,\nAlejandro Madrigal, Pietro Roversi, Susan M Lea, Jim Kaufman (2015) Expression levels of\nMHC class I molecules are inversely correlated with promiscuity of peptide binding eLife\n4:e05345\nChowell, D. et al. Evolutionary divergence of HLA class I genotype impacts efficacy of cancer\nimmunotherapy. Nat Med 25, 1715–1720 (2019).\nRadwan, J., Babik, W., Kaufman, J., Lenz, T. L., & Winternitz, J. (2020). Advances in the\nevolutionary understanding of MHC polymorphism. Trends in Genetics, 36(4), 298-311.\nWang, Q.-L. et al. Association of HLA diversity with the risk of 25 cancers in the UK\nBiobank. eBioMedicine 92, (2023)."}, {"round": 3, "reviewer_comments": {"reviewer_3": "I appreciate the substantial additional work the authors have put into responding to my continuous concerns. I acknowledge that the robustness of the main associations to these additional analyses is quite convincing. I also appreciate the clarifications towards the conceptual background and evolution of MHC diversity in general. My main concerns have been addressed, and I do not want to hold this manuscript up any further. I do think that this will be an interesting and thought- provoking contribution to the literature on cancer and MHC evolution."}, "author_response": ""}]

Decision letter