Complex fitness landscape shapes variation in a hyperpolymorphic species

It is natural to assume that patterns of genetic variation in hyperpolymorphic species can reveal large-scale properties of the fitness landscape that are hard to detect by studying species with ordinary levels of genetic variation. Here, we study such patterns in a fungus Schizophyllum commune, the most polymorphic species known. Throughout the genome, short-range linkage disequilibrium (LD) caused by attraction of minor alleles is higher between pairs of nonsynonymous than of synonymous variants. This effect is especially pronounced for pairs of sites that are located within the same gene, especially if a large fraction of the gene is covered by haploblocks, genome segments where the gene pool consists of two highly divergent haplotypes, which is a signature of balancing selection. Haploblocks are usually shorter than 1000 nucleotides, and collectively cover about 10% of the S. commune genome. LD tends to be substantially higher for pairs of nonsynonymous variants encoding amino acids that interact within the protein. There is a substantial correlation between LDs at the same pairs of nonsynonymous mutations in the USA and the Russian populations. These patterns indicate that selection in S. commune involves positive epistasis due to compensatory interactions between nonsynonymous alleles. When less polymorphic species are studied, analogous patterns can be detected only through interspecific comparisons.

Perhaps, the fitness landscape is complex macroscopically but is more smooth microscopically or, in other words, epistasis is genuinely more pronounced at a macroscopic scale (Ochs and Desai, 2015). If so, studying epistasis in hyperpolymorphic populations, where differences between genotypes can be as high as those between genomes of species from different genera or even families, holds a great promise because variation within such a population can cover multiple fitness peaks or a sizeable chunk of a curved ridge of high fitness (Bateson, 1909;Dobzhansky, 1937;Gavrilets, 1997;Kondrashov et al., 2002;Muller, 1942; Appendix 1). eLife digest Changes to DNA known as mutations may alter how the proteins and other components of a cell work, and thus play an important role in allowing living things to evolve new traits and abilities over many generations. Whether a mutation is beneficial or harmful may differ depending on the genetic background of the individual -that is, depending on other mutations present in other positions within the same gene -due to a phenomenon called epistasis.
Epistasis is known to affect how various species accumulate differences in their DNA compared to each other over time. For example, a mutation that is rare in humans and known to cause disease may be widespread in other primates because its negative effect is canceled out by another mutation that is standard for these species but absent in humans. However, it remains unclear whether epistasis plays a significant part in shaping genetic differences between individuals of the same species.
A type of fungus known as Schizophyllum commune lives on rotting wood and is found across the world. It is one of the most genetically diverse species currently known, so there is a higher chance of pairs of compensatory mutations occurring and persisting for a long time in S. commune than in most other species, providing a unique opportunity to study epistasis.
Here, Stolyarova et al. studied two distinct populations of S. commune, one from the USA and one from Russia. The team found that -unlike in humans, flies and other less genetically diverse species -epistasis maintains combinations of mutations in S. commune that individually would be harmful to the fungus but together compensate for each other. For example, pairs of mutations affecting specific molecules known as amino acids -the building blocks of proteins -that physically interact with each other tended to be found together in the same individuals.
One potential downside of having pairs of compensatory mutations in the genome is that when the organism reproduces, the process of making sex cells may split up these pairs so that harmful mutations are inherited without their partner mutations. Thus, epistasis may have helped shape the way S. commune and other genetically diverse species have evolved.

