Alphabet cardinality and adaptive evolution

One of the most fundamental characteristics of a fitness landscape is its dimensionality, which is defined by genotype length and alphabet cardinality—the number of alleles per locus. Prior work has shown that increasing landscape dimensionality can promote adaptation by forming new ‘uphill’ mutational paths to the global fitness peak, but can also frustrate adaptation by increasing landscape ruggedness. How these two topographical changes interact to influence adaptation is an open question. Here, we address this question in the context of alphabet cardinality, using theoretical fitness landscapes with tuneable fitness correlations, as well as three empirical fitness landscapes for proteins. We find that the primary effect of increasing alphabet cardinality is the introduction of a new global fitness peak. Controlling for this effect, we find that increasing alphabet cardinality promotes adaptation on uncorrelated fitness landscapes, but frustrates adaptation on correlated fitness landscapes. The primary explanation is that the increased ruggedness that accompanies alphabet expansion is characterized by an increase in mean peak height on uncorrelated fitness landscapes, but a decrease in mean peak height in correlated fitness landscapes. Moreover, in two of the empirical fitness landscapes we observe no effect of increasing alphabet cardinality on adaptation, despite an increase in the number of peaks and a decrease in mean peak height, calling into question the utility of these common measures of landscape ruggedness as indicators of evolutionary outcomes.


