Heterozygote advantage cannot explain MHC diversity, but MHC diversity can explain heterozygote advantage

  1. Joshua L Cherry

  • Curated by eLife

    eLife logo

    eLife Assessment

    This valuable study re-evaluates a published simulation model on the role of heterozygote advantage in shaping MHC diversity. By modifying key modeling assumptions, the author argues that the original conclusions depend on a narrow and potentially unrealistic parameter range. While the work is in principle solid, the robustness of this claim is viewed differently by the reviewers. The manuscript further proposes an alternative modeling framework in which expansion of the MHC gene family allows homozygotes to outperform heterozygotes, thereby challenging the idea that heterozygote advantage alone can account for high allelic diversity at MHC loci. The topic is highly relevant for eco-immunology and evolutionary genetics, although a clearer delineation of the model's scope would help readers assess its broader implications.

Reviewed by

Read the full article See related articles

Discuss this preprint

Start a discussion What are Sciety discussions?

Listed in

Log in to save this article

Abstract

Several theoretical studies have concluded that heterozygote advantage makes at most a minor contribution to MHC diversity. Siljestam and Rueffler (2024) recently presented models in which heterozygote advantage alone can lead to realistically high diversity. Here I argue that heterozygote advantage cannot by itself explain MHC diversity, and that its contribution to diversity is unlikely to be large in most species. I first show that the high diversity reported by Siljestam and Rueffler is so sensitive to parameter values that the underlying phenomenon cannot explain the widespread diversity of MHC genes. I then consider a fundamental problem with explaining MHC diversity by heterozygote advantage alone: selective forces that favored heterozygotes would lead to the evolution of haplotypes having much higher fitness when homozygous, diminishing or eliminating heterozygote advantage. Diversity maintained by another force, however, might bring about adaptation to the more common heterozygous state at the expense of homozygous fitness. Thus, substantial heterozygote advantage may arise as a consequence of MHC diversity.