Elevated LD between nonsynonymous polymorphisms
In a vast majority of species, nucleotide diversity π, the evolutionary distance between a pair of randomly chosen genotypes, is, at selectively neutral sites, of the order of 0.001 (as in Homo sapiens) or 0.01 (as in Drosophila melanogaster) (Leffler et al., 2012). Still, a few hyperpolymorphic species with π>0.1 are known, of which the wood-decaying fungus Schizophyllum commune is the most extreme, where π=0.20 or 0.13 in the USA or the Russian populations, respectively (Baranova et al., 2015;Appendix 3- figure 1). The two populations of S. commune are highly divergent (dS between populations ≈ 0.34, F st = 0.58), but there is essentially no structure within either of them (Appendix 3-figure 2). We studied 34 haploid genotypes from the USA and 21 from Russia, each obtained by sequencing and de novo assembly of a haploid culture originated from a single haplospore. The use of haploid samples and de novo assembly of each sample ensures robust identification of haplotypes. We then compared the LD between nonsynonymous SNPs (LD nonsyn ) to that between synonymous SNPs (LD syn ).
In both S. commune populations, at sites with minor allele frequency (MAF) >0.05, LD nonsyn is much higher than LD syn at the same nucleotide distance ( Figure 1A, Figure 1-figure supplement 1A). This excess of LD nonsyn is much stronger for pairs of SNPs located within the same gene, compared to pairs of SNPs from adjacent genes at the same distance. By contrast, the excess of LD nonsyn is independent of whether the two SNPs are located within the same or in different exons of a gene (Figure 1-figure supplement 2). In S. commune, the recombination rate is higher within exons (Seplyarskiy et al., 2014), which may affect the patterns of LD; however, this factor could only reduce within-gene LD, and in any case cannot explain the difference between LD nonsyn and LD syn . For S. commune, the excess of LD nonsyn over LDsyn holds when we explicitely control for MAFs ( Figure 1-figure supplement 3), indicating that differences in MAFs between synonymous and nonsynonymous polymorphisms cannot explain the excess of LD nonsyn .
A much weaker excess of LD nonsyn over LD syn for MAF >0.05 is also observed in the less genetically diverse D. melanogaster population ( Figure 1B). In the still less polymorphic human populations, LD nonsyn is indistinguishable from LD syn at the same distances ( Figure 1C, Figure 1 The excess of LD nonsyn over LD syn corresponds to the attraction between minor nonsynonymous alleles. This attraction can only appear due to positive epistasis between such alleles -higher-thanexpected fitness of their combinations (Appendix 2). Positive epistasis can be expected to cause stronger LD in more polymorphic populations ( Figure 1D-F, Appendix 1) and must be more common for pairs of sites located within the same gene, which are more likely to interact with each other.
For rare SNPs with MAF <0.05 taken alone, LD nonsyn is similar or lower to LD syn for all three species, consistent with the effects of random drift, Hill-Robertson interference, and/or negative epistasis (Figure 1-figure supplements 4-6; Appendix 2). Decreased LD between negatively selected polymorphisms is expected due to Hill-Robertson interference between deleterious alleles (Hill and Robertson, 1966;Roze and Barton, 2006); this effect has been described previously for H. sapiens (Garcia and Lohmueller, 2021) and D. melanogaster (Sandler et al., 2021) and is observed in our simulations (Appendix 2- figure 4). In addition, LD nonsyn can be reduced by negative epistasis between deleterious alleles (Garcia and Lohmueller, 2021), similarly to the negative LD detected among lossof-function polymorphisms in humans, flies and plants (Sandler et al., 2021;Sohail et al., 2017).

Elevated LD between interacting sites
Natural selection acting on physically interacting amino acids that are located close to each other within the three-dimensional structure of a protein is characterized by strong epistasis which leads to their coevolution at the level of between-species differences (Marks et al., 2011;Ovchinnikov et al., 2014;Sjodt et al., 2018). Genome-wide elevated LD between amino acid sites within structural domains was recently demonstrated in human populations (Ragsdale, 2021). Extraordinary diversity of S. commune makes it possible to observe an analogous phenomenon at the level of individual genes in within-population variation.
To test this, we aligned S. commune proteins to the PDB database of protein structures. In the obtained set of 5188 genes with a good match to a protein with known structure, we identified pairs of codons in the S. commune genome encoding amino acid residues positioned near each other (within 10 Å) in the protein structures, and calculated the average LD between SNPs in such pairs of codons. Naturally, pairs of physically interacting sites are more likely to be closely spaced in the gene sequence and therefore to be under a higher LD than non-interacting ones. To account for this, we discarded pairs of SNPs within five amino acids from each other, and used a controlled permutation test (see Materials and methods) to compare the LD between physically close pairs of sites to that between distant pairs of sites.
In both S. commune populations, pairs of nonsynonymous SNPs are in stronger LD when they are located at codons encoding physically close than distant amino acids (      p-value <1e-3). This is not the case for pairs of synonymous SNPs (Figure 2A; permutation test p-value = 0.58). Moreover, it is possible to identify individual proteins with significant associations between the patterns of LD and of physical interactions between sites. At a 5% FDR, we found 22 such proteins in the USA population, and 87 proteins in the Russian population (Appendix 3-table 1); three examples are shown in Figure 2B-D (see also Figure 2- figure supplements 1 and 2). The alignment of ADAT2 protein contains two segments (teal and red in Figure 2B) characterized by high withinsegment LD. The boundaries of these segments match those of structural units of the protein, but not the exon structure of its gene. In RadB protein, a similar pattern is observed, and LD is also elevated between pairs of SNPs from different segments on the interface of the corresponding structural units ( Figure 2C). The alignment of 4CL protein can be naturally split into four high-LD segments, which also match its structure ( Figure 2D).

Distinct regions of high LD
The magnitude of LD varies widely along the S. commune genome. Visual inspection of the data shows a salient pattern of regions of relatively low LD, alternating with mostly short regions of high LD (haploblocks, Figure 3-figure supplement 1). We calculated LD along the genome in a sliding window of 250 nucleotides and regarded as a haploblock any continuous genomic region with LD values that belong to the heavy tail of its distribution (see Materials and methods).
In the USA population, 8.4% of the genome is occupied by 5,316 such haploblocks, 56% consist of regions with background LD level, and the rest cannot be analyzed due to poor alignment quality or low SNP density. Eighty-eight percent of the haploblocks are shorter than 1000 nucleotides, although the longest haploblocks spread for several thousand nucleotides (Figure 3-figure supplement 2). In the Russian population, there are 10,694 haploblocks, occupying 15.9% of the genome, and regions of background LD cover 39% of it. There is only a modest correlation between the USA and Russian haploblocks: the probability that a genomic position belongs to a haploblock in both populations is 2.3% instead of the expected 1.3%, indicating their relatively short persistence time in the populations (examples shown in Figure 3-figure supplement 1).
LD within a haploblock is usually so high that most genotypes can be attributed to one of just two distinct haplotypes, which carry different sets of alleles ( Figure 3-figure supplement 3). This results in a bimodal distribution of the fraction of minor alleles in a genotype within a haploblock, because some genotypes belong to the major haplotype and, thus, carry only a small fraction of minor alleles, and other genotypes belong to the minor haplotype and, thus, possess a high fraction of minor alleles ( Figure 3A). Polymorphic sites within haploblocks are characterized by higher MAF than that at sites that reside in non-haploblock regions (t-test p-value <2e-16), and in the USA population MAFs within a haploblock are positively correlated with its strength of LD ( Figure 3B, Pearson correlation estimate = 0.07, p-value <2e-6).
There is no one-to-one correspondence between haploblocks and genes, which are, on average, longer. Still, different genes are covered by haploblocks to different extent, which leads to wide variation in the strength of LD and other characteristics among them. Genes with high LD, that is those that contain haploblocks, have the largest excess of LD nonsyn over LD syn ( Figure 3C). Positive correlation between the overall LD within the gene and the excess of LD nonsyn in this gene indicates that the attraction between nonsynonymous variants, driven by epistasis, is stronger if combinations of epistatic alleles are persisting within population for a long time, comprising haplotypes within a haploblock. Since both haplotypes tend to be common in a haploblock (Figure 3), this excess is much stronger for loci with MAF >0.05. three-dimensional structures. Heatmaps show the physical distance between pairs of sites in the protein structure; only positions carrying biallelic SNPs are shown. Black dots correspond to pairs of sites with high LD (>0.9 quantile for the gene). Dashed lines in (c) structure show high LD between physically close SNPs from different segments of high LD. In these examples, LD is calculated in the Russian population of S. commune.
The online version of this article includes the following figure supplement(s) for figure 2:   LD between alleles of all kinds is higher within genes with large ratios of nonsynonymous and synonymous polymorphisms p n /p s (Spearman correlation p-value <2e-16, Figure 3D). Genes with elevated p n /p s also have a stronger excess of LD nonsyn over LD syn ( Figure 3E, Spearman correlation p-value = 4.4e-17). This excess is the strongest for genes with high overall LD, but its correlation with p n /p s holds even when the overall LD is controlled for (Figure 3-figure supplement 5).
There can be multiple non-exclusive mechanisms by which epistasis could lead to the observed positive associations between pn/ps, overall LD, and excess LD nonsyn . First, genes under weaker selection, and therefore higher pn/ps, could be characterized by a higher overall amount and/or strength of epistasis. Second, epistasis, as estimated by excess LD nonsyn , can contribute to increased pn/ps by allowing nonsynonymous polymorphisms to segregate in the population when maintained in        coadapted combinations, therefore weakening negative selection against them. Third, epistasis can be more potent in genes with lower overall recombination rate due to competition between epistasis and recombination: recombination breaks positively interacting combinations of alleles, disrupting linkage between them and interfering with epistasis. Fourth, existence of cosegregating combinations of mutually beneficial alleles could select for reduced local recombination rate.

Excess of LD nonsyn requires stable polymorphism
Simulations show that positive epistasis alone cannot lead to the observed large excess LD nonsyn over LD syn , for which two extra conditions need to be satisfied. The general reason for this is simple: in order for a substantial LD between not-too-rare alleles to appear, these alleles must persist in the population for a long enough time.
First, positive epistasis must lead to a full compensation of deleterious effects of individual alleles. In other words, the fitnesses of at least two most-fit genotypes that are present in the population at substantial frequencies must be (nearly) the same (Figure 3-figure supplement 6). If this is not the case, selection favoring the only most-fit genotype leads to a too low level of genetic variation, which persists only due to recurrent mutation. It is natural to assume that the two major haplotypes that are common within a haploblock correspond to high-fitness genotypes. High-fitness genotypes can represent either isolated fitness peaks of equal heights (corresponding to a situation when two out of the four allele combinations confer high fitness) or a flat, curved ridge of high fitness (corresponding to a situation when three out of four combinations confer high fitness). The available data are insufficient to distinguish between these two options. Of course, with complete selective neutrality of all allele combinations there is no reason for LD nonsyn >LD syn , so that at least some mixed genotypes, carrying alleles from different high-fitness genotypes, must be maladapted.
Second, there must be some kind of balancing selection that specifically works to maintain variation, because otherwise random drift does not allow genetic variation to persist for a long enough time even if some, or even all, genotypes are equally fit (Figure 3-figure supplement 6). Here, there are at least two options. On the one hand, a 'real' negative frequency-dependent selection (NFDS) can act either directly at loci that display high LD or at some other tightly linked loci (Charlesworth, 2006;Olendorf et al., 2006). On the other hand, variation can be maintained due to associative overdominance (AOD), resulting from selection against recurrent deleterious mutations at linked loci (Gilbert et al., 2020;Ohta, 1971;Zhao and Charlesworth, 2016).
Balancing selection is also neccessary for the presence of haploblocks, because a pair of divergent haplotypes can evolve in a panmictic population only if they coexist for a considerable time. A single locus under NFDS is enough to maintain a haploblock comprising the region of the genome around it. By contrast, if variation is maintained by AOD, it is more likely that selection against recessive mutations acts at a number of tightly linked loci (Gilbert et al., 2020). Long coexistence of diverged haplotypes that comprise a haploblock enables accumulation of co-adapted combinations of nonsynonymous alleles within them. Thus, it is not surprising that a pronounced excess of LD nonsyn over LD syn in S. commune is observed primarily within haploblocks and that this excess is higher in genes with higher p n /p s .

Correlated LDs in two populations
Although a high excess of LD nonsyn is observed only within haploblocks, a signature of epistasis can also be seen outside of them in the form of a correlation between LDs in the two populations. This correlation can be high even if LDs by themselves are low.
The USA and the Russian populations share a large proportion of their SNPs. Given the high divergence between the two populations, few such shared SNPs are expected to have common origin in the ancestral population, and instead they are likely to have arisen from recurrent mutation. Since the haploblocks show little correlation between the two populations, we assume that they arose after their divergence. The high prevalence of coincident SNPs is not surprising because SNPs comprise 0.28 and 0.13 of all the aligned nucleotide sites in the USA and Russian populations, respectively (Baranova et al., 2015, Appendix 3- figure 2). We identified pairs of shared biallelic SNPs located within 2 kb from one another and calculated the LD between them in both populations. To avoid the effects of strong within-population linkage and the occasional co-ocсurrence of haploblocks between populations, we excluded SNPs located within haploblocks or within genes under high LD (>0.8 LD quantile for the corresponding population) in either population.
The values of LD in the two populations are strongly correlated only for pairs of nonsynonymous SNPs located within the same gene, and only if both populations carry the same pairs of amino acids in the same sites ( Figure 4). The correlation of LDs is the strongest if shared SNPs carry the same pairs of nucleotides, but is also observed if they encode the same amino acids by different nucleotides (Figure 4-figure supplement 1). The contrast between correlations within pairs of sites that reside in the same vs. different genes and the correlation of LDs observed for different nucleotides encoding the same amino acid cannot be explained by inheritance of LD from the common ancestral population. Moreover, synonymous SNPs are expected to be on average older than nonsynonymous ones, so that this mechanism should lead to a higher correlation of LDs for pairs of synonymous mutations. Thus, the observed pattern indicates that epistatic selection is shared between the two populations.
The correlation of LDs between SNPs located within haploblocks in both populations is high regardless of whether they reside in the same or different genes, apparently because of occasional coincidence of haploblocks between populations (Figure 4-figure supplement 2).

Discussion
On top of its most salient property, an exceptionally high π, genetic variation within S. commune possesses two other pervasive features. The first is a high prevalence of mostly short haploblocks, genome segments comprising two or occasionally three distinct haplotypes, which is a signature of  balancing selection. The overall fraction of the genome covered by haploblocks is ~10%, which is about an order of magnitude higher than the fraction covered by detectable signatures of balancing selection in genomes of other species (DeGiorgio et al., 2014;Leffler et al., 2013;Rasmussen et al., 2014). The second feature is the excessive attraction between nonsynonymous alleles polarized by frequency. This pattern is much stronger within haploblocks, indicating that they were shaped by both balancing and epistatic selection, so that amino acids common within a haplotype together confer a higher fitness. Polymorphisms that involve haplotypes that comprise many interacting genes, such as inversions (Charlesworth and Charlesworth, 1973;Dobzhansky and Pavlovsky, 1958;Singh, 2008;Sturtevant and Mather, 1938) and supergenes (Joron et al., 2011;Mather, 1950), are known from the dawn of population genetics, but here we are dealing with an analogous phenomenon at a much finer scale, because haploblocks are typically shorter than genes. Thus, instead of coadapted gene complexes (Dobzhansky and Pavlovsky, 1958), haplotypes represent coadaptive site complexes within genes.
In our simulations, equally high fitnesses of two or more genotypes was a necessary condition for a large excess of LD nonsyn , because otherwise the polymorphism did not live long enough for any substantial LD to evolve. However, epistasis between loci responsible for real or apparent balancing selection and those involved in compensatory interactions probably abolished the need for this finetuning of fitnesses. For example, if each haploblock carries its own complement of partially recessive deleterious mutations, together with alleles engaged in compensatory interactions with each other which also make these recessive mutations less deleterious, AOD can be expected to cause stable coexistence of these alleles.
Why are haploblocks and positive LD between minor nonsynonymous alleles so common in S. commune, but not in other, less polymorphic, species? There may be several, not mutually exclusive, reasons for this. Regarding haploblocks, real or apparent balancing selection may be more common in S. commune due to its higher polymorphism. Also, the same balancing selection may protect polymorphism in a huge population of S. commune, but not in populations with lower N e . Finally, an excess of haploblocks in S. commune may be at least due to better detection of signatures of balancing selection in a species with an extraordinary density of SNPs.
Excessive LD nonsyn in S. commune is also likely to be due to its hyperpolymorphism which increases the probability that mutually compensating alleles at a pair of interacting sites achieve high frequency and encounter each other in the same haplotype before being eliminated by selection. In other words, even if the fitness landscape remains the same, it results in more epistatic selection and, thus, in stronger LD in a species whose genetic variation covers a larger chunk of this landscape (Figure 1).
In a vast majority of species, π is a small parameter <<1. This imposes a severe constraint on operation of selection and obscures signatures of its particular modes. Thus, hyperpolymorphic species where π is ~1 provide a unique opportunity to study phenomena which are traditionally viewed as belonging to the domain of macroevolution through data on within-population variation.

S. commune sampling, sequencing, and assembly
Haploid cultures of 24 isolates, each originated from a single haplospore, were obtained from fruit bodies collected in Ann Arbor, MI, USA by T. James and A. Kondrashov table 2. To obtain isolates, wild fruit bodies were hung on the top lid of a 10 cm petri dish with agar medium. Petri dish was set at an angle of 60-70 degrees to the horizontal surface for 32 hr. A germinated spore was excised together with a square-shaped fragment (approximately 0.7 × 0.7 mm) of the medium from the maximally rarefied area of the obtained spore print under a stereomicroscope with 100 x magnification. The obtained isolates were cultured in Petri dishes on 2% malt extract agar for a week. For storage, cultures were subcultured into 1.5 ml microcentrifuge tubes with 2% malt extract agar. To obtain sufficient biomass for DNA isolation, isolates were cultured in 20 ml 0.5% malt extract liquid medium in 50 ml microcentrifuge tubes in a horizontal position on a shaker at 100 rpm in daylight for 5-10 days. The tubes with the cultures were then centrifuged at 4000 rpm, and the supernatant was decanted. The resulting mycelium was lyophilized. DNA was extracted using Diamond DNA kit according to the manufacturer's recommendations.
DNA libraries were constructed using the NEBNext Ultra II DNA Library Prep Kit kit by New England Biolabs (NEB) and the NEBNext Multiplex Oligos for Illumina (Index Primers Set 1) by NEB following the manufacturer's protocol. The samples were amplified using 10 cycles of PCR. The constructed libraries were sequenced on Illumina NextSeq500 with paired-end read length of 151. The genomes were assembled de novo using SPAdes (v3.6.0) (Bankevich et al., 2012); possible contaminations were removed using blobology (Kumar et al., 2013). Average N50 was ~165 kb for USA samples and ~70 kb for Russian samples (assembly statistics are provided in Appendix 3-table 2).
Together with the 30 samples sequenced previously (Baranova et al., 2015;Bezmenova et al., 2020), the obtained haploid genomes were aligned with TBA and multiz (Blanchette et al., 2004) and projected onto the reference scaffolds (Ohm et al., 2010). Ortholog sequences were extracted on the basis of the reference genome annotation (Ohm et al., 2010) and realigned using macse codonbased aligner (Ranwez et al., 2011). The alignments are available at https://makarich.fbb.msu.ru/ astolyarova/schizophyllum_data/. Only the gap-free columns of the whole-genome alignment and the orthologs that were found in all 55 genomes were used for analysis. The total number of detected SNPs was 5.8 million for the USA population (82% of them biallelic) and 2.7 million for the Russian population (93% biallelic). 25% of the USA SNPs were shared with the Russian population (11% with the same major and minor alleles), and 53% of the Russian SNPs were shared with the USA population (23% with the same major and minor alleles, Appendix 3- figure 1).
The phylogeny of the sequenced genomes was reconstructed with RAxML (Stamatakis, 2014; Appendix 3-figure 2). Nucleotide diversity (π) was estimated as the average frequency of pairwise nucleotide differences; π for different classes of sites is shown in Appendix 3-figure 1B. Two samples from Florida (USA population) were excluded from the further analysis to minimize the possible effect of population structure.

Data on H. sapiens and D. melanogaster populations
We used polymorphism data from 347 phased diploid human genomes from African and 301 genomes from European super-populations sequenced as part of the 1000 Genomes project (1,000 1000Genomes Project Consortium et al., 2015. If several individuals from the same family were sequenced, we included only one of them. As a D. melanogaster dataset, we used 197 haploid genomes from the Zambia population (Lack et al., 2015). Only autosomes were analyzed in both datasets.

Estimation of LD
As a measure of linkage disequilibrium between two biallelic sites, we used r 2 , calculated as follows: , where p(A) and p(B) are the minor allele frequencies at these sites, and p(AB) is the frequency of the genotype that carries minor alleles at both sites. Singletons (sites with minor allele present only in one genotype) were excluded from the analysis if not stated otherwise.

Haploblocks annotation
In order to annotate the haploblocks, we calculated LD along the S. commune genome in a sliding window of 250 nucleotides with a step of 20 nucleotides (only non-singleton SNPs are analyzed; the windows with less than 10 SNPs were excluded). Any continuous sequence of overlapping windows with r 2 larger than the threshold value was merged together in a haploblock. The LD threshold value was defined independently for each S. commune population as the heavy tail of the within-window LD distribution, as compared with the lognormal distribution with the same mean and variance as in the data (Figure 3-figure supplement 7).

Estimation of LD between physically interacting amino acid sites
Of 16,319 annotated protein-coding genes of S. commune (Ohm et al., 2010) 9,941 were found in all 55 aligned genomes. We blasted the protein sequences of these orthologous groups against the PDB database of protein structures. About 52% of them (5,188) had a match (e-value threshold = 1e-5) amongst the proteins with the known structure. We realigned the sequences of S. commune protein and the matching PDB protein with clustal and calculated within-population LD and physical distance (Å) for each pair of aligned positions. A pair of amino acid sites was considered physically adjacent if they were located within 10 Å from each other.
To compare LD between pairs of physically close and distant sites, we used the controlled permutation test (Figure 2A): for each pair of physically close amino acid sites (within 10 Å) we sampled a pair of physically distant amino acids on the same genetic distance (measured in aa). Pairs of sites closer than 5 aa were excluded from the analysis.
To examine LD patterns within individual protein structures, we calculated contingency tables of pairs of SNPs being located in codons encoding physically close amino acids and having high LD (no less than 90% quantile for a given gene). Pairs of amino acid sites located closer than 30 aa or more distant than 100 aa from each other were excluded; genes with less than five pairs of physically close sites under high or low LD were also excluded. From these contingency tables, we calculated the odds ratio (OR) and chi-square test p-value for each gene. p-values were adjusted using BH correction. We identified 22 genes with pairs of adjacent sites having significantly higher LD in the USA population (out of 1286 eligible genes in total), and 87 genes in the Russian population (out of 967) under 5% FDR (Appendix 3-table 1). Examples of such genes are shown in Figure 2 and

Simulations of epistasis
To simulate evolution of populations with or without epistasis and balancing selection (Figure 3figure supplement 6), we used an individual-based model implemented by SLiM (Haller and Messer, 2019). Simulations were performed with diploid population size N=1000 and recombination rate 0. To achieve the level of genetic diversity π similar to S. commune, mutation rate μ was scaled as μ=π/2N=5e-5. The length of the simulated sequence was 100 nt. Each simulation started with a monomorphic population and proceeds for 100 N generations. For calculations of synonymous and nonsynonymous LD, random 100 haploid genotypes were sampled from the population. Only SNPs with minor allele frequency >5% in the sample were analyzed.
We modelled two types of mutations, depending on whether they are neutral (with selection coefficient s syn = 0) or weakly deleterious (s nonsyn ≤0), representing synonymous and nonsynonymous variants correspondingly. There were twice as many nonsynonymous as synonymous sites. Under the nonepistatic model, s was independent of the genetic background. We assumed s nonsyn = −0.01 with the dominance coefficient h of 0.5.
Under the pairwise positive epistasis model, we assumed that one nonsynonymous mutation can be partially or fully compensated by a mutation at another site. In this model, all nonsynonymous sites were split into pairs. Each mutation of a pair individually occurring within a genotype was assumed to be deleterious, with selection coefficient s nonsyn = −0.01; however, the fitness of the double mutant is larger than expected under the additive (non-epistatic) model. We used several models of epistasis, with different strengths of epistasis strength and landscape shapes (Figure 3-figure supplement 6).
In the NFDS model of balancing selection, a single mutation at a random position was subjected to frequency-dependent selection (so that it is positively selected at frequencies below 0.5, and negatively selected at frequencies above 0.5). In the AOD model, mutations in 10 random positions were fully recessive (h=0) and weakly deleterious (s=−0.0025).
To simulate evolution of populations with different levels of genetic diversity under epistasis (Appendix 1), we used FFPopSim (Zanini and Neher, 2012). To achieve different levels of genetic diversity π, mutation rate μ was scaled as μ=π/2N. The calculations were performed the same way as in SLiM, but In this case, we used haploid population size N=2000, population-scaled recombination rate 0.01 and the simulated sequence length of 300 nucleotides.

Funder
Grant reference number Author

14-50-00150 Georgii A Bazykin
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Additional files
The following dataset was generated: Author ( Appendix 1 Epistatic selection is more efficient in genetically diverse populations Genetic interactions affect operation of selection only in a sufficiently variable population. The potency of any kind of selection increases with the amount of variation; for epistatic selection, however, this increase is expected to be faster than linear, because it depends on the number of possible allele combinations. In a highly polymorphic population, a particular allele is more likely to co-occur in the same haplotype with an interacting, for example compensatory, allele (Appendix 1- figure 1A) which should increase the impact of epistasis on linkage disequilibrium. To illustrate this point, we modelled the evolution of a genome region in the presence and in the absence of positive epistasis in a panmictic population. We assumed that all mutations at a set of sites are individually deleterious, and that all these sites are involved in pairwise positive (i.e. antagonistic) sign epistasis; specifically, each deleterious mutation can be fully compensated by another mutation at exactly one site elsewhere in the genome, which is also deleterious when present alone. In the non-epistatic simulations, the effects of mutations were independent; however, at the end of the simulation we randomly assigned the 'interacting' pairs of sites to account for the random coincidence of deleterious alleles. We found that in this model a higher polymorphism increases the probability that a deleterious mutation is compensated before being eliminated by selection (Appendix 1-figure 1B). This probability increases with genetic diversity even for the nonepistatic simulations, because increased diversity elevates the likelihood of randomly encountering a compensating allele in the same haplotype. For epistatic simulations, however, this increase is more radical, reflecting the effect of epistatic selection favoring compensated haplotypes.
After the mutation-selection equilibrium was reached, we measured the strength of epistatic selection between all segregating polymorphisms, asking to what extent the mutational load is reduced by epistasis maintaining combinations of compensatory mutations. As shown in Appendix 1- figure 1C, the ability of epistatic selection to reduce the mutation load (i.e. to increase the mean fitness) strongly depends on π. In less variable populations (π<0.01), epistasis is practically inefficient and does not affect LD (Wilcoxon test -values >0.33); this is because the probability of occurrence of the favorable combination of alleles in the population for selection to act upon is low. In more diverse populations, however, such combinations may arise and be favored by epistatic selection, which increases LD between them (Wilcoxon test p-value <0.01 for π≥0.01). First, random drift, which affects nearly-neutral synonymous sites more than nonsynonymous sites which are mostly under negative selection, leads to attraction between minor alleles (Sandler et al., 2021). Second, negative selection at nonsynonymous sites causes repulsion between rare, deleterious alleles, due to Hill-Robertson interference, even if this selection does not involve any epistasis (Comeron et al., 2008;Garcia and Lohmueller, 2021;Hill and Robertson, 1966; Appendix 2- figure 4). Third, there are data on negative epistasis in this selection, which also should lead to repulsion of deleterious alleles and, thus, negative LD between rare nonsynonymous alleles (Garcia and Lohmueller, 2021;Sandler et al., 2021;Sohail et al., 2017). The first and the second factors are weak and can produce noticeable LD only between tightly linked loci, while the third factor may generate even long-range LD. By contrast, LD nonsyn >LD syn can be explained only by positive epistasis in selection at nonsynonymous sites.
Although negative selection generally results in LD nonsyn <LD syn , our simulations demonstrated that Hill-Robertson interference without epistasis can produce attraction between minor alleles under a rather restrictive set of conditions. In these simulations, deleterious polymorphisms can achieve high frequency only in regions of low recombination, leading to LD nonsyn >LD syn for extremely high MAF (Appendix 2- figure 5A). However, this effect does not hold if assuming unequal fitness effects of deleterious mutations (Appendix 2-figure 5B) or while merging SNPs of different frequencies together (Appendix 2- figure 5C).