Introduction
The metaphor of the fitness landscape has framed evolutionary thought for nearly a century [1].In it, genotypes represent coordinates in a multi-dimensional genotype space, and fitness represents the 'elevation' of each coordinate in this space.This defines the surface upon which populations evolve, with natural selection favouring 'uphill' mutational steps, driving populations toward peaks of high fitness.This metaphor has inspired a diversity of theoretical models (e.g. the House-of-Cards, NK and rough Mt.Fuji (RMF) landscapes) [2][3][4], has been used to predict and direct adaptive evolution [5][6][7][8], and more recently has become the focus of many experimental studies, in systems ranging from macromolecules to metabolisms [9][10][11][12][13][14][15][16][17].
How an evolving population navigates a fitness landscape depends in part on population genetic conditions [18][19][20].In the strong selection, weak mutation regime, any beneficial mutation that appears and escapes stochastic loss will go to fixation prior to the arrival of a subsequent mutation, such that the population remains essentially monomorphic at all times [21].Adaptive evolution can then be thought of as a random adaptive walk, where populations traverse mutational paths in which fitness increases at each mutational step, ending on a local or global adaptive peak.Such paths are called accessible paths [14].
Many theoretical studies have focused on adaptive walks in fitness landscapes, revealing the important influence of topographical and topological properties, such as landscape ruggedness and dimensionality, on the dynamics and outcomes of evolutionary processes [9,[22][23][24][25][26][27].For example, as landscape ruggedness increases, the formation of maladaptive valleys renders many mutational trajectories to high fitness peaks inaccessible.In contrast, as landscape dimensionality increases, so-called 'extradimensional bypasses' emerge, which circumvent maladaptive valleys by forming indirect, accessible paths to high fitness peaks [28].How landscape ruggedness and dimensionality interact to open or close mutational trajectories to adaptation remains poorly understood.One of the reasons is that past research has tended to focus on bi-allelic landscapes (e.g. the presence or absence of adaptive mutations at each locus), controlling landscape dimensionality by changing genotype length L [23,24,28,29].However, empirical adaptive landscapes are not bi-allelic, because their genotypes are DNA or amino acid sequences.The number of alleles per locus is therefore an additional, and important, facet of landscape dimensionality [9,26].
The set of all possible alleles at any of a genotype's loci makes up an alphabet A, with cardinality |A| denoting the number of its constituent elements (figure 1).For instance, for DNA sequences, A = {A, C, G, T} and |A| = 4; for amino acid sequences, A is the set of 20 standard amino acids and |A| = 20.Allowable mutational transitions between alleles are defined by an 'allele graph' [27], such that for DNA sequences the allele graph is complete, because any nucleotide can mutate into any other, and for amino acid sequences the allele graph is defined by the standard genetic code, under which missense mutations can change an amino acid into six alternative amino acids, on average [30].Prior work on alphabet cardinality has revealed the formation of extra-dimensional bypasses that circumvent maladaptive fitness valleys [9], thus increasing the number of accessible paths to local and global fitness peaks [9,26], as well as the sizes of their basins of attraction [9], and the probability that at least one accessible path exists from an initial genotype to the global peak [26,27].However, because these studies did not systematically fine-tune alphabet cardinality on landscapes with different levels of fitness correlations, several open questions remain, specifically pertaining to the interplay between cardinality, ruggedness, and the outcomes of adaptive walks.Whereas increasing cardinality can increase the prevalence of extradimensional bypasses to the global peak, it can also change other aspects of landscape topography, such as the height, location, and number of local peaks [31,32].For instance, for the block model of NK landscapes, where each genotype is divided into B non-interacting blocks of size K + 1, the expected number of peaks is B , which grows approximately polynomially with increasing alphabet cardinality (for B ≪ L) and approximately exponentially with decreasing fitness correlations (i.e.B) [32].Such an increase in landscape ruggedness may offset any evolutionary benefit of extra-dimensional bypasses to the global peak if the landscape's local peaks are of low or intermediate fitness.Conversely, if these local peaks are of high fitness, their emergence may prove no obstacle to adaptation, even if they do titrate evolving populations from the global peak.Relatedly, increasing cardinality not only increases the number of extra-dimensional bypasses to the global peak, but also to local peaks [9,33], so whether evolving populations converge on the global or local peaks will depend on how these bypasses influence the relative sizes of the basins of attraction of these peaks.Finally, because extra-dimensional bypasses are by definition longer than direct paths, they are less likely to be utilized by evolving populations [34], which may instead follow direct accessible paths to emerging local peaks.Thus, calling a fitness landscape accessible based on the existence of at least one accessible path to the global peak [24,26,27] may belie a landscape's true accessibility.For these reasons, it remains unclear how and to what extent increasing alphabet cardinality facilitates adaptive evolution.
Here, we study the effect of increasing alphabet cardinality on fitness landscapes of different topographies.First, we consider a spectrum of theoretical landscapes with tuneable levels of fitness correlations, ranging from completely correlated Mt.Fuji landscapes to completely uncorrelated House-of-Cards landscapes.Systematically varying alphabet cardinality and controlling for the emergence of the global adaptive peak, we study the impact of alphabet cardinality on landscape topography, specifically the number and heights of fitness peaks, as well as the number of extra-dimensional bypasses to these fitness peaks.Further, we measure the effect of increasing |A| on the mean fitness reached by random adaptive walks-a quantity we call the final mean fitness.We show that landscapes with different levels of correlation show characteristically different behaviour upon increasing cardinality-whereas correlated fitness landscapes show a declining final mean fitness with increasing cardinality, uncorrelated fitness landscapes show an increasing final mean fitness with increasing cardinality.
Motivated by this observation, we then analyse the effect of increasing cardinality on three combinatorially complete empirical fitness landscapes with intermediate levels of fitness correlation.Using evolutionary simulations and controlling for the emergence of the global adaptive peak, we find that increasing cardinality does not lead to higher final mean fitness, even though it increases the number of accessible paths to the global peak.Whereas the results for the empirical fitness landscapes are qualitatively in agreement with the results for theoretical fitness landscapes with intermediate levels of fitness correlation, the theoretical fitness landscapes do not fully capture the behaviour of two out of the three empirical landscapes, thereby highlighting the limits of theoretical studies and the need for more experimental studies on fitness landscapes.

Theoretical fitness landscapes
To understand the effect of alphabet cardinality on landscapes of different topographies, we first analysed RMF landscapes [35,36].We generated the RMF landscapes by adding a normally distributed random variable with mean zero and variance r 2 to an additive Mt.Fuji landscape (see section 4), where the roughness parameter r controls the relative contribution of the random component and serves as a measure of the level of fitness correlation in the landscape.Thus, a landscape with r = 0 corresponds to a completely correlated additive landscape, whereas a landscape with r = 1 corresponds to an uncorrelated House-of-Cards-like landscape.
For proteins, the most obvious benefit of increasing cardinality is the expansion of the chemical toolkit of the amino acid alphabet, which may endow proteins with enhanced or entirely new functions [37,38].In the context of an adaptive landscape, this is equivalent to introducing a new fitness peak, which may be the global peak.Indeed, in RMF landscapes, increasing alphabet cardinality causes a monotonic increase in the global peak height (supplementary figure 1, top panel).To study the influence of alphabet cardinality on landscape topography beyond the effect of introducing a new global peak, we hold the global peak constant throughout every alphabet expansion.Specifically, we confine our analyses to a subset of the 20! Possible expansions of the alphabet in which the global peak genotype is always present.For example, consider a fitness landscape for a protein with L = 4 variable sites, in which the amino acid sequence WWLA is the fitness peak.In this case, we truncate the full landscape (i.e.|A| = 20) to a landscape with A = {W, L, A} and analyse its topography.We then incrementally expand A by randomly adding an allele from the remaining set and analyse the resulting landscape, continuing this process until all 20 alleles are included in A. To understand the general consequences of increasing |A| in RMF landscapes, we studied an ensemble of 100 randomly generated RMF landscapes and simulated 200 expansions of A for each of those landscapes.To enable comparison with the empirical protein landscapes in the following section, we present results for RMF landscapes with L = 4 in the main text and L = 3 in the supplementary text (supplementary figure 2).

Topographical analysis.
We first studied the topographical properties of each class of landscapes, namely the number of peaks, their mean height, and the number of accessible paths that terminate at the global peak, out of all accessible paths (see section 4).Except for the purely additive case with r = 0, the number of local peaks increases monotonically with increasing |A| (figure 2, first column).Note that the number of peaks also depends upon the size of the landscape and therefore, it is important to ascertain whether the increase in the number of peaks is due to the landscapes truly becoming more rugged or due to their increasing size alone.To do so, we also computed a landscape-size agnostic measure of ruggedness, the roughnessto-slope ratio, which measures the deviations of a fitness landscape from a linear model [39] (supplementary figure 3).The roughness-to-slope ratio also increases with increasing cardinality, confirming that the increase in the number of peaks is not solely due to the increasing landscape size.Further, regardless of the cause of the increase in the number of peaks, more peaks imply more sub-optimal endpoints for the adaptive walks, which frustrates adaptation.Thus, the number of peaks is indicative of the outcome of random adaptive walks, so we focus on this measure of ruggedness throughout the rest of the manuscript.
While the number of peaks increases monotonically regardless of the roughness parameter r, we observe that the trend in the mean fitness of these peaks depends on the level of correlation in these landscapes.For small values of r, i.e. strongly correlated landscapes, we observe monotonically declining mean peak height (figure 2(A), second column).In contrast, for large values of r, i.e. strongly uncorrelated landscapes, the mean peak height tends to increase with |A| (figure 2(D), second column).In line with prior literature [9,26], the number of accessible paths to the global peak increases with increasing alphabet cardinality (figures 2(A)-(D), third column).However, the fraction of accessible paths that terminate on the global peak, out of all accessible paths, decreases with |A|, regardless of the value of the roughness parameter r (except for the corner case of r = 0, in which case the fraction is independent of cardinality) (supplementary figure 4).The reason is that more local peaks emerge as alphabet cardinality increases, such that accessible paths are more likely to terminate on local peaks than on the global peak.In other words, although the absolute number of accessible paths to the global peak increases, the probability that a randomly chosen accessible path reaches the global peak decreases as more local peaks emerge in the landscape.
Why do the trends in mean peak height depend on the level of fitness correlation?For completely correlated landscapes (i.e.r = 0), the mutational effects of alleles at the various loci are completely independent of each other, so increasing |A| has no effect on the mean peak height because there is only one peak in the landscape (i.e. the global peak).By slightly decreasing the level of fitness correlation (e.g. when r = 0.05 or r = 0.125), we begin to introduce interactions between loci, yet the truncated landscape with only the global peak alleles remains very smooth, comprising only one or two peaks.Expanding the alphabet increases the number of sub-optimal peaks (figures 2(A) and (B), first column) and thus, the mean peak height decreases monotonically (figures 2(A) and (B), second column).In contrast, for landscapes with larger r, the truncated landscapes are already quite rugged and thus already have a low mean peak height.However, for intermediate levels of fitness correlation (e.g. when r = 0.4), we see a nearly constant mean peak height (figure 2(C), second column), while for uncorrelated landscapes (e.g. when r = 1.0) we find that it increases with |A| (figure 2(D), second column).This is because in uncorrelated landscapes, all the new mutational neighbours introduced due to the increased cardinality also compete to be peaks.Therefore, to be a peak, a genotype must be larger than all its (|A| − 1) • L neighbours.It is well known that the expected value of the maximum of (|A| − 1) • L + 1 normally distributed random variables increases with |A|, thus leading to higher mean peak height.The same constraint does not hold for landscapes with intermediate levels of fitness correlation, because a potential peak genotype does not have to compete with all its mutational neighbours; the existing correlations already ensure that some of them have lower fitness than the potential peak genotype.In other words, for uncorrelated landscapes, the average threshold fitness value to be a peak increases with increasing |A|, while it does not for landscapes with intermediate levels of fitness correlation.

Random adaptive walks.
In the previous section, we observed that the number of peaks and the number of accessible paths terminating at the global peak both increase with |A|, for all r (figure 2, first and third column), but that the trends in mean peak height depend on r.Whereas an increase in landscape ruggedness may frustrate adaptation by trapping evolving populations on local peaks, especially if these peaks are of low fitness, increasing the accessibility of the global peak may promote adaptation by opening new mutational paths to high fitness.Which of these two topographical changes matters more for adaptation?To find out, we simulated random adaptive walks in the RMF landscapes and recorded the final mean fitness for each value of |A|.We initialized the walks on all genotypes in the landscape, i.e. we simulated 20 4 = 160 000 walks.In each step of a walk, the current genotype was replaced by a randomly chosen neighbour of higher fitness.This process was repeated until the walk reached a peak (section 4).We then averaged over the fitness values reached by each of the walks to calculate the final mean fitness.We chose to simulate random adaptive walks, because they are more likely to employ indirect paths than greedy adaptive walks [9].
Figure 2 (fourth column) shows that the effect of |A| on the final mean fitness depends on the roughness parameter r.For the purely additive landscapes (r = 0), the final mean fitness is independent of |A|.For landscapes with small r, the final mean fitness decreases as a function of |A|; however, as r increases, the trend flips and increasing |A| leads to increasing final mean fitness (figure 2, fourth column).Interestingly, for each r, the final mean fitness curve resembles the mean peak height curve, except that it is shifted upwards (figure 2, second and fourth column).This suggests that the emergence of local peaks may have a stronger influence on adaptation than the formation of extra-dimensional bypasses to the global peak.But why are the curves for final mean fitness higher than those for mean peak height?

Basins of attraction
We reasoned that the curves for final mean fitness may be higher than those for mean peak height because the random walks preferentially converge on higher peaks.We therefore computed the fraction of random walks terminating on each peak in the landscape and used it as a proxy for the basin of attraction of that peak.
These data are shown in figure 3, which reveals four key trends.First, higher peaks tend to have larger basins of attraction, for all r and |A|.Indeed, the fitness rank of a peak and the rank of the size of its basin of attraction are always strongly correlated (even at |A| = 20, the Pearson correlation coefficient ranged between 0.7 and 1 for all r, with lower values for landscapes with larger r and perfect correlation for landscapes with r = 0.05 and r = 0.125).This is in line with prior work showing that genotypes are more likely to have an accessible path to a higher peak than to a lower peak [27].Second, for any |A|, the difference in the basin of attraction is larger for landscapes with small r than for landscapes with large r.Third, this difference decreases with |A| for all r.Fourth, for landscapes with small r, the basins of attraction of emerging local peaks grows, rather than shrink, as |A| increases.Together, these four trends explain why the curves for final mean fitness are higher than those for mean peak height.Had the basins of attraction been of equal size, these curves would be nearly identical.Indeed, the difference in the basins of attraction is smallest for landscapes with large r, where the final mean fitness curve most closely resembles the mean peak height curve, although it is still shifted slightly upwards due to the segregation of the basins of attraction based on peak height.Moreover, while for all r, the basin of attraction of the global peak decreases with |A| [31], the basins of attraction of the local peaks increase for landscapes with small r (figures 3(A) and (B)), but decrease for landscapes with large r (figures 3(C) and (D)).This implies that as |A| increases, local peaks become more important in landscapes with small r by attracting a higher proportion of the random walks, thus causing the final mean fitness to decline.
To further illustrate this point, we computed the entropy of the basins of attraction of all peaks in each landscape.This measure takes on its minimum value of zero when the landscape is single-peaked, and its maximum value of 1 when the landscape is multi-peaked and all peaks have basins of attraction of equal size.Supplementary figure 5 shows that entropy increases with |A| for all r.However, the rate of increase is highest for landscapes with intermediate r (e.g.r = 0.125), meaning that local peaks rapidly become significant attractors.Because these peaks tend to be of lower fitness than the global peak in landscapes with intermediate r, their emergence brings down the final mean fitness.In contrast, for landscapes with larger r, entropy is high and only increases marginally with |A|.Thus the basins of attractions are more uniform in size.However, entropy is always less than one, due to the persistent correlation of a peak's height and its basin of attraction.Taken together, these trends explain why the curves for final mean fitness track those for mean peak height, yet remain shifted upwards.

Empirical fitness landscapes
In the previous section, we saw that the effect of increasing cardinality on adaptation depends upon the level of fitness correlation in the landscape, which is controlled by the roughness parameter r.Empirical landscapes usually exhibit intermediate levels of fitness correlation [9,14,[40][41][42], and we focus on three such landscapes in the following section.

Data sets.
We used three empirical fitness landscapes built from combinatoriallycomplete data sets, i.e. data that characterize protein phenotype (e.g.binding affinity) for all possible 20 L combinations of amino acids at a small number L of protein sites [9,40].Following the protein evolution literature [9,43], we assume a direct mapping from phenotype to fitness, although this assumption is often violated [44].
The first data set pertains to streptococcal protein G domain B1.Wu et al [9] created a library of all possible combinations of amino acids at four sites in GB1 (20 4 = 160, 000 variants), specifically V39, D40, G41, and V54, which interact epistatically and influence the protein's binding affinity to immunoglobulin [15].They assayed protein phenotype-binding affinityfor all sequence variants by measuring the relative frequency of each variant before and after selection for binding immunoglobulin, defining protein phenotype as a log-enrichment ratio relative to the wild type (section 4).In the following, we refer to this data set as GB1.
The second and third data sets pertain to the bacterial ParD-ParE toxin-antitoxin system.Lite et al [40] created a library of all possible combinations of amino acids at three sites in the antitoxin ParD (20 3 = 8000 variants), specifically D61, K64, and E80, which influence the protein's specificity for ParE toxins.They assayed protein phenotype-the ability to antagonize toxin-for two toxins, ParE2 and ParE3, by measuring the relative frequency of each variant before and after selection for antagonizing toxin, defining protein phenotype as a logenrichment ratio relative to the wild type (see section 4).In the following, we refer to these data sets as ParE2 and ParE3, respectively.
For each of the three data sets we construct an adaptive landscape as follows: each protein sequence of a given length (L = 4 for the GB1 data set, and L = 3 for the ParE2 and ParE3 data sets) is represented as a vertex in a mutational network.Two vertices are connected by an edge if their corresponding sequences differ in exactly one position, i.e. we assume that each amino acid can mutate into any other [9,45].Each vertex is labelled with the phenotype (binding affinity or growth rate) of its sequence, defining its 'elevation' in the landscape.In order to reduce experimental noise, as well as impute missing values (6.6% of the GB1 variants failed to assay), we smoothed the data using empirical variance component regression [46] (section 4).The resulting landscapes comprised 17, 5, and 5 peaks for the GB1, ParE2, and ParE3 data sets, respectively.In the GB1 data set, the protein phenotype ranged from −8.28 (for genotype WPGI) to 2.52 (for genotype WWLA), while in the ParD datasets, the protein phenotype ranged from −7.78 (for genotype PIP) to 0.19 (for genotype ELK) in the presence of ParE2, and from −7.74 (for genotype APP) to 0.13 (for genotype DWE) in the presence of ParE3.

Topographical analysis and adaptive walks on empirical fitness landscapes.
The results for the topographical analysis of the three empirical landscapes are shown in figure 4 (columns 1-3).As expected based on the analysis of the RMF landscapes, the average number of peaks increases monotonically with increasing |A| for all three landscapes.Whereas the truncated GB1 landscape has approximately 2.5 peaks on average, the truncated ParE2 and ParE3 landscapes were almost always single-peaked.This reflects in the mean peak height, which shows little variation for the GB1 landscape, whereas it decreases for the ParE2 and ParE3 landscapes due to the emergence of local peaks.Finally, in accordance with the analyses on RMF landscapes, whereas the number of accessible paths to the global peak increases with |A| (figure 4, third column), the fraction of accessible paths that terminate at the global peak, out of all accessible paths, decreases (see supplementary figure 6).
Also in line with our analysis of RMF landscapes with intermediate r, we observed that alphabet cardinality has little impact on the final mean fitness for all three empirical fitness landscapes.This occurs despite a decrease in mean peak height in the ParE2 and ParE3 fitness landscapes (figures 4(B) and (C), second column), a trend we did not observe in the RMF landscapes.To explain this discrepancy, we next look at the basins of attraction of each of the peaks.

Basins of attraction
In figure 5, we show the basins of attraction of the top five peaks as a function of the alphabet cardinality, averaged over 3000 alphabet expansions.The results for the GB1 landscape (figure 5(A)) are similar to those obtained for RMF landscapes with intermediate level of correlation (figure 3(B))-fitter peaks have larger basins of attraction and the difference in the sizes decreases with increasing cardinality.The results for the ParE2 and ParE3 landscapes also resemble the L = 3 RMF landscape results (supplementary figure 7(A), but with one crucial difference: the basins of attraction of the lowest three peaks in the ParE2 and ParE3 landscapes  do not increase with increasing cardinality (figures 5(B) and (C)), a trend that we do not see in the RMF landscapes.This explains why despite showing very different trends in the mean peak height (figure 4, second column), all three empirical landscapes show similar trends in the variation of the final mean fitness (figure 4, fourth column).While in the GB1 landscape, this is due to the final mean fitness tracking the mean peak height, in the ParE2 and ParE3 landscapes, it is due to the sub-optimal peaks with very low fitness not affecting the final mean fitness, even though they significantly decrease the mean peak height, and the second highest peak having a reasonably high fitness.
Figure 6.The final mean fitness as a function of the alphabet cardinality for empirical fitness landscapes of protein (A) ParE2 and (B) ParE3 when the number and location of peaks were preserved.The black lines show the variation in the mean fitness averaged over 1000 possible expansions of the alphabet and the coloured lines show the variation for 50 randomly chosen expansions.We confined our analysis to only those alphabet expansions where the number of peaks was preserved.This was true for all ParE3 alphabet expansions and for approximately 68% of the ParE2 alphabet expansions.

Indirect paths promote adaptation when the number and location of peaks is preserved.
Our analyses of the RMF and empirical landscapes have revealed that increasing cardinality increases both the number of local peaks and the number of accessible paths to the global peak.However, our evolutionary simulations were more strongly influenced by the emergence of local peaks, such that final mean fitness tracked mean peak height.This led us to wonder if there exists a scenario in which the emergence of indirect paths to the global peak clearly promotes adaptation.
To find out, we considered an artificial scenario in which we controlled not only for the emergence of the global peak, but also for the emergence of local peaks.To do so, we began our alphabet expansions from a truncated alphabet that contained all of the amino acids in all of the peaks in the full landscape (i.e. with |A| = 20).This was not possible for GB1, because the peaks in the full landscape comprise all 20 amino acids.However, it was possible for both the ParE2 and ParE3 landscapes, because their peaks comprise only nine amino acids.
Figure 6 shows the outcomes of our evolutionary simulations on these landscapes.We observe a monotonic and non-saturating increase in final mean fitness with increasing alphabet cardinality.Thus, there exists a scenario in which the emergence of indirect paths to the global peak promotes adaptation, but this scenario is highly artificial.Under more realistic conditions, the emergence of local peaks that accompanies alphabet expansion has a stronger influence on adaptation than does the emergence of indirect paths to the global peak.

Discussion
Is dimensionality a blessing or a curse?For fitness landscapes, the answer may be a bit of both.On the one hand, increasing the length of a sequence or the cardinality of an alphabet can lead to the formation of new, high-fitness peaks or to the emergence of new mutational paths to high-fitness peaks.On the other hand, such increases in dimensionality can do the opposite; they can lead to the formation of new low-fitness peaks or to the emergence of new mutational paths to low-fitness peaks.Moreover, any increase in dimensionality necessarily increases the size of sequence space, which may slow adaptation if mutational paths to high-fitness peaks become very long.
For proteins, the most obvious advantage of increasing the size of the amino acid alphabet is expanding the chemical lexicon of the proteome.Different amino acids have different chemistries and adding to this chemical toolkit expands the space of functions available to proteins.Indeed, our analyses revealed that increasing cardinality caused a monotonic increase in the average height of the global fitness peak, revealing that the primary influence of alphabet cardinality on protein function is the introduction of protein variants with novel chemistries and improved functionalities.To probe the effects of alphabet cardinality on landscape topography beyond this primary effect, we studied an artificial scenario in which the global peak did not change during alphabet expansion.This has the advantage that it allows us to directly study secondary effects of alphabet expansion on landscape topography, but it has the disadvantage that it puts distance between our analyses and biological reality.An additional limitation of our study is that we did not take into consideration the genetic code, which influences both landscape topography [47] and the likelihood of various alphabet expansions [48][49][50][51].Moreover, we restricted our analyses to static fitness landscapes, although real fitness landscapes are often dynamic, for example due to fluctuating environmental conditions [44,52] or due to frequency-dependent selection [53].Finally, we confined our evolutionary simulations to the strong selection weak mutation regime and did not allow for downhill steps and crossing of fitness valleys.How relaxing this constraint affects our results is a direction for future work.With these caveats in mind, our results have four key implications for fitness landscape research.
The first implication concerns landscape accessibility, which is usually defined based on the existence of a single accessible path between a source and a target genotype [26,27,29].Accessibility, defined as such, has been shown to increase with increasing alphabet cardinality in both theoretical and empirical fitness landscapes [9,26,27], due to the emergence of so-called extra-dimensional bypasses [28].However, because these bypasses are necessarily indirect and therefore require more mutations than direct paths, their probability of realization can be very low.Moreover, there are typically many accessible paths emanating from any source genotype, and these often lead to non-target genotypes, such as those atop local fitness peaks.These observations indicate that landscape accessibility is more nuanced than existing metrics let on, and call for more holistic metrics that take into account a greater diversity of source-target pairs and the realization probabilities of individual accessible paths.
The second implication concerns the interaction between epistasis and alphabet cardinality, and its influence on adaptation.Past work has suggested that increasing alphabet cardinality may promote adaptation by increasing the number of accessible paths to the global adaptive peak [9,26], but may frustrate adaptation by increasing the number of local peaks in the landscape [54].How these two topographical changes come together to influence adaptation remained unclear.Our analyses of the RMF landscapes suggest that the latter topographical change has a stronger influence on adaptation, because the outcomes of our evolutionary simulations more closely tracked the average height of fitness peaks than they did the number of accessible paths to the global fitness peak.Importantly, how the average height of fitness peaks changed with alphabet cardinality depended on the level of epistasis amongst loci.When epistasis was low, the average height of fitness peaks decreased upon alphabet expansion; when epistasis was high, the average height of fitness peaks increased upon alphabet expansion.These results suggest that increased alphabet cardinality can either promote or frustrate adaptation, depending primarily upon its influence on the average heights of fitness peaks, a topographical property that is modulated by the extent to which loci epistatically interact.
Relatedly, the number and heights of fitness peaks are often used to characterize the ruggedness of a fitness landscape, and to make predictions about how adaptation may play out on a particular landscape [6,24,55,56].Indeed, in the RMF landscapes we found that the average heights of fitness peaks was a reasonable predictor of the outcomes of our evolutionary simulations.In contrast, in the ParD landscapes, we found that the average heights of fitness peaks decreased substantially upon alphabet expansion, but that this had little impact on the outcomes of our evolutionary simulations.The reason is that the three lowest local peaks that emerged during alphabet expansion had negligible basins of attraction, such that our evolutionary simulations preferentially converged on the global peak and the second highest peak.Therefore, the third implication of our study is that the number and heights of fitness peaks can be an insufficient proxy for landscape ruggedness, particularly in the context of predicting the outcomes of evolutionary processes.Alternative proxies that take into account the basins of attraction of fitness peaks and their entropy are likely to be more useful.
Finally, our results may have practical implications.Recent advances in biotechnology have enabled the design and construction of living organisms with expanded genetic codes, i.e. genetic codes that include a 21st, non-standard amino acid [38,[57][58][59].These expanded codes are commonly utilized in biophysical studies to elucidate the structure, function, and cellular localization of proteins (reviewed in [57]), and they are a promising tool in evolving and designing proteins with novel properties [59][60][61][62], including those with therapeutic applications [63][64][65].To date, nearly 250 non-standard amino acids have been explored for these purposes [66].By increasing the chemical diversity of the amino acid alphabet, these non-standard amino acids can introduce new high-fitness peaks in a protein's fitness landscape [67][68][69][70][71], such as those representing dramatic increases in thermostability [72,73].However, an important lesson from our study is that alphabet expansion can induce multiple topographical changes to a protein's fitness landscape, and some of these make finding the highest-fitness peaks more difficult.

RMF landscapes
To generate the RMF landscapes, we first generated an additive (Mt.Fuji) landscape of cardinality 20 by sampling the fitness effect of each allele at each locus from a uniform distribution between 0 and 1.The additive fitness value of any genotype σ was then just the sum of the contributions of the alleles at each locus i.e. f 0 σ = Σ L i =1 f i σi , where f i σi is the contribution of the allele σ i at locus i.Then to each genotype σ's additive fitness f 0 σ , we added a random component x that was scaled by the roughness parameter r, causing its total fitness to become with x ∼ N (0, 1).By tuning the parameter r, we could generate the entire spectrum of fitness landscapes with varying levels of correlation between fitness values and varying levels of epistatic interactions, with r = 0 leading to a purely additive landscape (with no epistatic interactions) and r ≫ 1 leading to a House-of-Cards landscape (with many epistatically interacting loci).All RMF landscapes were normalized to ensure that the fittest genotype had unit fitness.We used the same Mt.Fuji landscape (f 0 σ ) and added the random component x to it, for different values of r, in order to generate landscapes with varying levels of correlation that could still be compared with each other.

Data processing
For GB1, we downloaded the raw data from the supplementary tables in [9]; for ParE2 and ParE3, we obtained the data directly from the study authors.
We computed the fitness of each variant following Rubin et al [74].In particular, for GB1, the fitness of variant v is equal to For ParE2 and ParE3, two replicate measurements were reported for each variant.For each variant and each replicate, the fitness and variance were computed as described above for GB1.The fitness values of the two replicates were then combined to compute the fitness of variant v as a weighted average of the two replicates, with weights given by the inverse of the corresponding variance: and the variance of variant v was computed as , where f v,i and σ 2 v,i denote the fitness and variance of the ith replicate, respectively.These values of fitness and variance were used as an input for empirical variance component regression [46], a method to reduce experimental noise and impute missing variants.The smoothed landscapes were obtained as maximum a posteriori estimates.

Enumerating all accessible paths
To enumerate all accessible paths in the landscapes, we employed dynamic programming.In particular, the number of accessible paths starting in an antipodal sequence (relative to the global peak) and ending at sequence v was represented as vector P(v).The ith element of P(v), P(v) i , stores the number of accessible paths from any antipodal sequence to variant v that are of length i, i = 0, . . ., M, where M is the maximum possible length of the walks.The vectors were initialized to P(v) = (1, 0, . . ., 0) for the antipodal sequences, and to P(v) = (0, . . ., 0) for all other sequences.Iterating through the sequences in ascending order of their fitness, we updated P(v) using the following rule: where by N (v) we denote the set of all neighbours of sequence v, i.e. the set of sequences that differ from v in exactly one position, f (v) is the fitness of variant v, and I [ f(n) < f(v)] is an indicator function which is 1 if f(n) < f(v) and 0 otherwise.The number of accessible paths of length i ending at the global peak GP is then given by P(GP) i , the number of accessible paths of length i ending at a local peak is given by ∑ p∈local peaks P(p) i .For computational reasons, we computed the number of accessible paths only up to length M = 10.

Adaptive walk simulations
For the adaptive walk simulations, adaptation under strong selection and weak mutation can be treated as a Markov chain [18].We initialized the walks on each genotype in the landscape and used the equal fixation model [9,75] wherein, the next step of the walk was to a randomly chosen fitter neighbour of the initial genotype.This same process was iterated until no fitter neighbours could be found and a local or global peak had been reached.

Entropy of the basins of attraction
The entropy of the basins of attraction of all peaks in each landscape was computed by normalizing the Shannon entropy by the maximum possible Shannon entropy for a given number of peaks (n p ) i.e. log 2 (n p ).For landscapes that only had one peak, we omitted the normalization and let the entropy be zero.

Figure 1 .
Figure 1.Alphabet cardinality, extra-dimensional bypasses, and landscape topography.A two-locus genotype space with (A) A = {A, B}, |A| = 2; (B) A = {A, B, C}, |A| = 3; and (C) A = {A, B, C, D}, |A| = 4.In the middle row, the motifs show the underlying complete allele graphs for different A. The bottom row shows a fitness landscape for each |A|, in which nodes represent genotypes and edges connect genotypes that differ at a single locus, with arrows pointing from low-to high-fitness genotypes.Direct paths from the source genotype AA to the target global peak genotype BB are shown in red and are both inaccessible.Accessible indirect paths are shown in the remaining colors.As alphabet cardinality increases, the number of accessible paths to the global peak increases, as extra-dimensional bypasses open up, and the number of peaks (shown with underlined genotypes) increases from 2 to 3. Note that the location and the height of peaks can change as well, because the new peak genotypes may have different fitness values.

Figure 2 .
Figure 2. Topographical properties and evolutionary outcomes in relation to alphabet cardinality, for landscapes with different levels of fitness correlations.Each row corresponds to a different value of the roughness parameter r and each column to a topographical property or evolutionary outcome, which are indicated by column headers (these also define the corresponding y-axes).The number of accessible paths to the global peak only include paths starting from genotypes antipodal to the global peak.The coloured lines show results for 100 RMF landscapes and the black lines show the average over these landscapes.Data pertain to |A| ⩾ 5.

Figure 3 .
Figure 3. Basins of attraction of the (A), (B) top five and (C), (D) top ten peaks in RMF landscapes with varying values of the roughness parameter (A) r = 0.05, (B) r = 0.125, (C) r = 0.4, and (D) r = 1.0.Lines are rank ordered from top to bottom based on peak height.The results are averaged over 30 RMF landscapes.The insets in (C), (D) pertain to |A| ⩾ 15, showing how even at high cardinalities, the basins of attraction segregate based on fitness rank.Note the change in the y-axis limits across panels.

Figure 4 .
Figure 4. Topographical properties and evolutionary outcomes for three empirical fitness landscapes.Each row corresponds to a different landscape and each column to a topographical property or evolutionary outcome, which are indicated by column headers (these also define the corresponding y-axes).The coloured lines show results for approximately 10% of the alphabet expansions, which were randomly chosen, and the black lines show the results averaged over all 5000 expansions.

Figure 5 .
Figure 5.The basins of attraction of the top five peaks as a function of alphabet cardinality for empirical fitness landscapes of protein (A) GB1, (B) ParE2, and (C) ParE3.Lines are rank ordered from top to bottom based on peak height.For each alphabet cardinality, we generated 3000 random alphabet expansions.

,
where c v,sel is the count of variant v in the sample after selection for binding immunoglobulin, c wt,sel the count of the wild type in the sample after selection for binding immunoglobulin, c v,inp the count of variant v in the input sample, and c wt,inp the count of the wild type in the input sample.The variance of the estimate is given by