Article activity feed

  1. eLife

    eLife Assessment

    This valuable study re-evaluates a published simulation model on the role of heterozygote advantage in shaping MHC diversity. By modifying key modeling assumptions, the author argues that the original conclusions depend on a narrow and potentially unrealistic parameter range. While the work is in principle solid, the robustness of this claim is viewed differently by the reviewers. The manuscript further proposes an alternative modeling framework in which expansion of the MHC gene family allows homozygotes to outperform heterozygotes, thereby challenging the idea that heterozygote advantage alone can account for high allelic diversity at MHC loci. The topic is highly relevant for eco-immunology and evolutionary genetics, although a clearer delineation of the model's scope would help readers assess its broader implications.

  2. eLife

    Reviewer #1 (Public review):

    The manuscript "Heterozygote advantage cannot explain MHC diversity, but MHC diversity can explain heterozygote advantage" explores two topics. First, it is claimed that the recently published by Mattias Siljestam and Claus Rueffler conclusion (in the following referred to as [SR] for brevity) that heterozygote advantage explains MHC diversity does not withstand an even very slight change in ecological parameters. Second, a modified model that allows an expansion of MHC gene family shows that homozygotes outperform heterozygotes. This is an important topic and could be of potential interest to the readership of eLife if the conclusions are valid and non-trivial.

    The resubmitted manuscript addresses several questions from my previous review. In particular, there is a more detailed description of how the code of Siljestam and Rueffler ([SR]) was used for the simulations and the calculation of the factor 2.7 x 10^43 that is the key to the alleged breakdown of the numerical reasoning presented by in [SR].

    Yet I think that important aspects of my critique of the first statement of the manuscript about the flaws of [SR] model remain unanswered. I guess the discussion becomes rather general about the universality and robustness of various types of models to parameter changes. My point is that none of the models is totally universal. The model in [SR] is not phenomenological as none of the parameters or functional forms were derived empirically. Instead, it is a proof of principle demonstration that inevitably grossly simplifies the actual immune response. The choice of constants and functions used in Eqs. (1-5) is dictated by the mathematical convenience and works in a limited range of parameter values. It is shown in [SR] that for 3 pathogens and reasonable "virulence " \nu, the alleles branch. These conclusions are supported by the analytically derived Adaptive Dynamics branching criteria (7), which, contrary to the statement is the cover letter (" It is clear from Fig. 4 of Siljestam and Rueffler that the branching condition is far from sufficient for high MHC diversity.") is perfectly confirmed by the simulation data shown in Fig. 4.

    The mathematical simplicity of the [SR] model generates various artifacts, such as the mentioned by the Author reduction of the "condition" by an enormous factor 2.7 x 10^43 and the resulting decrease in the "survival" induced by the addition of a new pathogen. This occurs at the very large value of \nu=20, whose effect is enormous due to the Gaussian form of (1), which, once again, was chosen for the mathematical convenience. In reality, a new pathogen cannot reduce the "survival" by such a factor as it would wipe out any resident population. So to compensate for such an artifact, the additional factor c_max was introduced to buffer such an excess. There is no reason to fix c_max once for an arbitrary number of pathogens, because varying c_max basically reflects the observation that a well-adapted individual must have a reasonable survival probability. At the same time, there are many ways in which the numerical simulation may break down when the survival rates become of the order of 10^(-43) instead of one, so it comes to no surprise that the diversification, predicted by the adaptive dynamics, does not readily occur in the scenario with an addition or removal of the 8th pathogen with a very high virulence \nu=20.

    I have doubts that the reported breakdown of the [SR] model with fixed c_max remains observable with less extreme values of m and \nu (say, for \nu=7 and m=3 plus or minus 1 used in Fig. 3 in the manuscript).

    So I still find the claim that " the phenomenon that leads to high diversity in the simulations of Siljestam and Rueffler depends on finely tuned parameter values" is not well substantiated.

  3. eLife

    Reviewer #2 (Public review):

    Summary:

    This study addresses the population genetic underpinnings of the extraordinary diversity of genes in the MHC, which is widespread among jawed vertebrates. This topic has been widely discussed and studied, and several hypotheses have been suggested to explain this diversity. One of them is based on the idea that heterozygote genotypes have an advantage over homozygotes. While this hypothesis lost early on support, a reason study claimed that there is good support for this idea. The current study highlights an important aspect that allows us to see results presented in the earlier published paper in a different light, changing strongly the conclusions of the earlier study, i.e., there is no support for a heterozygote advantage. This is a very important contribution to the field. Furthermore, this new study presents an alternative hypothesis to explain the maintenance of MHC diversity, which is based on the idea that gene duplications can create diversity without heterozygosity being important. This is an interesting idea, but not entirely new.

    Strength:

    (1) A careful re-evaluation of a published model, questioning a major assumption made by a previous study.

    (2) A convincing reanalysis of a model that, in the light of the re-analysis-loses all support.

    (3) A convincing suggestion for an alternative hypothesis.

    Weakness:

    (1) The title of the study is catchy, but it is explained only in the very end of the paper.

  4. eLife

    Author response:

    The following is the authors’ response to the current reviews.

    Reviewer #1:

    Yet I think that important aspects of my critique of the first statement of the manuscript about the flaws of [SR] model remain unanswered.

    I believe that I have fully addressed the points in the earlier review. The reviewer had doubted that my results were correct, attributing them to “a poor setup of the model” on my part. The reviewer stated that if I were correct about the factor of >1043 change in cmax, this would “naturally break down all the estimates and conclusions made in Siljestam and Rueffler” (S&R).

    It appears that the reviewer is now convinced that my results represent a faithful analysis of the models on which S&R based their claims. The reviewer now contends that these results, including the factor of >1043, present no difficulties for the claims of S&R after all. In fact, this enormous factor of >1043 is now claimed to support the conclusions of S&R by invalidating my conclusions. I respond to these new and very different arguments in what follows.

    As I stated in the first round of review, the issue is not the enormity of this factor per se, but the fact that the compensatory adjustment of cmax conceals the true effects of changes in other parameters. These effects are large; small changes to the parameter values mostly eliminate the diversity that the model is claimed to explain.

    The model in [SR] is not phenomenological as none of the parameters or functional forms were derived empirically. Instead, it is a proof of principle demonstration that inevitably grossly simplifies the actual immune response.

    The hidden sensitivity of the results of S&R to paramater values is sufficient to invalidate them as a proof of principle. The manuscript goes further and explains how the problem "is not specific to the details of the models of Siljestam and Rueffler, but is inherent in the phenomenon invoked to allow high diversity" because "any change that affects condition by as much as the difference between MHC heterozygotes and homozygotes will eliminate high equilibrium diversity". This general principle addresses all of the reviewer's points.

    In reality, a new pathogen cannot reduce the "survival" by such a factor as it would wipe out any resident population. So to compensate for such an artifact, the additional factor cmax was introduced to buffer such an excess. There is no reason to fix cmax once for an arbitrary number of pathogens, because varying cmax basically reflects the observation that a well-adapted individual must have a reasonable survival probability.

    This is not a legitimate reason for making compensatory, diversity-promoting adjustments to cmax when evaluating sensitivity to other parameters. If the number of pathogens or their virulence changes, cmax obviously does not automatically change along with it. If the population or species consequently goes extinct, then it goes extinct. If it persists, it does so with the same value of cmax.

    The possibility of extinction arguably puts a minimum value on cmax, but it does not restrict it to a range of values that conveniently leads to high MHC diversity. In the examples that I analyzed, slightly decreasing the number of pathogens or their virulence, which increases survivability, eliminates diversity. This phenomenon obviously cannot be dismissed on the grounds that survivability would be too low for the species to exist.

    S&R in effect assume that the condition of the most fit homozygote remains fixed, regardless of the number of pathogens, their virulence, and myriad other differences between species. It is this assumption that is without justification.

    At the same time, there are many ways in which the numerical simulation may break down when the survival rates become of the order of 10^(-43) instead of one

    I am not sure what is meant by “the numerical simulation may break down”. Numerical error is not a tenable explanation of the lack of diversity observed in that simulation. The outcome is exactly what is expected from purely theoretical considerations: conditions of all genotypes fall on the steep part of the curve, making the mechanism proposed by S&R largely inoperative, so a pair of alleles forming a fit heterozygote comes to predominate. The numerical simulation is actually superfluous.

    Low survival rates are completely irrelevant to the effect of decreasing the number of pathogens or their virulence, which does not lower survival rates, but does eliminate diversity.

    so it comes to no surprise that the diversification, predicted by the adaptive dynamics, does not readily occur in the scenario with an addition or removal of the 8th pathogen with a very high virulence \nu=20.

    Whether or not it surprising, the lack of diversity is a problem for the claims of S&R, as there is no reason to expect the number of pathogens to have just the right value to produce high diversity. Furthermore, for many combinations of values of the other parameters (e.g., my v=19.5 and 20.5 examples), no number of pathogens leads to high diversity.

    Again, the general principle mentioned above makes the details that the reviewer refers to irrelevant. Nonetheless, some additional remarks are in order:

    (1) This comment ignores the fact that removal of a pathogen, or a slight decrease in “virulence”, eliminates diversity without lowering survival rates.

    (2) Small increases or decreases in v (virulence) eliminate diversity without having such large effects on condition.

    (3) In the example emphasized by the reviewer, mean survival rates are nowhere near as low as 10-43. Only homozygotes have such low fitness.

    (4) The adaptive dynamics predict the low diversity seen in the simulations, contrary to what the reviewer seems to suggest. Elimination of diversity is not an artifact of the simulation.

    (5) v=20 was chosen because it is most favorable to the model of S&R in that it yields the highest diversity. Indeed, S&R only observed realistically high diversity with the narrow gaussians that the reviewer objects to. With lower values of v, diversity is much lower, but even this meager diversity is eliminated by small changes in parameter values (see below). If narrow gaussians and large effects of pathogens somehow invalidate results, then they invalidate the high-diversity results of S&R.

    I have doubts that the reported breakdown of the [SR] model with fixed cmax remains observable with less extreme values of m and \nu (say, for \nu=7 and m=3 plus or minus 1 used in Fig. 3 in the manuscript).

    These doubts are unwarrented. With the suggested parameter values, for example, increasing or decreasing m by 1 reduces the effective number of alleles to around 1 or 2. This can easily be checked using the simulation code of S&R, as detailed in my initial response and now in a Supplementary Text. Even without this result, the general principle mentioned above tells us that considering other regions of parameter space cannot rescue the conclusions of S&R.

    So I still find the claim that " the phenomenon that leads to high diversity in the simulations of Siljestam and Rueffler depends on finely tuned parameter values" is not well substantiated.

    What is unsubstantiated is the claim of S&R that “For a large part of the parameter space, more than 100 and up to over 200 alleles can emerge and coexist”. As my manuscript illustrates, this is an illusion created by the adjustment of one parameter to compensate for changes in others.

    The reviewer even acknowledges that “the choice of constants and functions...works in a limited range of parameter values”. Furthermore, the manuscript explains why this problem is inherent to the general phenomenon, not specific to the details of the model or parameter values.

    The following is the authors’ response to the original reviews.

    Public Reviews:

    Reviewer #1 (Public review):

    It appears obvious that with no or a little fitness penalty, it becomes beneficial to have MHC-coding genes specific to each pathogen. A more thorough study that takes into account a realistic (most probably non-linear in gene number) fitness penalty, various numbers of pathogens that could grossly exceed the self-consistent fitness limit on the number of MHC genes, etc, could be more informative.

    The reviewer seems to be referring to the cost of excessively high presentation breadth. Such a cost is irrelevant to the inferior fitness of a polymorphic population with heterozygote advantage compared to a monomorphic population with merely doubled gene copy number. It is relevant to the possibility of a fitness valley separating these two states, but this issue is addressed explicitly in the manuscript.

    An addition or removal of one of the pathogens is reported to affect "the maximum condition", a key ecological characteristic of the model, by an enormous factor 10^43, naturally breaking down all the estimates and conclusions made in [RS]. This observation is not substantiated by any formulas, recipes for how to compute this number numerically, or other details, and is presented just as a self-standing number in the text.

    It is encouraging that the reviewer agrees that this observation, if correct, would cast doubt on the conclusions of Siljestam and Rueffler. I would add that it is not the enormity of this factor per se that invalidates those conclusions, but the fact that the automatic compensatory adjustment of cmax conceals the true effects of removing a pathogen, which are quite large.

    I am not sure why the reviewer doubts that this observation is correct. The factor of 2.7∙1043 was determined in a straightforward manner in the course of simulating the symmetric Gaussian model of Siljestam and Rueffler with the specified parameter values. A simple way to determine this number is to have the simulation code print the value to which cmax is set, or would be set, by the procedure of Siljestam and Rueffler for different parameter values. I have in this way confirmed this factor using the simulation code written and used by Siljestam and Rueffler. A procedure for doing so is described in the new Supplementary Text S1. In addition, I now give a theoretical derivation of this factor in Supplementary Text S2.

    This begs the conclusion that the branching remains robust to changes in cmax that span 4 decades as well.

    That shows at most that the results are not extremely sensitive to cmax or K. They are, nonetheless, exquisitely sensitive to m and v. This difference in sensitivities is the reason that a relatively small change to m leads to such a large compensatory change in cmax. It is evident from Fig. 4 of Siljestam and Rueffler that the level of diversity is not robust to these very large changes in cmax, which include, as noted above, a change of over 43 orders of magnitude.

    As I wrote above, there is no explanation behind this number, so I can only guess that such a number is created by the removal or addition of a pathogen that is very far away from the other pathogens. Very far in this context means being separated in the x-space by a much greater distance than 1/\nu, the width of the pathogens' gaussians. Once again, I am not totally sure if this was the case, but if it were, some basic notions of how models are set up were broken. It appears very strange that nothing is said in the manuscript about the spatial distribution of the pathogens, which is crucial to their effects on the condition c.

    I did not explicitly describe the distribution of pathogens in antigenic space because it is exactly the same as in Siljestam and Rueffler, Fig. 4: the vertices of a regular simplex, centered at the origin, with unity edge length.

    The number in question (2.7∙1043) pertains to the Gaussian model with v=20. As specified by Siljestam and Rueffler, each pathogen lies at a distance of 1 from every other pathogen, so the distance of any pathogen from the others is indeed much greater than 1/v. This condition holds, however, for most of the parameter space explored by Siljestam and Rueffler (their Fig. 4), and for all of the parameter space that seemingly supports their conclusions. Thus, if this condition indicates that “basic notions of how models are set up were broken”, they must have been broken by Siljestam and Rueffler.

    ...the branching condition appears to be pretty robust with respect to reasonable changes in parameters.

    It is clear from Fig. 4 of Siljestam and Rueffler that the branching condition is far from sufficient for high MHC diversity.

    Overall, I strongly suspect that an unfortunately poor setup of the model reported in the manuscript has led to the conclusions that dispute the much better-substantiated claims made in [SD].

    The reviewer seems to be suggesting that my simulations are somehow flawed and my conclusions unreliable. I have addressed the reasons for this suggestion above. Furthermore, I have confirmed the main conclusion—the extreme sensitivity of the results of Siljestam and Rueffler to parameter values--using the code that they used for their simulations, indicating that my conclusions are not consequences of my having done a “poor setup of the model”. I now describe, in Supplementary Text S1, how anybody can verify my conclusions in this way.

    Reviewer #2 (Public review):

    (1) The statement that the model outcome of Siljestam and Rueffler is very sensitive to parameter values is, in this form, not correct. The sensitivity is only visible once a strong assumption by Siljestam and Rueffler is removed. This assumption is questionable, and it is well explained in the manuscript by J. Cherry why it should not be used. This may be seen as a subtle difference, but I think it is important to pin done the exact nature of the problem (see, for example, the abstract, where this is presented in a misleading way).

    I appreciate the distinction, and the importance of clearly specifying the nature of the problem. However, as I understand it, Siljestam and Rueffler do not invoke the implausible assumption that changes to the number of pathogens or their virulence will be accompanied by compensatory changes to cmax. Rather, they describe the adjustment of cmax (Appendix 7) as a “helpful” standardization that applies “without loss of generality”. Indeed, my low-diversity results could be obtained, despite such adjustment, by combining the small change to m or v with a very large change to K (e.g., a factor of 2.7∙1043). In this sense there is no loss of generality, but the automatic adjustment of cmax obscures the extreme sensitivity of the results to m and v.

    (2) The title of the study is very catchy, but it needs to be explained better in the text.

    I have expanded the end of the Discussion in the hope of clarifying the point expressed by the title.

    Recommendations for the authors:

    Reviewer #1 (Recommendations for the authors):

    I would like to suggest to the author that they provide essential details about their simulations that would justify their claims, and to communicate with Mattias Siljestam and Claus Rueffler whether claims of the lack of robustness could be confirmed.

    The models simulated were modified versions of those of Siljestam and Rueffler. Thus, only the modifications were described in my manuscript. I have added a more detailed description of how cmax was set in the simulations concerned with sensitivity to parameter values. In addition, the new Supplementary Text S1, which describes confirmation of the lack of robustness using the code of Siljestam and Rueffler, should remove any doubt about this conclusion.

    Reviewer #2 (Recommendations for the authors):

    I have no further recommendations. The manuscript is well written and clear.

    Thank you.

    Reviewer #3 (Recommendations for the authors):

    (1) Since this is a full report and not just a letter to the editor, it would benefit from a bit more introduction of what the MHC actually is and what the current understanding of its evolution is. Currently, it assumes a lot of knowledge about these genes that might not be available to every reader of eLife.

    I have added some more information to the opening paragraph. I would also note that this report was submitted as a “Research Advance”, which may only need “minimal introductory material”.

    (2) Some more recent literature on MHC evolution should be added, e.g., the review by Radwan et al. 2020 TiG, a concrete case of MHC heterozygote advantage by Arora et al. 2020 MolBiolEvol, and a simulation of MHC CNV evolution by Bentkowski et al. 2019 PLOSCompBiol.

    I have cited some additional literature.

    (3) Since much of the criticism hinges on the cmax parameter, its biological meaning or role (or the lack thereof) could be discussed more.

    I am not sure what I can add to what is in the first paragraph of the Discussion.

    (4) I find it difficult to grasp how the v parameter, which is intended to define pathogen virulence, if I understand it correctly, can be used to amend the breadth of peptide presentation. Maybe this could be illustrated better.

    I have attempted to make this clearer. The parameter v actually controls the breadth of peptide detection conferred by an allele, which, if not identical to the breath of presentation, is certainly affected by it. The basis of the “virulence” interpretation seems to be that narrower detection breadth can, according to the model, only decrease peptide detection probability, which increases the damage done by pathogens.

    (5) Please check sentences in lines 279ff on peptide detection and cost of . There seem to be words missing.

    There was an extraneous word, which I have removed. Thank you for pointing this out.

  5. Version published to 10.7554/elife.107256.2 on eLife
  6. Version published to 10.7554/elife.107256 on eLife
  7. Version published to 10.7554/elife.107256.1 on eLife
  8. eLife

    eLife Assessment

    The is a valuable evaluation of a previously published simulation model on the role of heterozygote advantage in shaping MHC diversity, showing that the conclusions from this model hold only within a narrow parameter range that might be unrealistic. The author presents an alternative model, in which MHC homozygotes with duplicated MHC genes outperform heterozygotes with single genes, thereby challenging the explanation that heterozygote advantage will lead to high allelic variation at a given MHC gene. The topic is highly relevant for eco-immunology and evolutionary genetics, but several major aspects of the author's claim need to be clarified to make the model interpretable. While the work has the potential to improve our understanding of the question of how the extraordinary diversity at the MHC locus evolves, without this addition, the conclusions remain incomplete.

  9. eLife

    Reviewer #1 (Public review):

    The manuscript "Heterozygote advantage cannot explain MHC diversity, but MHC diversity can explain heterozygote advantage" explores two topics. First, it is claimed that the recently published conclusion by Mattias Siljestam and Claus Rueffler (in the following referred to as [SR] for brevity) that heterozygote advantage explains MHC diversity does not withstand even a very slight change in ecological parameters. Second, a modified model that allows an expansion of the MHC gene family shows that homozygotes outperform heterozygotes. This is an important topic and could be of potential interest to readers if the conclusions are valid and non-trivial.

    Let me first comment on the second part of the manuscript that describes the fitness advantage of the 'gene family expansion'. I think this, by itself, is a totally predictable result. It appears obvious that with no or a little fitness penalty, it becomes beneficial to have MHC-coding genes specific to each pathogen. A more thorough study that takes into account a realistic (most probably non-linear in gene number) fitness penalty, various numbers of pathogens that could grossly exceed the self-consistent fitness limit on the number of MHC genes, etc, could be more informative. Yet, as I understood the narrative of the manuscript, the expansion of the gene family serves as a mere counter-example to the disputed finding of [SR], rather than a systematic study of the eco-evolutionary consequences of this process.

    Now to the first part of the manuscript, which claims that the point made in [RS] is not robust and breaks down under a small change in the parameters. An addition or removal of one of the pathogens is reported to affect "the maximum condition", a key ecological characteristic of the model, by an enormous factor 10^43, naturally breaking down all the estimates and conclusions made in [RS]. This observation is not substantiated by any formulas, recipes for how to compute this number numerically, or other details, and is presented just as a self-standing number in the text. The only piece of information given in the manuscript is that, unlike in [SR], the adjustable parameter c_{max} is kept constant when the number of pathogens is changed.

    In my opinion, the information provided in the manuscript does not allow one to conclude anything about the relevance and the validity of its main claim. At the same time, the simulations done in [SR] are described with a fair amount of detail. Which allows me to assume that the conclusions made in [SR] are fairly robust and, in particular, have been demonstrated not to be too sensitive to changes in the main "suspect', c_{max}. Let me briefly justify my point.

    First, it follows from Eqs (4,5) in the main text and (A12-A13) in the Appendix that c_{max} and K do not independently affect the dynamics of the model, but it's rather their ratio K/c_max that matters. It can be seen by dividing the numerator and denominator of (5) by c_max. Figure 3 shows the persistent branching for 4 values of K that cover 4 decades. As it appears from the schemes in the top row of Figure 3, those simulations are done for the same positions and widths/virulences of pathogens. So the position of x* should be the same in all 4 cases, presumably being at the center of pathogens, (x*,x*) = (0,0). According to the definition of x* given in the Appendix after Eqs (A12-A13), this means that c_max remains the same in all 4 cases. So one can interpret the 4 scenarios shown in Figure 3 as corresponding not to various K, but to various c_max that varied inversely to K. That is, the results would have been identical to those shown in Figure 3 if K were kept constant and c_max were multiplied by 0.1, 1, 10, and 100, or scaled as 1/K. This begs the conclusion that the branching remains robust to changes in c_max that span 4 decades as well.

    Naturally, most, if not all, the dynamics will break down if one of the ecological characteristics changes by a factor of 10^43, as it is reported in the submitted manuscript. As I wrote above, there is no explanation behind this number, so I can only guess that such a number is created by the removal or addition of a pathogen that is very far away from the other pathogens. Very far in this context means being separated in the x-space by a much greater distance than 1/\nu, the width of the pathogens' gaussians. Once again, I am not totally sure if this was the case, but if it were, some basic notions of how models are set up were broken. It appears very strange that nothing is said in the manuscript about the spatial distribution of the pathogens, which is crucial to their effects on the condition c. In [SP], it is clearly shown where the pathogens are.

    Another argument that makes me suspicious in the utility of the conclusions made in the manuscript and plays for the validity of [SP] is the adaptive dynamics derivation of the branching conditions. It is confirmed by numerics with sufficient accuracy, and as it stands in its simple form of the inequality between two widths, the branching condition appears to be pretty robust with respect to reasonable changes in parameters.

    Overall, I strongly suspect that an unfortunately poor setup of the model reported in the manuscript has led to the conclusions that dispute the much better-substantiated claims made in [SD].

  10. eLife

    Reviewer #2 (Public review):

    Summary:

    This study addresses the population genetic underpinnings of the extraordinary diversity of genes in the MHC, which is widespread among jawed vertebrates. This topic has been widely discussed and studied, and several hypotheses have been suggested to explain this diversity. One of them is based on the idea that heterozygote genotypes have an advantage over homozygotes. While this hypothesis lost early on support, a reason study claimed that there is good support for this idea. The current study highlights an important aspect that allows us to see results presented in the earlier published paper in a different light, changing strongly the conclusions of the earlier study, i.e., there is no support for a heterozygote advantage. This is a very important contribution to the field. Furthermore, this new study presents an alternative hypothesis to explain the maintenance of MHC diversity, which is based on the idea that gene duplications can create diversity without heterozygosity being important. This is an interesting idea, but not entirely new.

    Strengths:

    (1) A careful re-evaluation of a published model, questioning a major assumption made by a previous study.

    (2) A convincing reanalysis of a model that, in the light of the re-analysis-loses all support.

    (3) A convincing suggestion for an alternative hypothesis.

    Weaknesses:

    (1) The statement that the model outcome of Siljestam and Rueffler is very sensitive to parameter values is, in this form, not correct. The sensitivity is only visible once a strong assumption by Siljestam and Rueffler is removed. This assumption is questionable, and it is well explained in the manuscript by J. Cherry why it should not be used. This may be seen as a subtle difference, but I think it is important to pin done the exact nature of the problem (see, for example, the abstract, where this is presented in a misleading way).

    (2) The title of the study is very catchy, but it needs to be explained better in the text.

  11. eLife

    Reviewer #3 (Public review):

    This manuscript describes a careful and thorough evaluation of an evolutionary simulation model published previously. The model and this report address the question, whether heterozygote advantage (HA) by itself as a selection mechanism can explain a substantial level of allelic diversity as it is often seen in MHC immune genes. Despite decades of research on the topic of pathogen-mediated selection for MHC diversity, it remains an open question by which specific selection mechanisms this exceptional allelic diversity is maintained.

    The previously published paper posits, in contrast to various previous studies, that HA is, in fact, able to maintain a level of allelic diversity as seen in many populations, just by itself, given certain conditions. The current manuscript now challenges this conclusion by highlighting that the previous model results only hold under very narrow parameter ranges.

    Besides criticizing some of the conceptual points of the previous paper, the author carefully rebuilt the previously published model and replicated their results, before then evaluating the robustness of the model results to reasonable variation in different parameters. From this evaluation, it becomes clear that the previously reported results hinge strongly on a certain scaling or weighing factor that is adjusted for every parameter setting and essentially counteracts the changes induced by changing the parameters. The critical impact of this one parameter is not clearly stated in the previous paper, but raises serious doubts about the generalizability of the model to explain MHC allelic variation across diverse vertebrate species.

    Given the fact that the MHC genes are among the most widely studied genes in vertebrates, and that understanding their evolution will shed light on their association with various complex diseases, the insights from this report and the general discussion of how MHC diversity evolved are of interest to at least some of the community. The manuscript is very well written and makes it easy to follow the theoretical and methodological details of the model and the arguments. I have only a few minor comments that I am detailing below. Furthermore, I would be very interested to read a response by the previous authors, especially on the relevance of this scaling/weighing factor that they introduced into their model, as it is possible that I might have missed something about its meaning.

  12. eLife

    Author response:

    Reviewer #1 (Public review):

    It appears obvious that with no or a little fitness penalty, it becomes beneficial to have MHC-coding genes specific to each pathogen. A more thorough study that takes into account a realistic (most probably non-linear in gene number) fitness penalty, various numbers of pathogens that could grossly exceed the self-consistent fitness limit on the number of MHC genes, etc, could be more informative.

    The reviewer seems to be referring to the cost of excessively high presentation breadth. Such a cost is irrelevant to the inferior fitness of a polymorphic population with heterozygote advantage compared to a monomorphic population with merely doubled gene copy number. It is relevant to the possibility of a fitness valley separating these two states, but this issue is addressed explicitly in the manuscript.

    An addition or removal of one of the pathogens is reported to affect "the maximum condition", a key ecological characteristic of the model, by an enormous factor 10^43, naturally breaking down all the estimates and conclusions made in [RS]. This observation is not substantiated by any formulas, recipes for how to compute this number numerically, or other details, and is presented just as a self-standing number in the text.

    It is encouraging that the reviewer agrees that this observation, if correct, would cast doubt on the conclusions of Siljestam and Rueffler. I would add that it is not the enormity of this factor per se that invalidates those conclusions, but the fact that the automatic compensatory adjustment of cmax conceals the true effects of removing a pathogen, which are quite large.

    I am not sure why the reviewer doubts that this observation is correct. The factor of 2.7∙1043 was determined in a straightforward manner in the course of simulating the symmetric Gaussian model of Siljestam and Rueffler with the specified parameter values. A simple way to determine this number is to have the simulation code print the value to which cmax is set, or would be set, by the procedure of Siljestam and Rueffler for different parameter values. In another section of this response I will describe how to do this with the simulation code written and used by Siljestam and Rueffler; doing so confirms the value that I obtained with my own code. Furthermore, I will now give a theoretical derivation of this factor.

    As specified by Siljestam and Rueffler, the positions of the m pathogens in (m-1)-dimensional antigenic space correspond to the vertices of a regular simplex centered at the origin, with distance between vertices equal to 1. The squared distance from the origin to each of the m vertices of such a simplex is (m-1)/2m (https://polytope.miraheze.org/wiki/Simplex). Thus, the sum of the m squared distances is (m-1)/2. For the (0, 0) homozygote, condition is multiplied by a factor of exp(-(vr)2/2) for each pathogen, where r is the distance from the origin. It follows that, with v=20, all the pathogens together decrease condition by a factor of exp(202∙(m-1)/4) = exp(100∙(m-1)). Thus, increasing or decreasing m by 1 changes this value by a factor of exp(100) = 2.7∙1043.

    This begs the conclusion that the branching remains robust to changes in c_max that span 4 decades as well.

    That shows only that the results are not extremely sensitive to cmax or K. They are, nonetheless, exquisitely sensitive to m and v. This difference in sensitivities is the reason that a relatively small change to m leads to such a large compensatory change in cmax a change large enough to have a major effect on the results.

    As I wrote above, there is no explanation behind this number, so I can only guess that such a number is created by the removal or addition of a pathogen that is very far away from the other pathogens. Very far in this context means being separated in the x-space by a much greater distance than 1/\nu, the width of the pathogens' gaussians. Once again, I am not totally sure if this was the case, but if it were, some basic notions of how models are set up were broken. It appears very strange that nothing is said in the manuscript about the spatial distribution of the pathogens, which is crucial to their effects on the condition c.

    I did not explicitly describe the distribution of pathogens in antigenic space because it is exactly the same as in Siljestam and Rueffler, Fig. 4: the vertices of a regular simplex, centered at the origin, with unity edge length.

    The number in question (2.7∙1043) pertains to the Gaussian model with v=20. As specified by Siljestam and Rueffler, each pathogen lies at a distance of 1 from every other pathogen, so the distance of any pathogen from the others is indeed much greater than 1/v. This condition holds, however, for most of the parameter space explored by Siljestam and Rueffler (their Fig. 4), and for all of the parameter space that seemingly supports their conclusions. Thus, if this condition indicates that “basic notions of how models are set up were broken”, they must have been broken by Siljestam and Rueffler.

    Overall, I strongly suspect that an unfortunately poor setup of the model reported in the manuscript has led to the conclusions that dispute the much better-substantiated claims made in [SD].

    The reviewer seems to be suggesting that my simulations are somehow flawed and my conclusions unreliable. I will therefore describe how my conclusions about sensitivity to parameter values can be verified using the simulation code provided by Siljestam and Rueffler themselves, with only small, easily understood modifications. I will consider adding this description as a supplement when I revise the manuscript.

    The starting point is the Matlab file MHC_sim_Dryad.m, available at https://doi.org/10.5061/dryad.69p8cz98j. First, we can add a line that prints the value of the variable logcmax, which represents the natural logarithm of cmax determined and used by the code. Below line 116 (‘prework’), add the line ‘logcmax’ (with no semicolon).

    Now, at the Matlab prompt, execute MHC_sim_Dryad(false, 8, 20, 1) to run the simulation for the Gaussian model with m=8, v=20, and K=1. The output will indicate that logcmax=700, in accord with the theoretical factor exp(100*(m-1)) derived above. The allelic diversity, ne, will rise to a steady state-level of about 140, as in the red curve of my Fig. 2.

    Now lower m to 7, i.e, run MHC_sim_Dryad(false, 7, 20, 1). The output will indicate that logcmax=600. This confirms that lowering m by 1 causes the code to lower the value of cmax by a factor exp(100)=2.7∙1043, which must also be the factor by which the condition of the most fit homozygote would increase without this adjustment.

    With the change of m to 7 and the compensatory change in cmax, steady-state allelic diversity remains high. But what if m changes but cmax remains the same, as it would in reality?

    To find out, we can fix the value of cmax to the value used with m=8 by adding the following line below the line previously added: ‘logcmax = 700’. With this additional modification in place, executing MHC_sim_Dryad(false, 7, 20, 1) confirms that without a compensatory change to cmax, lowering m from 8 to 7 mostly eliminates allelic diversity, in accord with the corresponding curve in my Fig. 2. Similarly, raising m from 8 to 9, or changing v from 20 to 19.5 or 20.5 (executing MHC_sim_Dryad(false, 8, 19.5, 1) or MHC_sim_Dryad(false, 8, 20.5, 1)), largely eliminates diversity, confirming the other results in my Fig. 2. Results for the bitstring model can also be confirmed, though this requires additional changes to the code.

    Thus, the extreme sensitivity of the results of Siljestam and Rueffler to parameter values can be verified with the code that they used for their simulations, indicating that my conclusions are not consequences of my having done a “poor setup of the model”.

    Response to Reviewer #2 (Public review):

    (1) The statement that the model outcome of Siljestam and Rueffler is very sensitive to parameter values is, in this form, not correct. The sensitivity is only visible once a strong assumption by Siljestam and Rueffler is removed. This assumption is questionable, and it is well explained in the manuscript by J. Cherry why it should not be used. This may be seen as a subtle difference, but I think it is important to pin done the exact nature of the problem (see, for example, the abstract, where this is presented in a misleading way).

    I appreciate the distinction, and the importance of clearly specifying the nature of the problem. However, Siljestam and Rueffler do not invoke the implausible assumption that changes to the number of pathogens or their virulence will be accompanied by compensatory changes to cmax. Rather, they describe the adjustment of cmax (Appendix 7) as a “helpful” standardization that applies “without loss of generality”. Indeed, my low-diversity results could be obtained, despite such adjustment, by combining the small change to m or v with a very large change to K (e.g., a factor of 2.7∙1043). In this sense there is no loss of generality, but the automatic adjustment of cmax obscures the extreme sensitivity of the results to m and v.

    (2) The title of the study is very catchy, but it needs to be explained better in the text.

    I had hoped that the final paragraph of the Discussion would make the basis for the title clear. I will consider whether this can be clarified in a revision.

  13. Version published to 10.1101/2025.05.27.656382 on bioRxiv