Major histocompatibility complex (MHC) fragment numbers alone – in Atlantic cod and in general - do not represent functional variability

The correspondence by Dijkstra & Grimholt¹ provides critical concerns about a publication by Malmstrøm et al in Nature Genetics in October 2016², concluding that their new model of MHC class I evolution, reflected in their title “Evolution of the immune system influences speciation rates in teleost fish”, is unsubstantiated. I concur with their three main criticisms “(1) that the authors did not provide sufficient evidence for presence or absence of intact functional MHC class I or MHC class II system genes, (2) that they did not discuss that an MHC subpopulation gene number alone is a very incomplete measure of MHC variance, and (3) that the MHC system is more likely to reduce speciation rates than to enhance them.”
 
All three critical points are well founded and stand-alone without too much need for further support. However, I have added the following 14-point commentary for Dijkstra & Grimholt¹, Malmstrøm et al (2016)² and others to consider and elaborate on if they would like to because they are important concerns in the field of MHC genomics, biological function and evolution.
 
1. According to Dijkstra & Grimholt¹ “the MHC counting strategy by Malmstrøm et al.² should be deemed incomplete and far too simplistic. For their number determination Malmstrøm et al.² solely relied on estimation of U plus Z lineage genomic α3 exon fragment numbers, despite that the typical “birth and death” mode of MHC evolution can produce many pseudogenes.”  I agree that this is a major problem with the Malmstrøm et al (2016)² paper, one that also omits the other categories of the MHC I such as the L, S and P lineages that might contribute to a much large number of MHC I-like genes. In this regard, it seems that Malmstrøm et al.² have taken only the genomic exon fragment numbers of the MHC I U and Z lineages to represent the entire immune system of their title.
 
2. Dijkstra & Grimholt¹are also right to point out that low coverage sequencing by next generation sequencing (NGS) can result in the artifactual loss of genes that in turn can lead to misleading or incorrect conclusions when counting for gene copy numbers or looking for gene losses. Malmstrøm et al (2016)² sequencing coverage was 9 to 39x and they recovered only about 75% of the conserved eukaryotic genes. Therefore, in this situation, the de novo low coverage sequencing data should not have been used as evidence for the absence of genes from the genome without providing a properly organized high coverage map of genomic assemblies to show where the sequences are missing in the genomes. The reviewers and editors of the Malmstrøm et al (2016)²paper should have been aware of this basic problem of using low coverage NGS, particularly with respect to looking for a few needles in a haystack. For a better understanding of the advantages and disadvantages of MHC genotyping and haplotyping by NGS see the review by Shiina et al (2016)³.
 
3. Malmstrøm et al (2016)² reported that there was no overall correlation between the combined MHC copy numbers and the HOX gene copy numbers that were used as a control. As already pointed out by Dijkstra & Grimholt¹, the U and Z genes should have been analysed separately. Nevertheless, it would have been interesting to see how these duplicated HOX development genes, which also have been implicated in driving speciation, compared with the properly classified duplicated MHC class I adaptive immune genes (separate Z, X, L, S and P lineages) at the classical and non-classical level in Fig 3 during speciation rate simulations².  In addition, it appears from Fig S3 that Malmstrøm et al (2016)²  might have missed an inverse relationship between MHC & HOX for the MHC copy numbers up to 25 and direct correlation between MHC & HOX for the MHC copy number from 25 to 50.
 
