BACKGROUND
the patterns and causes of genotypic variation in human genes have been a focus of great recent interest in evolutionary biology.
different processes such as natural selection, genetic recombination, genetic drift, demography, as well as physicochemical properties of cells, can influence this diversity.
various methods have been devised to represent and quantify genetic variation and to detect its causes .
here we use a novel approach based on genotype networks to represent and analyze genetic variation in human genes.
genotype networks are graphs that consist of nodes, which correspond to genotypes with the same phenotype, where sameness can be defined as narrowly as enzyme activity, or as broadly as viability.
nodes that differ in some minimal way from each other are adjacent, i.e., connected by an edge in such a graph.
the genotypes we consider are haploid genotypes  of human genes in a sample of the human population, and we call two genotypes adjacent if they differ in a single nucleotide.
genotype networks can be useful to address various evolutionary questions, such as how novel adaptations originate, and what role phenotypic robustness or plasticity play in adaptation  <cit> .
in the past, they have been mostly built from computational models of genotype-phenotype maps , but high-throughput genotyping allows us to build genotype networks from experimental data  <cit> .
representing such data in the form of a network makes the large analytical toolbox of graph theory available, which has been useful in fields as different as ecology and the social sciences .
a common form of representation for genetic variation data is the phylogenetic tree, which shows the evolutionary relationship among a set of taxa, individuals, or genes that constitute the leaves of the tree.
the common ancestors of these taxa form the interior nodes of such a tree.
in a gene tree, these ancestors can be reconstructed with the help of probabilistic models of sequence evolution .
phylogenetic trees are by definition acyclic graphs: they do not contain cycles – paths of edges that start from a node, pass through other nodes, and return to the same node.
the acyclic nature of phylogenetic trees implies one major limitation of such trees: they cannot easily accommodate evolutionary genealogies more complex than simple vertical descent with modification .
such genealogies can lead to reticulate networks of phylogenetic relationships.
thus, multiple mechanisms to create genetic diversity, such as hybridization, allopolyploidization, sexual reproduction, recombination, gene conversion, and homoplasy, which lead to mosaic patterns of relationships among nodes are not easily accommodated in tree-like structures.
genotype networks provide information complementary to phylogenetic trees that are not subject to this limitation, because they can accommodate cycles.
figure  <dig> shows a short cycle in a hypothetical genotype network involving four dna sequences.
edges reflect adjacent genotypes that differ in a single nucleotide.
assume, for example, that genotype  <dig> is ancestral to the other genotypes, and different substitutions  produce genotypes  <dig> and  <dig> from it.
genotype  <dig> then experiences an additional a10t substitution that creates genotype  <dig>  this mutational path leads to a closed cycle, where three of the four edges reflect a substitution event.
the fourth edge is a consequence of the first three events, because they render genotype  <dig> adjacent to genotype  <dig>  similar scenarios can be developed if a genotype different from genotype  <dig> is ancestral.
regardless of this ancestor, cycles require sequence changes that render the descendants of one  genotypes more similar rather than less similar.
in other words, cycles require some form of homoplasy, i.e., parallel or convergent evolution .
more generally, homoplasy is said to exist when two lineages display the same genetic or phenotypic characters, even though this similarity has not arisen through common ancestry  <cit> .fig.
<dig> a hypothetical example of a four-node cycle in a haplotype network.
the example indicates a hypothetical dna sequence where two nucleotide changes occur at position  <dig> and  <dig>  circles  correspond to genotypes.
an edge connects two nodes if they differ by a single mutation.
lettering next to each node indicates the nucleotides at which two genotypes differ.
edge labels show changes required to create a genotype from its neighbor, e.g., “a20g” indicates a change from a to g at position  <dig> of the hypothetical sequence.
in the example, mutations at positions  <dig> and  <dig> create genotypes  <dig> and 3 from genotype  <dig> .
then, either genotype  <dig> mutates at position  <dig> from a to g, or genotype  <dig> mutates at position  <dig> from a to t, or both of these mutations take place, and create genotype 4
homoplastic sequence evolution has been documented in a wide variety of molecules .
it can be caused by chance alone, which is expected to be rare in long evolving biopolymers with multiple kinds of monomers, because random mutations are more likely to cause such polymers to diverge than to converge.
mutational biases, strong selective constraints on sequence evolution  <cit> , positive selection , or genetic recombination  <cit>  can also cause homoplasy.
here we construct haploid genotype networks for each of  <dig>  genes in the human genome, based on single nucleotide variation data available for  <dig>  individuals from four continents  <cit> .
we analyze short cycles up to length eight in these networks, and find that the haploid genotype  networks of  <dig> genes show a significant excess of cycles that cannot be explained by chance alone.
after having excluded recombination as a prominent cause of these cycles, we focus on positive selection as a possible cause, and present evidence that in at least some of these genes positive selection may help explain the existence of cycles.
RESULTS
constructing and characterizing haplotype networks
to construct genotype networks for  <dig>  human individuals, we used haploid genotypes  with single nucleotide variation data available from the  <dig>  genomes project  <cit> .
thus, our genotype networks are haplotype networks, and from now on, we use the term haplotype network instead of genotype network.
for each human gene, we constructed one haplotype network.
two principal definitions of such a network are germane for this paper.
by the first definition, a haplotype network of a human gene is a graph whose nodes correspond to protein-coding dna sequences of the gene in different individuals.
two nodes  are adjacent if they differ in a single base pair .
by the second definition, two nodes are adjacent if their coding sequences differ by a single non-synonymous  change.
the second kind of network can thus also be viewed as a network of proteins or amino acid sequences, in which neighboring proteins differ in a single amino acid.
we first created both dna- and protein-based haplotype networks based on the above definition, collapsing those nodes with identical sequences into one .
networks can be made of one or more components.
each component is a subgraph in which any two nodes are connected to each other by a path of edges.
we found that the average size of the largest connected component – commonly referred to as the giant component – relative to total network size is significantly larger in protein-based networks  than in dna networks   .
because our statistical analyses focus on the giant component of each haplotype network and work best if this component comprises as many nodes as possible, we focus on protein-based haplotype networks for the rest of this contribution.
the  <dig>  genomes dataset we use contains information from  <dig>  genes, but we constructed haplotype networks only for those  <dig>  protein-coding sequences that showed at least one amino acid variant.fig.
<dig> haplotype networks vary greatly in structure among genes.
a haplotype network of the gene otog .
among all protein-based haplotype networks comprising more than  <dig> sequences, otog has the network with the largest giant component where all nodes fall into this component .
b haplotype network of the gene hla-b, which is the most fragmented network, with  <dig>  nodes in  <dig>  components.
circles in a) and b) correspond to different genotypes, while edges connect genotypes that differ by a single point mutation.
circle color and size correspond to the degree  of the node, where darker and larger nodes have a higher degree.
c number of components versus network size for dna-based  and protein-based haplotype networks .
circle size in c) corresponds to the relative size of the giant component within each haplotype network
figure 2a and b illustrate with two examples that haplotype networks for different genes can differ greatly in their topology.
the left network , derived from the gene otog, which encodes otogelin, comprises  <dig> nodes organized into a single component, whereas the right network, from gene hla-b  is highly fragmented and has  <dig>  nodes in  <dig>  different components
more generally, fig.
2c shows the distributions of the number of components and the size of the largest component for all genes we considered.
there are  <dig>  networks with only a single component, but most of these networks are small, with an average of  <dig>  sequences.
the network with the most components is the highly fragmented hla-b network with  <dig>  components.
hla-b is known to be under strong balancing  <cit>  and divergent selection  <cit> , which causes great genotypic diversity.
this diversity translates into high network fragmentation, i.e., a network with many components.
some haplotype networks have very large giant components with up to  <dig> sequences.
however, in most  networks, the largest component is very small, comprising a maximum of ten sequences.
the network with the largest giant component where all sequences fall into that component is that of otog .
cycles in haplotype networks
a cycle in a network is a series of edges starting from one node and ending with the same node, while passing other nodes along the cycle only once.
in haplotype networks constructed from biallelic gene variants, the simplest elementary cycle, i.e.
a cycle not decomposable into smaller cycles, is a square.
the reason is that cycles with an odd number of edges, e.g.
triangles or pentagons, are impossible when all snps are biallelic.
figure  <dig> shows a square that involves the mutation of a hypothetical dna sequence at two different sites .
next to each circle  the nucleotide residues at these positions are indicated, and along the edges, the specific nucleotide changes that occurred.
if genotype  <dig> is the most recent common ancestor of its neighbors, then its two neighbors have undergone two different mutations: specifically, genotype  <dig> has experienced a change from a to t at position  <dig> and genotype  <dig> has a change from a to g at position  <dig>  to form the single genotype  <dig> from its ancestors, i.e.
from either genotype  <dig> or  <dig>  either genotype  <dig> has undergone a change from a to g at position  <dig>  or genotype  <dig> has undergone a change from a to t at position  <dig>  so that the descendants of the two ancestral sequences  <dig> and  <dig> become not only more similar but identical to one another.
it is not necessary for both of sequences  <dig> and  <dig> to mutate to form genotype  <dig>  but a mutation in either of them can lead to the genotype and form a cycle.
regardless of whether genotype  <dig> or any other genotype is the common ancestor of the others, a square like this requires convergent sequence change.
in long biopolymer sequences with multiple monomers that evolve through random mutation alone, cycles should be rare, because it is unlikely that mutations become reversed to create sequences more similar to one another.
however cycles can be introduced by mutation biases that allow only certain residues to change, or by selection that causes only certain changes to survive, i.e., by evolutionary constraints.
another possibility is recombination, which might occur between genotypes  <dig> and  <dig>  which can result in genotype  <dig>  the same mechanisms can give rise to longer cycles .
figure 3a shows the distribution of the number of squares in all networks.
<dig>  of  <dig>  networks had no squares.
the network with the most squares is that of the gene dnah <dig> and contains  <dig>  squares.
the inset of fig.
3a shows the distribution of hexagons and octagons.
the networks with the largest number of hexagons  and octagons  are those of genes map2k <dig> and hla-b, respectively.
note the small numbers of hexagons and octagons compared to squares.
even though we enumerated elementary cycles up to length eight – beyond that, our analyses become computationally too costly – we focus most of the following analysis on squares, because they are by far the most abundant cycles.fig.
<dig> distribution of the number of cycles in all networks and in networks with an excess of squares.
a distribution of the number of elementary squares, as well as elementary hexagons and octagons  in all protein-based networks.
b distribution of the number of squares in  <dig>  randomized networks derived from the giant component of the hla-b network, whose number of squares  is indicated by a red arrow.
c number of squares  in the largest components of the haplotype networks of  <dig> genes with significantly more squares than expected by chance alone, together with the mean number of squares  found in  <dig>  randomly generated networks for each gene.
shaded areas depict the maximum and minimum number of squares in the randomized networks.
note the logarithmic scale on the vertical axis.
unconstrained or constrained mutation cannot explain the large number of cycles in many networks
because some amount of homoplasy can occur by chance alone, we wished to determine whether all squares we observed could have occurred by chance homoplasy.
to this end, we created randomized haplotype networks in which the same amount of evolutionary change occurred as in the actual networks.
in our first randomization procedure, we created a  dna sequence of the same length as the coding sequence of a gene, and created a haplotype network from it by simulating a pattern of mutations designed to yield a network with the same number of edges  and the same distribution of degrees  as the actual network .
specifically, we compared the number of squares in each haplotype network to  <dig>  such randomly generated networks, and found  <dig>  genes whose actual number of cycles was greater than all of the  <dig>  randomly generated networks.
thus, based on this criterion there are  <dig>  genes whose total number of squares cannot be explained by chance homoplasy alone  at  <dig> ) .
one can argue that this procedure does not take into consideration the actual patterns of variation observed in the data, namely that only a small subsets of sites in any one gene have been subject to mutation, and that all of the sites are biallelic, that is, only two variant nucleotides occur in them.
both patterns arise from the fact that the human population sample is not highly diverged, and that natural selection constrained the evolution of these sequences, i.e., it eliminated some mutations that occurred in them.
we thus modified our randomization procedure to reflect these facts .
with these more conservative criteria, we still found  <dig> genes  whose haplotype networks have significantly more squares in their networks than expected by chance alone .
that is, their number of squares cannot be explained by mutational patterns and purifying selection alone.
figure 3b shows, as an example, the number of squares  in the haplotype network of hla-b, which is  <dig>  standard deviations greater than the mean number of cycles  in  <dig>  randomized networks .
figure 3c shows the number of squares in all  <dig> networks , together with the mean , minimum, and maximum  number of squares for  <dig>  randomized networks created for each of the  <dig> haplotype networks.table  <dig> genes with an excess of squares in their giant component
ttn
mki67
obscn
pkd1l2
muc16
muc17
igfn1
gpr98
prune2
syne2
ahnak2
hla-dpb1
alpk2
hla-c
flg
pramef2
hrnr
muc5b
pclo
hla-a
muc12
lama5
cyp2a7
hla-b
poted
neb
muc4
pkd1l1
fbn3
dchs2
fras1
plin4
eys
fcgbp
tg
ush2a
lilrb3
lilra6
dnah17
hla-drb1
dnah5
dnah11
the number of squares in these genes cannot be explained by random homoplasy or mutational constraints.
the middle column cites studies that provide evidence for positive selection, wherever such evidence is available.
after fdr correction, the p-value of the statistical test comparing the actual number of cycles against that in  <dig>  randomized networks  is  <dig>  for all genes
additional file 4: figure s <dig> shows the distribution of elementary cycles with length four, six and eight among the  <dig> genes with an excess of squares, and fig.
<dig> shows the proportion of the sequences that form part of a square in the largest connected component of each gene network.
for some genes, such as poted  all sequences form part of a square, and even for genes where the proportion of sequences in a square is low, such as hla-c  and ttn , it exceeds 40 % .fig.
<dig> proportion of sequences that are part of a cycle.
proportion  and actual number of sequences  that are part of a square in the giant connected components of haplotype networks for those  <dig> genes  with a significant excess of squares
we note that the  <dig> networks with an excess of squares are otherwise very heterogeneous in their properties.
they range from the network of mki <dig>  where only  <dig> nodes lie in the largest connected component, to the network of dnah <dig> , where  <dig> nodes do .
some of the networks have very few components, such as that of poted with a single component, whereas others have many components, such as the highly fragmented hla-b network with  <dig>  components .
even properties within the largest connected components are heterogeneous.
for example, in some networks, such as that of pkd1l <dig>  the distribution of the numbers of neighbors of each sequence is highly right-skewed and dominated by sequences with few neighbors, while in others it is more symmetric .
assortativity coefficients, which quantify the tendency of each node to connect to other nodes with a similar number of neighbors, also vary broadly.
some networks are assortative , whereas others are disassortative .
gene ontology  enrichment analysis on biological processes shows several immune system-related processes that are enriched in the  <dig> genes, namely “antigen processing and presentation of endogenous peptide antigens” and “interferon-gamma-mediated signaling pathway” .
go enrichment analysis of molecular functions reveals the two enriched functions “calcium ion binding” and “peptide antigen binding”.
“peptide antigen binding” is again associated with the immune system.
given the strong representation of hla genes among genes with an excess of cycles, we also asked how the go enrichment analysis would change if we excluded the hla genes.
we found a single enriched biological process, namely “o-glycan processing”, and two enriched molecular functions, namely “calcium ion binding” and “extracellular matrix constituent, lubricant activity”.
we also asked whether genes with an excess of squares preferentially occurred in specific kegg  <cit>  or reactome  <cit>  pathways.
six genes were preferentially associated with kegg pathways.
they include tg  and the genes in the hla family.
the enriched pathways comprise “epstein-barr virus infection”, “autoimmune thyroid disease”, “htlv-i infection”, “viral myocarditis”, “allograft rejection”, “phagosome”, “antigen processing and presentation”, “graft-versus-host disease”, “cell adhesion molecules ”, “herpes simplex infection”, and “type i diabetes mellitus”.
for reactome pathways, we found twelve genes enriched in six pathways.
the genes include those encoding mucins, the hla family and lilr family genes .
the enriched pathways are “termination of o-glycan biosynthesis”, “interferon gamma signaling”, “endosomal/vacuolar pathway”, “immunoregulatory interactions between a lymphoid and a non-lymphoid cell”, “antigen presentation: folding, assembly and peptide loading of class i mhc” and “defective galnt <dig> causes colorectal cancer  <dig> ”.
we note that both enriched kegg and reactome pathways include several immunity-related pathways.
recombination cannot account for an excess of squares in most networks
to exclude the possibility that genetic recombination may account for the excess of squares in some networks, we performed two complementary analyses.
first, we simulated for each gene the effect of recombination on haplotype network structure by creating haplotype networks based on a set of sequences that were subject to approximately as many recombination events as occurred in the human population since their common ancestry, as well as to as many mutations as there are edges in the network .
we repeated this process  <dig>  times for each gene, creating  <dig>  simulated haplotype networks, and counted the number of squares in them.
for each of the  <dig> genes, the empirical network showed more squares than each of the  <dig>  simulated networks .
in the second analysis, we asked whether gene conversion, a process of unidirectional recombination in which only one of the recombining sequences changes, may have caused the excess of squares  <cit> .
to this end, we used the program geneconv   <cit>  to detect gene conversion among the sequences in the giant components of the  <dig> haplotype networks.
we used sequences comprising both synonymous and non-synonymous changes to give the program more power in finding gene conversion events.
only one gene showed any sign of gene conversion, and it did so for only two of  <dig> sequences in cyp2a <dig> .
in sum, based on these analyses, it seems unlikely that recombination can explain the excess of squares we observe in the haplotype networks of  <dig> genes.
positive selection as a potential cause of squares
positive selection can be a driver for homoplastic or convergent evolution, where two separate lineages evolve the same trait independently  <cit> .
because such adaptive homoplasy can occur not only at the phenotypic level  <cit> , but also at the sequence level  <cit> , we wished to find out whether positive selection can help explain the excess of squares we observed in the haplotype networks of  <dig> genes.
previous studies had indeed indicated positive selection for at least  <dig> of the  <dig> genes  .
in addition, we used results from a branch-site likelihood test  <cit>  which indicates positive selection based on a ratio dn / ds >  <dig> observed along one or more branches of a phylogenetic tree.
this test has been applied to vertebrate genes in the selectome database  <cit> , which indicates that  <dig> of our  <dig> genes with abundant squares show patterns of positive selection, either in primates or in the bony vertebrates .
this number –  <dig> of  <dig> – is unlikely to be observed by chance alone ).
in addition, we used the xp-clr  test for neutrality  <cit>  .
the test compares different populations to identify rapid changes in a locus’ allele frequency that cannot be explained by random drift alone.
in applying this test, we used a test statistics  <cit>  pre-computed over 2 kb sliding windows that covered the human genome, and asked for each of our  <dig> genes whose haplotype network showed an excess of squares, whether two or more of the windows where the test-statistic indicated the action of positive section  overlapped with the gene .
by this criterion, six of our  <dig> genes showed evidence of positive selection in at least one population .
overall,  <dig> of our  <dig> genes with an excess of squares showed signs of positive selection by at least one of these criteria or by previous work.
we also analyzed patterns of synonymous and non-synonymous changes in more detail.
a commonly used indicator of positive selection for two protein-coding dna sequences is dn / ds, i.e.
the ratio of nonsynonymous changes dn per nonsynonymous site to synonymous changes ds per synonymous site.
values of dn / ds >  <dig> can indicate positive selection  <cit> .
unfortunately, dn / ds can be computed only for sequences more distantly related than those in our haplotype networks.
the reason is that in these networks, adjacent sequence pairs differ only in a single nonsynonymous mutation, and many adjacent pairs do not even show a single synonymous change.
more specifically, in the giant component of our networks, up to 80 % of sequence pairs do not show a single synonymous mutation , and this incidence of synonymous mutations is similarly low in the entire network.
moreover, it has been suggested that for very closely related sequences, dn / ds is not a sensitive indicator of positive selection  <cit> .
for these reasons, we compared the incidence of nonsynonymous and synonymous changes among groups of edges , reasoning that groups of edges with very few synonymous changes might provide hints that some or all members of the group may have been subject to positive selection.
most edges show no synonymous changes at all in some networks, which hints that positive selection may have played a role in creating their pattern of diversity .
we specifically compared edges with no synonymous change inside squares and outside squares.
while the fractions of edges without synonymous changes inside squares was not significantly different from those outside squares , the average number of synonymous changes on edges inside squares was significantly smaller than that outside squares for 14 %  of the genes .
figure  <dig> shows the average number of synonymous changes per edge for edges inside squares divided by that for edges outside squares.
for genes where this ratio is below  <dig>  the average number of synonymous changes are smaller inside squares than outside squares.
overall, the distribution of synonymous changes among edges inside squares and outside squares does not suggest that all incidences of excessive squares are due to positive selection, but it suggests that positive selection may have contributed to this excess for at least some genes.fig.
<dig> ratio of the average number of synonymous changes per edge for edges inside squares relative to edges outside squares.
the red line corresponds to a value of this ratio that is equal to one, i.e, edges inside and outside squares have the same average number of synonymous changes.
bars that end to the left  of this line indicate genes in which the average number of synonymous changes per edge is lower  inside squares than outside squares.
*, **, and *** indicate that the difference between the average number of synonymous changes inside versus outside squares is significant at p-values of  <dig> ,  <dig> , and  <dig> , respectively .
the p-values are corrected following  <cit>
using a test based on the hypergeometric distribution  <cit> , we found no significant overlap between the genes that showed evidence of positive selection in the xp-clr test and those genes among our  <dig> focal genes that  have significantly fewer synonymous mutations inside the squares than outside the squares of their haplotype network  or  had been identified in several previous studies as being subject to positive selection .
balancing selection is not a likely cause of an excess of squares
in a final analysis, we also asked for evidence of balancing selection, which manifests itself as an elevated amount of heterozygosity and can in principle produce squares.
consider, for example, the square in fig.
<dig>  in which nodes represent hypothetical diploid genotypes.
next to each circle  the nucleotide residues at positions  <dig> and  <dig> are indicated, and along the edges, the specific nucleotide changes that occurred for the first of two haplotypes.
if genotype  <dig> is the most recent common ancestor of genotypes  <dig> and  <dig>  then a substitution at site  <dig> in the first haplotype of genotype  <dig> creates genotype  <dig>  and a substitution at site  <dig> of the first haplotype creates genotype  <dig>  if balancing selection is acting on both sites , individuals  <dig> and  <dig> will be favored over individual  <dig>  because they are heterozygous at one of the two sites under balancing selection.
a further substitution to genotype  <dig>  would create a double-heterozygous genotype – and a square – that is even more favored by balancing selection.fig.
<dig> balancing selection can produce cycles.
the example indicates a hypothetical diploid genotype where two nucleotide changes occur at position  <dig> and  <dig>  circles  correspond to genotypes.
an edge connects two nodes if they differ by a single mutation.
lettering next to each node indicates the nucleotides at which two genotypes differ.
edge labels show changes required to create a genotype from its neighbor, e.g., “a20g” indicates a change from a to g at position  <dig> of the first haplotype of the hypothetical genotype.
see text for details
we computed for each gene the fraction of heterozygous individuals averaged over all sites that experienced nonsynonymous changes in at least one individual of the sample population .
among our  <dig> genes with an excess of squares, we found no significant  correlation between the number of squares and heterozygosity.
for all  <dig>  genes, we found a very small  yet significant correlation  between heterozygosity and the number of squares in a gene’s haplotype network .
in sum, balancing selection is not a likely explanation for the prevalence of squares in some genes.
multiple genes whose haplotype networks show an excess of squares are implicated in immune functions
especially prominent among the  <dig> genes whose haplotype networks show an excess of squares are genes with immune functions.
such genes are also known to be subject to frequent positive selection  <cit> .
for example, five of the  <dig> genes belong to the human leukocyte antigen  family.
these are the genes hla-a, hla-b, hla-c, hla-dpb <dig>  and hla-drb <dig>  hla genes show the highest level of polymorphisms in the human genome  <cit> , and display multiple signatures of positive selection, including a high dn / ds in antigen-recognition sites, trans-species polymorphisms, high levels of heterozygosity, as well as long range haplotypes, a key signature of recent positive selection  <cit> .
five more among the  <dig> genes with an excess of squares encode mucins, which are important for the immune response, because they help form mucus that can prevent pathogen entry, and cooperate with antibodies to fight pathogens .
these are muc <dig>  muc5b, muc <dig>  muc <dig> and muc <dig>
two more among the  <dig> genes, lilrb <dig> and lilra <dig>  encode leukocyte immunoglobulin-like receptors, which cooperate with mhc proteins.
lilrb <dig>  another member of this family, has co-evolved with hla, which is under positive selection in sub-saharan population  <cit> .
another immune-relevant gene among the  <dig> genes is fcgbp, which may play an important role in immune protection and inflammation in the intestines of primates  <cit> .
discussion
we show that the haplotype networks of  <dig> genes display a significant excess of squares that cannot be explained by chance homoplasy, genetic recombination, or balancing selection alone.
this leaves constrained evolution as a prominent candidate cause, which limits the diversity of alleles that are generated or preserved in a sequence.
while such constrained evolution can have multiple causes  <cit> , strong purifying and positive selection are most relevant for the kind of data we analyze.
strong purifying selection may play a role in the occurrence of squares, because we observed significantly fewer squares for many genes in our randomization tests when we allowed the whole protein coding sequence to change, and when we permitted substitutions to any nucleotide.
in addition, some of the genes with an excess of squares may have experienced positive selection.
first, up to 80 % of edges in the giant component of some of these genes do not have any synonymous mutations at all .
second, six of the genes with an excess of squares  have significantly more synonymous changes outside their squares than inside them .
third, six genes contained at least two adjacent windows with a significantly high value of the xp-clr test statistic that can indicate positive selection .
fourth, previous studies have suggested that  <dig> of the  <dig> genes with an excess of squares have been subject to positive selection .
finally, multiple genes among those with an excess of squares are involved in immune functions, which are frequently subject to positive selection  <cit> .
more generally, it is relevant that there is a mounting number of known genes where convergence at the sequence level has been caused by positive selection.
for example, sequence convergence occurred in the peptide-binding regions of human and mouse class ib genes in the major histocompatibility complex , the same gene family in which five members show an excess of squares in our study  <cit> .
the motor protein prestin which is involved in the mammalian auditory system has experienced adaptive sequence convergence between echolocating bats and echolocating dolphins  <cit> .
two other genes involved in the mammalian auditory system, tmc <dig> and pjvk, also have experienced convergence due to positive selection  <cit> .
in addition, whole genome sequencing of four bat species showed extensive genome-wide convergence among these taxa  <cit> .
moreover, extensive convergent evolution occurred between snake and agamid lizard mitochondrial genomes, much of which may be adaptive  <cit> .
our analysis is based on some  <dig>  human genomes, which raises the question how its results might be affected as the size of the available data set increases.
most importantly, a larger data set would lead to larger and more connected networks.
our analysis is focused on the largest connected component of each network, and increasing the size of the largest connected component could lead to more cycles just by chance alone.
indeed, larger connected components of a haplotype network in our data set also contain more cycles .
this pattern also extends to those networks with a significant excess of cycles.
specifically, giant component sizes are significantly larger for networks that have a significant excess of cycles than for the remainder of the haplotype networks .
conversely, a higher fraction of genes with an excess of cycles have large giant components .
these observations suggest that increasing the size of our data set might not just increase the overall number of cycles, but also the number of haplotype networks with an excess of cycles.
in other words, it would increase the sensitivity of our analysis.
a recent study  <cit>  has shown that hla genes show reference allele bias in the  <dig>  genomes data.
removing these alleles from the dataset could in principle lead to smaller giant components in the hla networks and hence to fewer cycles.
however, this is unlikely to materially affect our observations, because the largest components, with one exception, comprise a small fraction of the hla networks .
thus, most removed alleles would fall into other components, and their removal would thus not affect our giant-component-based analysis.
in sum, while we have not been able to explain the abundance of squares conclusively, we suggest that a mix of constrained evolution through purifying selection and positive selection may be responsible.
as data from more and more individuals from the global human population become available, it will be possible to disentangle these causes.
such data may also help explain the great differences in haplotype network structure among the human genes we characterized here.
CONCLUSIONS
we explored a novel way of representing human genetic variation data through a network-based approach whose strengths are complementary to phylogenetic trees.
despite the fact that the genes in the genomes we analyze have a shared phylogenetic history, they show very diverse properties in their haplotype networks.
specifically, these networks show different numbers of genotypes , different extents of fragmentation , different degree distributions , and different assortativity .
our analysis focuses on the feature of these networks that cannot be easily represented in phylogenetic trees, i.e., cycles.
phylogenetic trees are acyclic, and thus not ideally suited to represent evolutionary histories more complex than direct descent, such as allopolyploidization, convergent evolution, sexual reproduction, recombination and horizontal gene transfer.
such events can transform a tree-like evolutionary history into a reticulate network.
haplotype networks can represent such reticulation, and can thus complement phylogenetic trees in their ability to represent and describe evolutionary processes.

 the patterns and causes of genotypic variation in human genes have been a focus of great recent interest in evolutionary biology. different processes such as natural selection, genetic recombination, genetic drift, demography, as well as physicochemical properties of cells, can influence this diversity. various methods have been devised to represent and quantify genetic variation and to detect its causes. here we use a novel approach based on genotype networks to represent and analyze genetic variation in human genes. genotype networks are graphs that consist of nodes, which correspond to genotypes with the same phenotype, where sameness can be defined as narrowly as enzyme activity, or as broadly as viability. nodes that differ in some minimal way from each other are adjacent, i.e., connected by an edge in such a graph. the genotypes we consider are haploid genotypes of human genes in a sample of the human population, and we call two genotypes adjacent if they differ in a single nucleotide. genotype networks can be useful to address various evolutionary questions, such as how novel adaptations originate, and what role phenotypic robustness or plasticity play in adaptation. in the past, they have been mostly built from computational models of genotype-phenotype maps, but high-throughput genotyping allows us to build genotype networks from experimental data. representing such data in the form of a network makes the large analytical toolbox of graph theory available, which has been useful in fields as different as ecology and the social sciences. a common form of representation for genetic variation data is the phylogenetic tree, which shows the evolutionary relationship among a set of taxa, individuals, or genes that constitute the leaves of the tree. the common ancestors of these taxa form the interior nodes of such a tree. in a gene tree, these ancestors can be reconstructed with the help of probabilistic models of sequence evolution. phylogenetic trees are by definition acyclic graphs: they do not contain cycles – paths of edges that start from a node, pass through other nodes, and return to the same node. the acyclic nature of phylogenetic trees implies one major limitation of such trees: they cannot easily accommodate evolutionary genealogies more complex than simple vertical descent with modification. such genealogies can lead to reticulate networks of phylogenetic relationships. thus, multiple mechanisms to create genetic diversity, such as hybridization, allopolyploidization, sexual reproduction, recombination, gene conversion, and homoplasy, which lead to mosaic patterns of relationships among nodes are not easily accommodated in tree-like structures. genotype networks provide information complementary to phylogenetic trees that are not subject to this limitation, because they can accommodate cycles. figure  shows a short cycle in a hypothetical genotype network involving four dna sequences. edges reflect adjacent genotypes that differ in a single nucleotide. assume, for example, that genotype is ancestral to the other genotypes, and different substitutions produce genotypes and from it. genotype then experiences an additional a10t substitution that creates genotype this mutational path leads to a closed cycle, where three of the four edges reflect a substitution event. the fourth edge is a consequence of the first three events, because they render genotype adjacent to genotype similar scenarios can be developed if a genotype different from genotype is ancestral. regardless of this ancestor, cycles require sequence changes that render the descendants of one genotypes more similar rather than less similar. in other words, cycles require some form of homoplasy, i.e., parallel or convergent evolution. more generally, homoplasy is said to exist when two lineages display the same genetic or phenotypic characters, even though this similarity has not arisen through common ancestry.fig. a hypothetical example of a four-node cycle in a haplotype network. the example indicates a hypothetical dna sequence where two nucleotide changes occur at position and circles correspond to genotypes. an edge connects two nodes if they differ by a single mutation. lettering next to each node indicates the nucleotides at which two genotypes differ. edge labels show changes required to create a genotype from its neighbor, e.g., “a20g” indicates a change from a to g at position of the hypothetical sequence. in the example, mutations at positions and create genotypes and 3 from genotype. then, either genotype mutates at position from a to g, or genotype mutates at position from a to t, or both of these mutations take place, and create genotype 4 homoplastic sequence evolution has been documented in a wide variety of molecules. it can be caused by chance alone, which is expected to be rare in long evolving biopolymers with multiple kinds of monomers, because random mutations are more likely to cause such polymers to diverge than to converge. mutational biases, strong selective constraints on sequence evolution, positive selection, or genetic recombination can also cause homoplasy. here we construct haploid genotype networks for each of genes in the human genome, based on single nucleotide variation data available for individuals from four continents. we analyze short cycles up to length eight in these networks, and find that the haploid genotype networks of genes show a significant excess of cycles that cannot be explained by chance alone. after having excluded recombination as a prominent cause of these cycles, we focus on positive selection as a possible cause, and present evidence that in at least some of these genes positive selection may help explain the existence of cycles.  constructing and characterizing haplotype networks to construct genotype networks for human individuals, we used haploid genotypes with single nucleotide variation data available from the genomes project. thus, our genotype networks are haplotype networks, and from now on, we use the term haplotype network instead of genotype network. for each human gene, we constructed one haplotype network. two principal definitions of such a network are germane for this paper. by the first definition, a haplotype network of a human gene is a graph whose nodes correspond to protein-coding dna sequences of the gene in different individuals. two nodes are adjacent if they differ in a single base pair. by the second definition, two nodes are adjacent if their coding sequences differ by a single non-synonymous change. the second kind of network can thus also be viewed as a network of proteins or amino acid sequences, in which neighboring proteins differ in a single amino acid. we first created both dna- and protein-based haplotype networks based on the above definition, collapsing those nodes with identical sequences into one. networks can be made of one or more components. each component is a subgraph in which any two nodes are connected to each other by a path of edges. we found that the average size of the largest connected component – commonly referred to as the giant component – relative to total network size is significantly larger in protein-based networks than in dna networks. because our statistical analyses focus on the giant component of each haplotype network and work best if this component comprises as many nodes as possible, we focus on protein-based haplotype networks for the rest of this contribution. the genomes dataset we use contains information from genes, but we constructed haplotype networks only for those protein-coding sequences that showed at least one amino acid variant.fig. haplotype networks vary greatly in structure among genes. a haplotype network of the gene otog. among all protein-based haplotype networks comprising more than sequences, otog has the network with the largest giant component where all nodes fall into this component. b haplotype network of the gene hla-b, which is the most fragmented network, with nodes in components. circles in a) and b) correspond to different genotypes, while edges connect genotypes that differ by a single point mutation. circle color and size correspond to the degree of the node, where darker and larger nodes have a higher degree. c number of components versus network size for dna-based and protein-based haplotype networks. circle size in c) corresponds to the relative size of the giant component within each haplotype network figure 2a and b illustrate with two examples that haplotype networks for different genes can differ greatly in their topology. the left network, derived from the gene otog, which encodes otogelin, comprises nodes organized into a single component, whereas the right network, from gene hla-b is highly fragmented and has nodes in different components more generally, fig. 2c shows the distributions of the number of components and the size of the largest component for all genes we considered. there are networks with only a single component, but most of these networks are small, with an average of sequences. the network with the most components is the highly fragmented hla-b network with components. hla-b is known to be under strong balancing and divergent selection, which causes great genotypic diversity. this diversity translates into high network fragmentation, i.e., a network with many components. some haplotype networks have very large giant components with up to sequences. however, in most networks, the largest component is very small, comprising a maximum of ten sequences. the network with the largest giant component where all sequences fall into that component is that of otog. cycles in haplotype networks a cycle in a network is a series of edges starting from one node and ending with the same node, while passing other nodes along the cycle only once. in haplotype networks constructed from biallelic gene variants, the simplest elementary cycle, i.e. a cycle not decomposable into smaller cycles, is a square. the reason is that cycles with an odd number of edges, e.g. triangles or pentagons, are impossible when all snps are biallelic. figure  shows a square that involves the mutation of a hypothetical dna sequence at two different sites. next to each circle the nucleotide residues at these positions are indicated, and along the edges, the specific nucleotide changes that occurred. if genotype is the most recent common ancestor of its neighbors, then its two neighbors have undergone two different mutations: specifically, genotype has experienced a change from a to t at position and genotype has a change from a to g at position to form the single genotype from its ancestors, i.e. from either genotype or either genotype has undergone a change from a to g at position or genotype has undergone a change from a to t at position so that the descendants of the two ancestral sequences and become not only more similar but identical to one another. it is not necessary for both of sequences and to mutate to form genotype but a mutation in either of them can lead to the genotype and form a cycle. regardless of whether genotype or any other genotype is the common ancestor of the others, a square like this requires convergent sequence change. in long biopolymer sequences with multiple monomers that evolve through random mutation alone, cycles should be rare, because it is unlikely that mutations become reversed to create sequences more similar to one another. however cycles can be introduced by mutation biases that allow only certain residues to change, or by selection that causes only certain changes to survive, i.e., by evolutionary constraints. another possibility is recombination, which might occur between genotypes and which can result in genotype the same mechanisms can give rise to longer cycles. figure 3a shows the distribution of the number of squares in all networks. of networks had no squares. the network with the most squares is that of the gene dnah and contains squares. the inset of fig. 3a shows the distribution of hexagons and octagons. the networks with the largest number of hexagons and octagons are those of genes map2k and hla-b, respectively. note the small numbers of hexagons and octagons compared to squares. even though we enumerated elementary cycles up to length eight – beyond that, our analyses become computationally too costly – we focus most of the following analysis on squares, because they are by far the most abundant cycles.fig. distribution of the number of cycles in all networks and in networks with an excess of squares. a distribution of the number of elementary squares, as well as elementary hexagons and octagons in all protein-based networks. b distribution of the number of squares in randomized networks derived from the giant component of the hla-b network, whose number of squares is indicated by a red arrow. c number of squares in the largest components of the haplotype networks of genes with significantly more squares than expected by chance alone, together with the mean number of squares found in randomly generated networks for each gene. shaded areas depict the maximum and minimum number of squares in the randomized networks. note the logarithmic scale on the vertical axis. unconstrained or constrained mutation cannot explain the large number of cycles in many networks because some amount of homoplasy can occur by chance alone, we wished to determine whether all squares we observed could have occurred by chance homoplasy. to this end, we created randomized haplotype networks in which the same amount of evolutionary change occurred as in the actual networks. in our first randomization procedure, we created a dna sequence of the same length as the coding sequence of a gene, and created a haplotype network from it by simulating a pattern of mutations designed to yield a network with the same number of edges and the same distribution of degrees as the actual network. specifically, we compared the number of squares in each haplotype network to such randomly generated networks, and found genes whose actual number of cycles was greater than all of the randomly generated networks. thus, based on this criterion there are genes whose total number of squares cannot be explained by chance homoplasy alone at ). one can argue that this procedure does not take into consideration the actual patterns of variation observed in the data, namely that only a small subsets of sites in any one gene have been subject to mutation, and that all of the sites are biallelic, that is, only two variant nucleotides occur in them. both patterns arise from the fact that the human population sample is not highly diverged, and that natural selection constrained the evolution of these sequences, i.e., it eliminated some mutations that occurred in them. we thus modified our randomization procedure to reflect these facts. with these more conservative criteria, we still found genes whose haplotype networks have significantly more squares in their networks than expected by chance alone. that is, their number of squares cannot be explained by mutational patterns and purifying selection alone. figure 3b shows, as an example, the number of squares in the haplotype network of hla-b, which is standard deviations greater than the mean number of cycles in randomized networks. figure 3c shows the number of squares in all networks, together with the mean, minimum, and maximum number of squares for randomized networks created for each of the haplotype networks.table genes with an excess of squares in their giant component ttn mki67 obscn pkd1l2 muc16 muc17 igfn1 gpr98 prune2 syne2 ahnak2 hla-dpb1 alpk2 hla-c flg pramef2 hrnr muc5b pclo hla-a muc12 lama5 cyp2a7 hla-b poted neb muc4 pkd1l1 fbn3 dchs2 fras1 plin4 eys fcgbp tg ush2a lilrb3 lilra6 dnah17 hla-drb1 dnah5 dnah11 the number of squares in these genes cannot be explained by random homoplasy or mutational constraints. the middle column cites studies that provide evidence for positive selection, wherever such evidence is available. after fdr correction, the p-value of the statistical test comparing the actual number of cycles against that in randomized networks is for all genes additional file 4: figure s shows the distribution of elementary cycles with length four, six and eight among the genes with an excess of squares, and fig.  shows the proportion of the sequences that form part of a square in the largest connected component of each gene network. for some genes, such as poted all sequences form part of a square, and even for genes where the proportion of sequences in a square is low, such as hla-c and ttn, it exceeds 40 %.fig. proportion of sequences that are part of a cycle. proportion and actual number of sequences that are part of a square in the giant connected components of haplotype networks for those genes with a significant excess of squares we note that the networks with an excess of squares are otherwise very heterogeneous in their properties. they range from the network of mki where only nodes lie in the largest connected component, to the network of dnah, where nodes do. some of the networks have very few components, such as that of poted with a single component, whereas others have many components, such as the highly fragmented hla-b network with components. even properties within the largest connected components are heterogeneous. for example, in some networks, such as that of pkd1l the distribution of the numbers of neighbors of each sequence is highly right-skewed and dominated by sequences with few neighbors, while in others it is more symmetric. assortativity coefficients, which quantify the tendency of each node to connect to other nodes with a similar number of neighbors, also vary broadly. some networks are assortative, whereas others are disassortative. gene ontology enrichment analysis on biological processes shows several immune system-related processes that are enriched in the genes, namely “antigen processing and presentation of endogenous peptide antigens” and “interferon-gamma-mediated signaling pathway”. go enrichment analysis of molecular functions reveals the two enriched functions “calcium ion binding” and “peptide antigen binding”. “peptide antigen binding” is again associated with the immune system. given the strong representation of hla genes among genes with an excess of cycles, we also asked how the go enrichment analysis would change if we excluded the hla genes. we found a single enriched biological process, namely “o-glycan processing”, and two enriched molecular functions, namely “calcium ion binding” and “extracellular matrix constituent, lubricant activity”. we also asked whether genes with an excess of squares preferentially occurred in specific kegg or reactome pathways. six genes were preferentially associated with kegg pathways. they include tg and the genes in the hla family. the enriched pathways comprise “epstein-barr virus infection”, “autoimmune thyroid disease”, “htlv-i infection”, “viral myocarditis”, “allograft rejection”, “phagosome”, “antigen processing and presentation”, “graft-versus-host disease”, “cell adhesion molecules ”, “herpes simplex infection”, and “type i diabetes mellitus”. for reactome pathways, we found twelve genes enriched in six pathways. the genes include those encoding mucins, the hla family and lilr family genes. the enriched pathways are “termination of o-glycan biosynthesis”, “interferon gamma signaling”, “endosomal/vacuolar pathway”, “immunoregulatory interactions between a lymphoid and a non-lymphoid cell”, “antigen presentation: folding, assembly and peptide loading of class i mhc” and “defective galnt causes colorectal cancer ”. we note that both enriched kegg and reactome pathways include several immunity-related pathways. recombination cannot account for an excess of squares in most networks to exclude the possibility that genetic recombination may account for the excess of squares in some networks, we performed two complementary analyses. first, we simulated for each gene the effect of recombination on haplotype network structure by creating haplotype networks based on a set of sequences that were subject to approximately as many recombination events as occurred in the human population since their common ancestry, as well as to as many mutations as there are edges in the network. we repeated this process times for each gene, creating simulated haplotype networks, and counted the number of squares in them. for each of the genes, the empirical network showed more squares than each of the simulated networks. in the second analysis, we asked whether gene conversion, a process of unidirectional recombination in which only one of the recombining sequences changes, may have caused the excess of squares. to this end, we used the program geneconv to detect gene conversion among the sequences in the giant components of the haplotype networks. we used sequences comprising both synonymous and non-synonymous changes to give the program more power in finding gene conversion events. only one gene showed any sign of gene conversion, and it did so for only two of sequences in cyp2a. in sum, based on these analyses, it seems unlikely that recombination can explain the excess of squares we observe in the haplotype networks of genes. positive selection as a potential cause of squares positive selection can be a driver for homoplastic or convergent evolution, where two separate lineages evolve the same trait independently. because such adaptive homoplasy can occur not only at the phenotypic level, but also at the sequence level, we wished to find out whether positive selection can help explain the excess of squares we observed in the haplotype networks of genes. previous studies had indeed indicated positive selection for at least of the genes. in addition, we used results from a branch-site likelihood test which indicates positive selection based on a ratio dn / ds >  observed along one or more branches of a phylogenetic tree. this test has been applied to vertebrate genes in the selectome database, which indicates that of our genes with abundant squares show patterns of positive selection, either in primates or in the bony vertebrates. this number – of – is unlikely to be observed by chance alone ). in addition, we used the xp-clr test for neutrality. the test compares different populations to identify rapid changes in a locus’ allele frequency that cannot be explained by random drift alone. in applying this test, we used a test statistics pre-computed over 2 kb sliding windows that covered the human genome, and asked for each of our genes whose haplotype network showed an excess of squares, whether two or more of the windows where the test-statistic indicated the action of positive section overlapped with the gene. by this criterion, six of our genes showed evidence of positive selection in at least one population. overall, of our genes with an excess of squares showed signs of positive selection by at least one of these criteria or by previous work. we also analyzed patterns of synonymous and non-synonymous changes in more detail. a commonly used indicator of positive selection for two protein-coding dna sequences is dn / ds, i.e. the ratio of nonsynonymous changes dn per nonsynonymous site to synonymous changes ds per synonymous site. values of dn / ds >  can indicate positive selection. unfortunately, dn / ds can be computed only for sequences more distantly related than those in our haplotype networks. the reason is that in these networks, adjacent sequence pairs differ only in a single nonsynonymous mutation, and many adjacent pairs do not even show a single synonymous change. more specifically, in the giant component of our networks, up to 80 % of sequence pairs do not show a single synonymous mutation, and this incidence of synonymous mutations is similarly low in the entire network. moreover, it has been suggested that for very closely related sequences, dn / ds is not a sensitive indicator of positive selection. for these reasons, we compared the incidence of nonsynonymous and synonymous changes among groups of edges, reasoning that groups of edges with very few synonymous changes might provide hints that some or all members of the group may have been subject to positive selection. most edges show no synonymous changes at all in some networks, which hints that positive selection may have played a role in creating their pattern of diversity. we specifically compared edges with no synonymous change inside squares and outside squares. while the fractions of edges without synonymous changes inside squares was not significantly different from those outside squares, the average number of synonymous changes on edges inside squares was significantly smaller than that outside squares for 14 % of the genes. figure  shows the average number of synonymous changes per edge for edges inside squares divided by that for edges outside squares. for genes where this ratio is below the average number of synonymous changes are smaller inside squares than outside squares. overall, the distribution of synonymous changes among edges inside squares and outside squares does not suggest that all incidences of excessive squares are due to positive selection, but it suggests that positive selection may have contributed to this excess for at least some genes.fig. ratio of the average number of synonymous changes per edge for edges inside squares relative to edges outside squares. the red line corresponds to a value of this ratio that is equal to one, i.e, edges inside and outside squares have the same average number of synonymous changes. bars that end to the left of this line indicate genes in which the average number of synonymous changes per edge is lower inside squares than outside squares. *, **, and *** indicate that the difference between the average number of synonymous changes inside versus outside squares is significant at p-values of,, and, respectively. the p-values are corrected following using a test based on the hypergeometric distribution, we found no significant overlap between the genes that showed evidence of positive selection in the xp-clr test and those genes among our focal genes that have significantly fewer synonymous mutations inside the squares than outside the squares of their haplotype network or had been identified in several previous studies as being subject to positive selection. balancing selection is not a likely cause of an excess of squares in a final analysis, we also asked for evidence of balancing selection, which manifests itself as an elevated amount of heterozygosity and can in principle produce squares. consider, for example, the square in fig.  in which nodes represent hypothetical diploid genotypes. next to each circle the nucleotide residues at positions and are indicated, and along the edges, the specific nucleotide changes that occurred for the first of two haplotypes. if genotype is the most recent common ancestor of genotypes and then a substitution at site in the first haplotype of genotype creates genotype and a substitution at site of the first haplotype creates genotype if balancing selection is acting on both sites, individuals and will be favored over individual because they are heterozygous at one of the two sites under balancing selection. a further substitution to genotype would create a double-heterozygous genotype – and a square – that is even more favored by balancing selection.fig. balancing selection can produce cycles. the example indicates a hypothetical diploid genotype where two nucleotide changes occur at position and circles correspond to genotypes. an edge connects two nodes if they differ by a single mutation. lettering next to each node indicates the nucleotides at which two genotypes differ. edge labels show changes required to create a genotype from its neighbor, e.g., “a20g” indicates a change from a to g at position of the first haplotype of the hypothetical genotype. see text for details we computed for each gene the fraction of heterozygous individuals averaged over all sites that experienced nonsynonymous changes in at least one individual of the sample population. among our genes with an excess of squares, we found no significant correlation between the number of squares and heterozygosity. for all genes, we found a very small yet significant correlation between heterozygosity and the number of squares in a gene’s haplotype network. in sum, balancing selection is not a likely explanation for the prevalence of squares in some genes. multiple genes whose haplotype networks show an excess of squares are implicated in immune functions especially prominent among the genes whose haplotype networks show an excess of squares are genes with immune functions. such genes are also known to be subject to frequent positive selection. for example, five of the genes belong to the human leukocyte antigen family. these are the genes hla-a, hla-b, hla-c, hla-dpb and hla-drb hla genes show the highest level of polymorphisms in the human genome, and display multiple signatures of positive selection, including a high dn / ds in antigen-recognition sites, trans-species polymorphisms, high levels of heterozygosity, as well as long range haplotypes, a key signature of recent positive selection. five more among the genes with an excess of squares encode mucins, which are important for the immune response, because they help form mucus that can prevent pathogen entry, and cooperate with antibodies to fight pathogens. these are muc muc5b, muc muc and muc two more among the genes, lilrb and lilra encode leukocyte immunoglobulin-like receptors, which cooperate with mhc proteins. lilrb another member of this family, has co-evolved with hla, which is under positive selection in sub-saharan population. another immune-relevant gene among the genes is fcgbp, which may play an important role in immune protection and inflammation in the intestines of primates. discussion we show that the haplotype networks of genes display a significant excess of squares that cannot be explained by chance homoplasy, genetic recombination, or balancing selection alone. this leaves constrained evolution as a prominent candidate cause, which limits the diversity of alleles that are generated or preserved in a sequence. while such constrained evolution can have multiple causes, strong purifying and positive selection are most relevant for the kind of data we analyze. strong purifying selection may play a role in the occurrence of squares, because we observed significantly fewer squares for many genes in our randomization tests when we allowed the whole protein coding sequence to change, and when we permitted substitutions to any nucleotide. in addition, some of the genes with an excess of squares may have experienced positive selection. first, up to 80 % of edges in the giant component of some of these genes do not have any synonymous mutations at all. second, six of the genes with an excess of squares have significantly more synonymous changes outside their squares than inside them. third, six genes contained at least two adjacent windows with a significantly high value of the xp-clr test statistic that can indicate positive selection. fourth, previous studies have suggested that of the genes with an excess of squares have been subject to positive selection. finally, multiple genes among those with an excess of squares are involved in immune functions, which are frequently subject to positive selection. more generally, it is relevant that there is a mounting number of known genes where convergence at the sequence level has been caused by positive selection. for example, sequence convergence occurred in the peptide-binding regions of human and mouse class ib genes in the major histocompatibility complex, the same gene family in which five members show an excess of squares in our study. the motor protein prestin which is involved in the mammalian auditory system has experienced adaptive sequence convergence between echolocating bats and echolocating dolphins. two other genes involved in the mammalian auditory system, tmc and pjvk, also have experienced convergence due to positive selection. in addition, whole genome sequencing of four bat species showed extensive genome-wide convergence among these taxa. moreover, extensive convergent evolution occurred between snake and agamid lizard mitochondrial genomes, much of which may be adaptive. our analysis is based on some human genomes, which raises the question how its results might be affected as the size of the available data set increases. most importantly, a larger data set would lead to larger and more connected networks. our analysis is focused on the largest connected component of each network, and increasing the size of the largest connected component could lead to more cycles just by chance alone. indeed, larger connected components of a haplotype network in our data set also contain more cycles. this pattern also extends to those networks with a significant excess of cycles. specifically, giant component sizes are significantly larger for networks that have a significant excess of cycles than for the remainder of the haplotype networks. conversely, a higher fraction of genes with an excess of cycles have large giant components. these observations suggest that increasing the size of our data set might not just increase the overall number of cycles, but also the number of haplotype networks with an excess of cycles. in other words, it would increase the sensitivity of our analysis. a recent study has shown that hla genes show reference allele bias in the genomes data. removing these alleles from the dataset could in principle lead to smaller giant components in the hla networks and hence to fewer cycles. however, this is unlikely to materially affect our observations, because the largest components, with one exception, comprise a small fraction of the hla networks. thus, most removed alleles would fall into other components, and their removal would thus not affect our giant-component-based analysis. in sum, while we have not been able to explain the abundance of squares conclusively, we suggest that a mix of constrained evolution through purifying selection and positive selection may be responsible. as data from more and more individuals from the global human population become available, it will be possible to disentangle these causes. such data may also help explain the great differences in haplotype network structure among the human genes we characterized here.  we explored a novel way of representing human genetic variation data through a network-based approach whose strengths are complementary to phylogenetic trees. despite the fact that the genes in the genomes we analyze have a shared phylogenetic history, they show very diverse properties in their haplotype networks. specifically, these networks show different numbers of genotypes, different extents of fragmentation, different degree distributions, and different assortativity. our analysis focuses on the feature of these networks that cannot be easily represented in phylogenetic trees, i.e., cycles. phylogenetic trees are acyclic, and thus not ideally suited to represent evolutionary histories more complex than direct descent, such as allopolyploidization, convergent evolution, sexual reproduction, recombination and horizontal gene transfer. such events can transform a tree-like evolutionary history into a reticulate network. haplotype networks can represent such reticulation, and can thus complement phylogenetic trees in their ability to represent and describe evolutionary processes. 
