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

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].

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.

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.

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 c_max 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∙10⁴³ 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 c_max 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) 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(20²∙(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∙10⁴³.

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 c_max 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 c_max 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∙10⁴³) 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, n_e, 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 c_max by a factor exp(100)=2.7∙10⁴³, 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 c_max, steady-state allelic diversity remains high. But what if m changes but c_max remains the same, as it would in reality?

To find out, we can fix the value of c_max 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 c_max, 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 c_max. Rather, they describe the adjustment of c_max (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∙10⁴³). In this sense there is no loss of generality, but the automatic adjustment of c_max 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.