4. It seems absurd to count up only the short sequences of a3 fragments from a low coverage sequence library and extrapolate the numbers counted in Fig 2b to reconstruct an artificial model for duplicated MHC gene copies influencing speciation or evolution without first knowing their categories (classical versus non-classical presentations), functions, overall structure and coding ability, transcriptional activity and genomic locations.  Malmstrøm et al (2016)² provided no properly organized genomic assemblies or genomic gene maps and no information about genomic distribution of the MHC I or II sequenced fragments or the duplication mechanisms involved. If they had done so, they might have added important information to better assess and place the threshold MHC I copy numbers and gene distributions into some sort of genomic perspective⁴. More reliable models for the evolution of MHC class I genomic duplications might be achieved by providing duplication gene maps and the phylogenic relationships of the duplicated gene sequences showing likely duplication mechanisms, where and how these genes are located relative to each other, and how the genomic structures have changed in a comparative analysis between species. See Kulski et al. (2004)⁵ for an example of one such duplication model. Mapping with phylogeny is a more informative approach than just constructing phylogenetic trees using one or more single exonic sequences from a limited number of each species and then claiming that the changes in copy numbers influence the speciation rates of almost the entire number of extant fish species. Perhaps, the Malmstrøm et al (2016)² low coverage sequence libraries could still be used effectively to reconstruct full-length gene structures and targeted genomic regions that harbour multiple copies of the MH(C) genes in a comparative analysis.
 
5. The multiple-optima Ornstein Uhlenbeck (OU) model⁶. 
(a) According to Malmstrøm et al (2016)² the multiple-optima OU model vastly outperformed alternatives such as Brownian motion, white noise, single-peaked OU and early-burst models. This find­ing corroborated their hypothesis that MHC I copy number evolution is characterized by selection toward intermediate optima, resulting from a tradeoff between detection and elimination of pathogens. Presumably, the authors preferred the OU multiple-optima model to the other models because it supported rather than falsified their hypothesis. Of course, this highly artificial computing model did not detect a tradeoff between detection and elimination of pathogens, this would be the authors’ own biological hypothesis and bias.  The OU model intrinsically sets optima (biases) according to its built-in algorithm⁶, and this is one of the main objections to the use of this prediction model. The OU model artificially generates bias because its purpose is to find the trait optimum that stabilizes selection⁶. The misguided conclusion by Malmstrøm et al (2016)² in using the OU model is that the trait optimum influences speciation. 
 
(b) Interestingly, Malmstrøm et al (2016)² did not directly test the opposing hypothesis that the MHC I copy number evolution is not characterized by selection toward intermediate optima, and does not result from a tradeoff between detection and elimination of pathogens. Possibly, their best control in this regard was the simple Brownian model that did not work as well for them as the multiple-optima OU model that has extra parameters such as the addition of an overall optimum trait value to which all lineages are attracted. Other evolutionists often prefer the Brownian model for the reason that it is a simple, neutral model without the added bias of creating optimum trait value.
 
( c ) The OU multiple-optima model is not a fool-proof algorithm, and a number of evolutionists believe that it can be an unreliable or misleading model. According to Cooper et al (2016)⁷, although widely used, the properties of the OU model, and other direct extensions of the Brownian model, are poorly understood leading to the potential for inappropriate use and misinterpretation of results. In particular, Cooper et al (2016)⁷ used computer simulation studies to demonstrate that the single peaked OU model error rates are unacceptably high when tree size is small (< 200 species tips), when likelihood ratio tests or Akaike information criterion (AICc) are used to select the best model, and when measurement errors are introduced into the data. They also showed that when the alpha parameter of trait evolution was extremely small (<1) in the OU model it was indistinguishable from Brownian motion, and as the alpha value became larger it favoured OU prediction models, until the larger values of alpha were indistinguishable from white noise and it was therefore independent of phylogeny. The alpha values for the Malmstrøm et al (2016)² model selection analysis were markedly less than one (Supplementary Table 13), suggesting that they could have accepted the Brownian model over the OU model as the better model fit. 
 
6. The BiSSE threshold model. 
Malmstrøm et al (2016)² carried out binary state speciation and extinction (BiSSE) analysis to estimate differences in diversification rate between lineages with high and low MHC I copy numbers. They found that diver­sification rates based on correlation estimates differed most when the threshold was placed between 20 and 25 copies (Fig. 3c). With a threshold in this range, the model with two separate speciation rates for lineages with high and low copy numbers was statistically better supported than a model with a single speciation rate parameter.  On this basis, they concluded that, ‘These results suggest that the influence of MHC I genes on speciation rates is stronger in species that have already evolved at least 20 copies.’ In comparison, the number of MHC I gene copy numbers in humans (excluding haplotype differences) is approximately 18 genes; 6 classical and non-classical MHC I genes, 5 CD1 genes, and 5 PHFZ genes (MICA, MICB, MR1, HFE, Zn-A2-GP, etc). Thus, in comparison to some fish species, humans are diversifying along very nicely as a ‘diversified’ species approaching the ‘magical’ threshold of between 20 and 25 MHC I copies. 
 