we focus on those genes whose genotype networks show many cycles, which can indicate homoplasy, i.e., parallel or convergent evolution, on the sequence level. the structure of these networks varies widely among genes, indicating different patterns of variation despite a shared evolutionary history. when analyzing possible explanations, we discovered evidence for positive selection in of these genes and, in addition, a potential role for constrained variation and purifying selection. genotype networks are representations of genetic variation data that are complementary to phylogenetic trees. genotype networks are representations of genetic variation data that can help understand unusual patterns of genomic variation. genotype networksgenetic variationnatural selectionhuman genomehttp://dx.doi.org/ /501100001711schweizerischer nationalfonds zur förderung der wissenschaftlichen forschung31003a_146137wagner andreas issue-copyright-statement© the author 2016 for genes, the observed number of cycles is so large that it cannot be explained by either chance homoplasy or recombination. a genotype network is a graph whose nodes are genotypes with the same broadly defined phenotype. two nodes are connected if they differ in some minimal way, e.g., in a single nucleotide.
 genotype networks are representations of genetic variation data that are complementary to phylogenetic trees. a genotype network is a graph whose nodes are genotypes with the same broadly defined phenotype. two nodes are connected if they differ in some minimal way, e.g., in a single nucleotide.  we analyze human genome variation data from the genomes project, and construct haploid genotype networks for protein coding genes. the structure of these networks varies widely among genes, indicating different patterns of variation despite a shared evolutionary history. we focus on those genes whose genotype networks show many cycles, which can indicate homoplasy, i.e., parallel or convergent evolution, on the sequence level.  for genes, the observed number of cycles is so large that it cannot be explained by either chance homoplasy or recombination. when analyzing possible explanations, we discovered evidence for positive selection in of these genes and, in addition, a potential role for constrained variation and purifying selection. balancing selection plays at most a small role. the genes with excess cycles are enriched in functions related to immunity and response to pathogens. genotype networks are representations of genetic variation data that can help understand unusual patterns of genomic variation. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords genotype networksgenetic variationnatural selectionhuman genomehttp://dx.doi.org/ /501100001711schweizerischer nationalfonds zur förderung der wissenschaftlichen forschung31003a_146137wagner andreas issue-copyright-statement© the author 2016 