It is noteworthy that Maddison et al (2007)⁸ highlighted the following assumptions that need to be taken into account when using their BiSSE package. For the BiSSE model analysis none of the characters associated with speciation rates can be said to be causing or influencing evolution, even if Maddison et al (2007)⁸ write, ‘the correct conclusion given a significant result using our method is that the character examined or a codistributed character appear to be controlling diversification rates.’ At best, the binary character state is an association, at worse a misleading one. Maddison et al (2007)⁸ provided the following cautions and assumptions about the likelihood of the BiSSE (binary-state speciation and extinction) model:
 
(a) the transitions happen instantaneously over the time scales considered (i.e., ignore periods of time during which a species is polymorphic). 
(b) these events are independent of one another; in particular, the character state change does not, in and of itself, cause speciation (or vice versa).
(c) an accurate rooted phylogenetic tree with branch lengths is known (the "inferred tree") and the character state is known for each of the terminal taxa. 
(d) the tree is assumed complete: all extant species in the group have been found and included.  
(e) all terminal taxa are contemporaneous, and the tree is ultrametric (i.e., the total root-to-tip distance is the same for all tips). 
 
I’m not convinced that Malmstrøm et al (2016)² considered or accepted these constraining assumptions when using BiSSE modeling.
 
7. Speciation and diversification rates. 
(a) In the study by Santini et al (2009)⁹, the speciation rates within the Percomorpha clade were calculated to be at least ten times greater than in the Gadiformes order. Yet, according to Malmstrøm et al (2016)², there were fewer than 20 copies of the U genes for each of the 5 species in the Perciformes clade, compared to more than 20 copies and up to 100 copies of the U genes in 16 of 30 species in the Gadiformes order (Fig 2b, Malmstrøm et al 2016²). The use of only 5 species in the Perciformes order is a gross underestimate of the ten thousand or more species found in that order⁹. Moreover, in the Gadiformes order, there were closely related species (n=16) with > 20 U genes and different groups of closely related species (n=14) with < 20 U genes. Again, the number of species that were sequenced in the Gadiformes order is grossly under-represented. The Gadiformes comprises 10 families and more than 600 species⁹, whereas Malmstrøm et al (2016)² sequenced only 27, i.e.,  27 x 100/600 = 2,700/600 = 4.7% analysis of extant sequence, a percentage that is simply not good enough to support their extravagant conclusions.  Is < 5% of the 600 species really representative of the Gadiformes. Malmstrøm et al (2016)² have to be more temperate with their conclusions using such a small representative sample.There are clear inconsistencies with MHC I a3 fragment copy numbers in the Gadiformes order. The MHC I a3 frag copy numbers are low (<20) for Moridae and M. occidentalis in Macrourinae, and for Phycinae, Lotinae, and three species in Gadinae. Five species of Gadinae have between 20 and 40 copies. On the other hand, Bregmacerotidae, Merlucciidae, Melanonidae, Muraenolepidodae, and Trachyrincinae have MHC I a3 fragment copy numbers between 50 and 100. The threshold levels (20 to 25) are all over the place. Moreover, the lineage, genomic block duplication and hitchhiking (linkage) effects on MHC gene duplications (8 to >100 copies) in the Gadiformes have not been taken into account in the analysis of speciation rates (Fig 3, Malmstrøm et al  2016²), and, therefore, make the entire analysis unreliable.
 
(b) “Diversification rate analyses were calculated on the basis of the time-calibrated phylogeny and counts of species richness in each of the 37 mutually exclusive clades of teleost fishes"²; mostly from the Gadiformes order. The MHC I speciation model of Malmstrøm et al (2016)² appears contradictory for the Perciformes (10,033 species) that have speciation rates 18 times greater than Gadiformes (555 species)⁹ and yet the MHC I a3 fragment copy numbers are at least two times lower in Perciformes than Gadiformes (Fig 2b, Malmstrøm et al  2016²). Also, the Anabantiformes have 252 species – a species rate 40 times lower than Gadiformes⁹ and yet their MHC I a3 fragment copy numbers are at least two times higher than in Perciformes.
 
(c) Considering that there are more than 29,000 species of teleost fishes⁹, a highly limited analysis by Malmstrøm et al (2016)² of using a sample group of less than 0.2% of the available extant species cannot be considered to be statistically, taxonomically or biologically significant or sufficiently reliable to conclude that, “Evolution of the immune system influences speciation rates in teleost fish”². What does a species half-life of 25 million years mean in the context of 29,000 species of teleost fishes? If the multiple-regime OU model is wrong, highly biased or misinterpreted then does it validate or support the overall hypothesis of Malmstrøm et al (2016)²? Also, what does an optimal trait actually mean in the context of 29,000 species? If, a suboptimal number of MHC I copies are detrimental to a species, then how have divergent species managed to survive for so long with a half-life of 25 million years of adaptation?  Also, if, as Malmstrøm et al (2016)² say, ‘Such gene family expansions may promote biological diversification by introducing new raw genetic material, potentially resulting in sub- or neofunctionalization and thus novel immunological pathways.’, then which of the non-optimal (greater than or less than the threshold of 20 -25 copies) MHC I genes are detrimental to the species? In this regard, there must be a gradation of functionally good and bad MHC I genes as their copy numbers approach the threshold (optimally good) and then deteriorate beyond it. Is this assumption of a MHC I copy number functional trait value as a quantitative marker of speciation at all testable?
 
8. In their discussion, Malmstrøm et al (2016)²  referenced the hypothesis of T cell depletion and hybrid fitness by Eizaguirre et al (2009)¹⁰ and concluded that,"Our analyses identify this threshold at 20–25 MHC I copies, suggesting that the effect of T cell depletion on hybrid fitness becomes more pronounced in this range and that this might affect mate choice in species with copy numbers above this threshold, promoting inbreeding and reinforcement."  Eizaguirre et al (2009)¹⁰ suggested that, “Super-optimal individual MHC diversity should be a common disadvantage for species hybrids in vertebrates, resulting in elevated parasite loads.”  In this regard, if high copy numbers of the MHC class I genes leads to hybridization and loss of the immune system as inferred by Eizaguirre et al (2009)¹⁰, then this more than likely would lead to extinction of populations and species. Extinction would be the most extreme and bizarre form of immune system influence on speciation rates. Furthermore, it is extremely speculative for Malmstrøm et al (2016)² to say that high copy numbers of the MHC class I genes with copy numbers above the threshold of 20 to 25 copies, promotes inbreeding and reinforcement, because, in fact, there is no such evidence for it. A more reasonable hypothesis is that high copy numbers of linked MHC class I genes, such as in rhesus macaque, or the mouse¹¹, or the cod¹², might benefit the species to better adapt to microbial inhabitants in a greater variety of geographical environments, although the evidence for this is tenuous as well. Despite ongoing debates, the selective advantage of MHC diversity in host-pathogen coevolution might not be easily resolved (at the macroevolution level) because of the constant number of insults by large numbers of pathogens in the life-time of an individual organism, population or a species and the arms race or Red Queen effect. Studies on extant species will always discover an example of a pathogen associated with a polymorphic MHC gene that might favour selective advantage for host-pathogen coevolution, whereas the pathogen that caused the extinction of a species is rarely or never found. To conclude that the immune system (that is, different copy numbers of the class I MHC genes²) influences speciation rates, it would have to be shown that the immunity gene products can commonly create reproductive barriers or genetic incompatibilities among populations that permit the maintenance of the genetic and phenotypic distinctive­ness of these populations in geographical proximity¹³; and this was not done².

9. Malmstrøm et al (2016)² did not provide any reliable evidence to support their speculation that evolution of the immune system influences speciation rates in teleost fish or that increas­ing MHC I diversity facilitates speciation⁷.  Instead, Malmstrøm et al (2016)² used their limited data and analyses using speculative models to jump to highly unsupported conclusions and quickly position the cart before the horse. Dijkstra & Grimholt¹ pointed out that the Malmstrøm et al (2016)² title “Evolution of the immune system influences speciation rates in teleost fish”, is unsubstantiated, and that their hypothesis seems to be “the wrong way around”. It should have been, “Speciation (rate) influences the evolution of the immune system in teleost fish.” Or, “Speciation rates are associated with diversity of MHC class I genes in teleost fishes”, which perhaps is too obvious and underwhelming. This is not simply the chicken or the egg causality dilemma; in fact, the change in title is better supported by the literature and the established theories of MHC genomic evolution in vertebrates⁴. However, because it is less “sexy” and controversial than the original title, it might not have been so readily published. 

10. A large number and variety of genome-wide gene duplications have been associated with speciation¹³, that is, genomic gene duplications are not limited to only class I MHC genes. If MHC I gene duplications effect or affect speciation, how do the other hundreds of gene duplications contribute to speciation rates? Also, do sequence variants or mutations in non-duplicated genes have any influence on speciation rates? It seems absurd to pick on only one group of gene duplications (e.g., MHC class I genes²) as those that are responsible for speciation and ignore all the others as an inconvenience. For example, a relatively recent comparative genomic study revealed how genomes change with speciation in an examination of genomes from five cichlid fish species, an ancestral lineage from the Nile, and four species from the East Africa lakes, Tanganyika, Malawi, and Victoria¹⁴. Compared to the ancestral Nile lineage, the East African cichlid genomes had many alterations in regulatory elements, accelerated evolution of protein-coding elements in genes for pigmentation, an excess of gene duplications, and other distinct features that affect gene expression associated with transposable element insertions and novel microRNA. Each species also contains a reservoir of mutations different from the other species¹⁴. Much of the diversity between the cichlid fish species evolved in a nonparallel manner often rapidly due to sexual selection and genetic conflicts between males and females or between different regions of the genome at a regulatory level¹⁴ rather than by the slower and weaker forces of classical natural selection¹³. If sexual selection and genetic conflict at the genomic regulatory level are the prime movers of speciation rate, it is difficult to conclude that the variable diversity of a few MHC gene copies are responsible for speciation as well as for the many other associated genomic changes associated with speciation.
 
11. Malmstrøm et al (2016)² informed us in the introduction section of their publication that "Our results highlight the plasticity of the vertebrate adaptive immune system and support the role of MHC genes as ‘speciation genes’, promoting rapid diversification in teleost fishes."  MHC class I gene copy number variability occurs across many different species, families, orders and domains. Because there is such enormous variability in MHC class I gene copy number for hundreds⁴ or possibly even thousands of different chordate species, it is not possible to conclude meaningfully that the expansion of MHC class I genes provides an undefined advantage of one species over another. For example, the great apes (humans, chimpanzees, gorillas and orangutans) have about six functional MHC class I genes, whereas the old and new world monkeys often have up to 15 or more^4,11. Is this evidence that the MHC class I genes influence the rate of speciation in primates? And if so, what does that really mean in the whole scheme of things? How do the species with low copy numbers of MHC class I genes survive so well over millions of years without the presence of another 90 to 100 copies of MHC class I genes? This question is often neglected, and yet it is important for a better understanding of the function and evolution of MHC genes between and within the vertebrate species. 
 
12. Taxonomic and lineage markers.
Mutations, indels and duplications drive diversity and evolution. However, most mutated genes within species and their families do not create or influence speciation rates in the sense that Malmstrøm et al (2016)² use the term, ‘speciation genes’. In comparative genomics and their sequence relationships between different species, most genomic sequences range between newly derived genes and the ultraconserved or the essential core coding and noncoding genes with varying amounts of sequence differences. Some genic and nongenic sequences such as the MHC genes and retrotransposons are highly polymorphic and therefore are useful taxonomic markers at the individual, population, species and broader lineage levels. The MHC gene sequences clearly are one of these useful taxonomic or lineage markers along with olfactory receptors, immunoglobulins, globins, HOX, TOLLs, KIRs, mitochondrial DNA, ribosomal RNA sequences and thousands of others that can be used comparatively in the phylogeny to undertake an examination of the accuracy and reliability of current taxonomical rankings and sister lineages. However, because many thousands of coding and noncoding genes (or sequences) are variants (polymorphic) or vary in copy numbers, we do not immediately or easily imply that all or some of them are responsible for speciation without providing further concrete evidence. This kind of extrapolation without the burden of proof is absurd and wrong. Similarly, to say that the polymorphisms demonstrate natural selection as if natural selection was a biological or molecular mechanism is meaningless without showing experimentally how these polymorphisms benefit or disadvantage the organism over all the other different polymorphisms. 
 
13. On the basis of either a priori or a posteriori reasoning, the immune system obviously affects the wellbeing of individuals and populations, but whether it can be extrapolated to speciation events and speciation rates remains highly dubious and most probably unlikely. It seems too farfetched to blame MHC class I genes with high copy numbers over threshold levels of promoting inbreeding and reinforcement² because this in turn could create hybrid inviability or sterility resulting in postzygotic isolation. Although the population conditions in many models of rapid speciation do favour inbreeding and/or hybridization¹³, none of the teleost species tested by Malmstrøm et al (2016)² were shown to be either inbreeding or in postzygotic isolation. The factors responsible for either prezygotic or postzygotic isolation are likely to be independent of the adaptive immune system, although zealots might argue otherwise. Hybridization between diverging lineages in post-zygotic reproductive isolation can trigger genome instability. For most animals without an adaptive immune system and for plants without a MHC, speciation depends on the shrinkage, expansion and equilibrium (e.g., aneuploidization and dysploidy) of the genome and the containment and functionality of all the essential genomic information to develop an optimal balance between stability and plasticity within the organism in order to first survive and then propagate and expand itself as a new species¹³. In those rare and ‘traumatic’ transitional situations, there is no need for particular ‘speciation’ genes such as variable copies of the class I MHC genes to influence speciation. The rarely observed transition from population ‘trauma’ to a new speciation event depends on an array of totally different factors for creating postzygotic isolation events including interbreeding between semi-isolated populations and an elaboration on the existence of stress-induced changes in chromosomal and ploidy integrity both in hostile and non-hostile environments.
 
14. Finally, Malmstrom et al (2016)² admirably sequenced 66 teleost species by a next generation sequencing method and identified an array of MHC I and MHC II exonic fragments for phylogenetic and speciation analysis using the multiple-regime OU model to predict the optimal MHC I copy number as an evolutionary trait optimum affecting speciation. However, the conclusions of the paper by Malmstrom et al (2016)² especially for the MHC I gene copy numbers are unreliable because they are based on far too many assumptions, speculations, contradictions, incomplete or missing data and unproven predictive models with little or no empirical evidence in support. Nevertheless, their simple, but controversial hypothesis is published, and now it is up to them and others to test its validity and "consider plausible alternative hypotheses in a firm hypothesis-testing framework in which alternative hypotheses make clear [and sensible] predictions of emerging patterns that can be unambiguously associated with particular models."⁷

The critique of Malstrøm et al. presented here makes some very valid points that are well supported by the literature.

It has long been true in fish MHC research that the fact that a gene has not been reported to be present in a particular species does not mean that is not present. Modern genomics techniques has presented better proof for this assertion, with the lack of an MHCII/CD4 pathway in gadids being the most prominent example, but even modern genomics techniques are not iron clad 100% proof and should be checked very carefully before definitive statements are made. Thus the comments about verifying the presence or absence of specific genes in the numerous  species by other means is valid.

Additionally, the treatment of all U and Z genes as identical units while ignoring the allelic diversity of each gene within those classes is indeed a serious flaw in the reasoning of Malstrøm et al. There is significant variability in diversity U gene families which will have differing effects on T cell selection that simply counting gene numbers will not address.

Dijkstra and Grimholt's critique should be carefully read and addressed.