Structural properties of genotype–phenotype maps

1. Introduction

Reproduction, as a fundamental property of biological systems, depends on the storage, processing and transfer of biological information. That information is typically stored in the form of sequences, such as DNA, RNA or amino acid sequences, and is more generally referred to as the genotype. In abstract models of biological evolution, the genotype can take more general forms, such as for genetic algorithms [1], where the genotypes are often binary strings. Some have chosen network topologies [2] or Boolean update rules [3] as genotypes. But genotypes are almost always linear and discrete representations of biological information. By contrast, the definition of a phenotype is almost impossibly broad. At every level of resolution of biological structure or function, any higher-level outcome resulting from a sequence could be viewed as a phenotype. For example, the amino acid to which a given triplet codon maps can already be considered a phenotype. The structure, interactions and functions of RNA and proteins are phenotypes, as is the function of a metabolic circuit or the state of a biological cell. On a longer time scale, the development of a type of tissue or of a whole organism can also be viewed as a phenotype, as can the behaviour of an organism and its interaction with its ecosystem. This complexity explains the wide range of genotype–phenotype (GP) maps found in the literature. It also explains why many GP maps that are studied are abstract theoretical or computational models, as the space of biological possibilities is often intractably vast, and the hierarchical levels of biological complexity are numerous. An illustration of the difficulties associated with the study of real biological GP maps can be found in the form of the protein folding problem [4], which has occupied biologists for decades, and which represents one of the most immediate connections between genotype and phenotype—the spatial rearrangement of an amino acid sequence into a protein structure.

It has been known since Mendel [5] that genetic mutations can cause phenotypic changes. Building on Mendel's and Darwin's [6] work, Fisher, Haldane and Wright produced the modern evolutionary synthesis [7–9] in the early twentieth century, which also introduced the idea of the fitness landscape [9]. A fitness landscape relates the space of genotypic variation to survival. It, therefore, contains a GP map implicitly. But the fitness landscape really consists of two distinct mappings: one is the GP map, and the other is the mapping from phenotypes to fitness values. A seminal contribution to the unravelling of these two mappings was made by Kimura [10], almost 50 years ago, who postulated that many mutations that are important for evolution must be neutral, meaning that they do not affect the fitness of the phenotype. While neutrality refers to selection rather than phenotypic change, and while it is possible for different phenotypes to be equally fit, the ubiquity of neutrality proposed by Kimura made it reasonable to suppose that many mutations are not only neutral in terms of selection, but in fact leave the phenotype entirely unchanged.

This idea was substantially extended by Maynard Smith in 1970 [11], who addressed the apparent contradiction [12] between the vast number of possible amino acid sequences and the tiny fraction of these sequences that give rise to the proteins observed in nature. If the evolution of these proteins was driven by natural selection, how would this evolutionary process have found these solutions through random mutations? Maynard Smith's answer was to postulate that ‘functional proteins must form a continuous network which can be traversed by unit mutational steps without passing through non-functional intermediates’ [11], which is similar to the definition of what is now commonly referred to as a neutral network. Such networks mean that functional proteins occupy connected subsets of genotype space, which makes their discovery through an evolutionary process, driven by random mutations, feasible. Maynard Smith arguably laid the foundations for the modern study of GP maps in this paper by proposing the concept of a ‘protein space’ of all possible amino acid sequences, in which neighbours are defined by single amino acid substitutions. This construct is almost exactly the same framework that has been used to study GP maps and their properties in RNA, proteins and many other systems, mainly from the late 1980s onwards. An explicit discussion of GP mappings can be found in R. C. Lewontin's 1974 work The genetic basis of evolutionary change [13], but Maynard Smith's earlier work already marks the separation of the fitness landscape from the GP map, as it considers the set of genotypes that map to functional proteins in general. While fitness is very hard to quantify in a biological system, properties of GP maps are much more easily quantifiable. The barrier that remains is the vast size of the space, but the exponential increase in both computational power [14] and the availability of biological data [15] means that ever more complex GP maps can now be studied quantitatively.

2. The RNA secondary structure genotype–phenotype map

While the earliest theoretical work on GP maps [10,11], and the earliest abstract models of such maps [16] considered protein tertiary structure, one of the most widely studied GP maps is that of RNA secondary structure, which is the specific configuration of base pair interactions between different parts of a given RNA strand. This is because the folded structure for a given RNA sequence can be predicted computationally to a sufficiently high degree of accuracy [17], which makes large-scale surveys of this GP map possible. Furthermore, the phenotype can be characterized precisely in terms of the base pairs that are formed in the folded structure, and can be denoted succinctly by using the ‘dot-bracket’ notation. In this notation, a pair of parentheses represents a base pair, while a period represents an unbound base. The GP map therefore can be viewed as the map from an RNA sequence, composed of the four bases A, C, G and U, such as:

to a dot-bracket sequence, such as:

The structure that corresponds to this sequence is shown in figure 1. Schuster et al. [18] in 1994 were the first to study global properties of the RNA secondary structure GP map using computational methods. They established that the number of RNA sequences, or genotypes, vastly outnumbers the number of RNA secondary structures, or phenotypes. In addition, they found that the number of sequences per structure followed a generalized form of Zipf's law, meaning that many structures mapped to only a few sequences each, while a few structures occupied very large proportions of sequence space. Furthermore, it was found that these common phenotypes could be reached in only a small number of mutational steps, starting from any random genotype. Lastly, this investigation also showed the existence of long neutral paths in genotype space by demonstrating that 21.7% of neutral paths ended on genotypes that had no base in common with the reference sequence used as the starting point. This result implies the presence of extensive neutral networks. It has been proposed [19] that these neutral networks arise as a percolation effect, meaning that a connected network emerges spontaneously for sufficiently high link densities, purely based on the probabilities of individual links. As some phenotypes have a large number of genotypes, such an explanation of neutral networks in terms of percolation would imply that neutral networks result from the Zipfian distribution of the number of genotypes per phenotype. Random graph theory calculations based on this argument yield a theoretical percolation threshold of 0.37 [19] for four-letter alphabets (such as in RNA), meaning that a phenotype occupying 37% of genotypes will percolate across the entire genotype space. Very recent work by Greenbury et al. [23], however, has established that the genotypes belonging to a given phenotype are in fact much more correlated than random graph theory would suggest, and that neutral networks are therefore even more ubiquitous and pervasive. These results will be discussed in more detail in later sections of this paper.

3. The HP model genotype–phenotype map

The most obvious map between sequence and structure in modern biology, and the one that Maynard Smith addressed in his seminal work [11], is that from an amino acid sequence to a protein structure. The alphabet of amino acid sequences is much larger than that of RNA sequences, as there are 20 amino acids that appear in protein sequences, compared to only four bases in RNA. The length of amino acid sequences tends to also be considerably larger than the length of RNA sequences that form functional secondary structures in the cell. The former typically consist of several hundred residues, while many examples of the latter are less than 50 bases in length [24]. In addition, the folding process in protein tertiary structure relies on a complex network of residue–residue interactions, whereas the folding of RNA secondary structure is primarily governed by hybridization, which is a much simpler mechanism. For these reasons, the prediction of protein structures from amino acid sequences is a much greater challenge than the prediction of RNA secondary structure. In 1985, Dill [16] considered an abstract model of protein sequences that focused on a distinction between hydrophilic and hydrophobic residues. The latter would form the core of globular proteins in this model, which led Dill to the conclusion that ‘the number of accessible conformations in the globular state’ was an ‘exceedingly small fraction of the number accessible to the random coil’. In other words, the phenotypes that correspond to globular proteins occupy a small fraction of the total space of possible phenotypes. Building on this work, Dill & Lau [20] introduced the HP model of protein folding, which distinguished hydrophobic (H) and polar (P) residues that fold on a two-dimensional lattice (figure 1). This level of abstraction made it possible to study the space of conformations exhaustively for short chain lengths, and resulted in the work of Lipman & Wilbur [25], published shortly afterwards, who studied the GP map of the HP model extensively. Their findings included the discovery of the large neutral networks, similar to the connected ‘protein space’ Maynard Smith originally envisaged, which led them to conclude that ‘perhaps the most striking observation is the critical role of neutral mutations in traversing this evolutionary space’, and to raise the question ‘whether this particular phenomenon is a general property of evolutionary spaces’ [25]. A more general version of this question, namely to what extent the structural properties of GP maps are universal in biology, has received increasing attention in the recent literature on GP maps [22,23,26,27], and perhaps the most central question of this review. As will be discussed below, the likely answer to Lipman & Wilbur's question is positive—that properties such as extensive neutral networks are indeed universal properties of evolving systems. The reason for this may lie in the fundamental way in which biological information is organized [27].

The tractability of the HP model enabled further studies of its GP map, such as the work by Li et al. [28] who found that the number of genotypes per phenotype in the HP model displayed a highly skewed distribution, similar to the Zipfian distribution found for RNA [18]. This observation was described using the notion of designability, meaning that HP model proteins with a large number of genotypes are highly ‘designable’, which makes them more easily accessible in an evolutionary process.

It has been proposed [29] that the neutral networks of different protein phenotypes are largely disconnected from each other, resembling a ‘plum pudding’ [29] in contrast with the RNA GP map, whose interconnected phenotypes could be described with a ‘spaghetti’ metaphor [29]. Recent work [23], discussed in more detail below, has shown, however, that this strongly depends on the type of HP model used. The compact HP model, which only considers folded structures, exhibits interconnected phenotypes much like RNA.

4. The Polyomino genotype–phenotype map

In RNA secondary structure and protein tertiary structure, the most studied GP maps, the transformation from genotype to phenotype is achieved through a folding process that immediately links the genotype sequence to the final phenotype. Above individual proteins, the next level in the hierarchical organization of biological structures is protein quaternary structure, which describes the binding of multiple proteins to each other, for example, in the form of protein complexes [30,31]. Such structures differ from RNA secondary structure and single proteins because they require several independent subunits to come together and self-assemble into a larger structure. Recent work by Ahnert et al. [21] introduced a two-dimensional lattice model of self-assembly that can also be used as a GP map to model protein quaternary structure [22,32]. In this model, two-dimensional square tiles can bind to each other, forming larger connected structures called ‘Polyominos’. In order to use this model as an abstract GP map for protein complexes, we can encode the configuration of binding interactions on a set of tiles as a string, or genotype, and define the final assembled ‘Polyomino’ shape on the two-dimensional lattice as the corresponding phenotype (figure 1). The Polyomino GP map model has been found to exhibit a rich and complex structure while still being tractable enough to explore GP maps and fitness landscapes exhaustively [22,32]. Furthermore, it exhibits many of the same properties as the RNA and HP model GP maps, such as a large number of genotypes relative to the number of phenotypes, a highly skewed distribution of the number of genotypes per phenotype, and extensive neutral networks [22], among several other properties that will be discussed in more detail below.

5. Deleterious phenotypes in RNA, the HP model and Polyominoes

In all three GP maps discussed so far—RNA secondary structure, the HP model of protein structure and the Polyomino model—a large proportion of genotype space maps to phenotypes that are in some sense not viable, and therefore deleterious. In RNA, this corresponds to the case when the secondary structure without any paired bases has the lowest free energy, meaning that the RNA molecule remains unfolded. This is the case for 85% of RNA sequences of length L = 12 [23]. This proportion drops with increasing sequence length, but even for sequences of length L = 20, a third of RNA sequences still remain completely unfolded [23].

For the non-compact HP model, which considers all HP sequences of a given length, 98% of sequences of length L = 24 do not have a unique ground state structure. The absence of a unique ground state is conventionally regarded as a deleterious phenotype. For the compact HP model, specifically for the sequences on a 5 × 5 lattice, the fraction of such phenotypes is 82% [23]. Despite this high proportion of deleterious phenotypes, the neutral spaces of the 5 × 5 compact HP model are highly connected to each other, thereby conforming more to the ‘spaghetti’ than the ‘plum pudding’ metaphor [23] discussed earlier. Unlike RNA secondary structure, the HP model retains a high, and possibly constant [33] fraction of unfolded phenotypes with increasing sequence length.

In the Polyomino phenotype, certain configurations of interfaces on the square lattice tiles can lead to unbound growth. Furthermore, some configurations are non-deterministic, meaning that different structures are built in a stochastic assembly process, because the same pair of building blocks can interact in multiple distinct ways. Unbound and non-deterministic (UND) building block sets are categorized as a single deleterious phenotype, similar to the RNA strands that remain unfolded and the HP sequences without a unique ground state. In the Polyomino models most commonly studied, the proportion of the UND phenotype is 54% (two tiles, eight interaction types) and 80% (three tiles, eight interaction types).

In all of these three GP maps—RNA, HP 5 × 5 and Polyominoes—the deleterious phenotype is underrepresented in the neighbourhood of non-deleterious phenotypes, which illustrates that the neutral networks of different phenotypes are connected to each other even when a large proportion of genotype is occupied by the deleterious phenotype [23].

6. Network genotype–phenotype maps

The previous sections have given an overview of GP maps between sequences and biological structures. The definition of a phenotype is of course much broader than this. Biological interactions, such as gene regulatory interactions, metabolic interactions and signalling networks, offer an alternative way to characterize phenotypes, and are discussed in this section.

7. Direct measurements of biological genotype–phenotype maps

All GP maps discussed so far are to a greater or lesser extent abstractions of real biological systems. The RNA secondary structure GP map is the closest to a real biological GP map, as small RNAs fulfil a wide range of functions in the cell, and the map has been studied in the context of this naturally occuring variety [24]. Nevertheless, the phenotype representation in terms of the base pairs only captures part of the complex three-dimensional structure of folded RNAs. A level of abstraction is inevitable if the goal is to study general properties of GP maps, which at the level of protein tertiary structure remain intractable due to the enormous complexity of the folded amino acid chain, and at the level of interaction networks remains challenging due to the difficulty of reliably measuring and verifying the large number of interactions between proteins, genes, DNA and RNA in the cell. There have, however, been attempts to measure ‘real’ GP maps of protein domains [46] and of transcription factor binding sequences [47] to see whether they, at least within the limited view that existing sequence data affords us, share the same properties as the more abstract, comprehensive GP map models. The GP map of single-domain protein sequences uses the enzymatic function classification as the phenotype [46]. Because this is a functional classification (in terms of the enzyme commission classification [48]), and not the domain structure itself, there is an additional mapping between structure and function that complicates the GP map. Several different structures can fulfil the same function—a phenomenon that the authors of [46] refer to as structural promiscuity. While 86% of functions are only carried out by one structure, some functions (such as DNA-directed DNA polymerase) can be carried out by up to 14 different structures. In this GP map, the sets of enzymatic functions in the neighbourhoods of two genotypes diverge rapidly with increasing distance between the genotypes. This is even true when considering genotypes that map to the same protein structure. The multitude of functions exhibited by some protein families (which group together proteins of the same structure) is the reason for this rapid growth in functional diversity with increasing genotype distance, or in other words ‘Functional neighbourhood diversity emerges from the multifunctionality of structures’ [46]. Another prominent example of an attempt to directly measure the properties of a real biological GP map is the work by Payne & Wagner [47] on the transcription factor binding sites for 104 transcription factors in mouse, and 89 in yeast. By examining the space of all binding sites of eight nucleotides in terms of an enrichment score [49] that represents the binding affinity, one can construct the genotype network and observe where the sequences that bind to a given transcription factor fall on this network. In almost all cases, the majority of bound sites for a particular transcription factor form a single connected component, and, for a majority of transcription factors, all bound sites are part of a single connected genotype network. The binding site sequences are significantly more robust than expected by chance, and also significantly more evolvable (measured by considering how many mutations were required to reach a binding sequence for a different transcription factor). Both robustness and evolvability increase with the size of the binding site repertoire. These findings mirror the results regarding robustness and evolvability found in more abstract model GP maps of RNA secondary structure, the HP model, Polyominoes and network phenotypes.

The observation that many biological phenotypes are highly robust has been also made on a more general level, supporting the more detailed investigations of GP maps discussed so far. For example, robustness may explain why similar gene expression patterns are observed across a wide range of species, even when the cis-regulatory regions are highly divergent [50]. A similar result was found across a range of Drosophila species [51].

8. Structural properties of genotype–phenotype maps

The previous sections provide an overview of research into the properties of GP maps, and show that certain characteristics of the distribution of phenotypes over the space of genotypes have been observed across many different phenotypes and modelling approaches. In this section, I explore these ‘structural’ properties in more detail, particularly in the context of relatively tractable GP map models that can be studied on a global level, such as the RNA secondary structure, HP model and Polyomino GP maps. The first three of these, redundancy, bias and robustness, are illustrated in figure 2.

8.3. Robustness

Robustness is typically measured in terms of the fraction of possible mutations that leave a phenotype unchanged. This measurement can be applied at the genotypic level or at the phenotypic level [52]. In the latter case, the robustness of a given phenotype can be calculated as the average of genotypic robustness over the genotypes that map to that phenotype [52]. Robustness as a local measurement of neutrality is a key property of neutral networks. A comparison of the number of genotypes that map to a phenotype—in other words the frequency of a phenotype—and its robustness shows that a phenotype's robustness ρ_p scales as the logarithm of its frequency f_p in RNA secondary structure [22,23,53,54], Polyominoes [22,23] and the HP model [23], as illustrated schematically in figure 3. Formally, we therefore have

where a and b are positive constants. If genotypes of the same phenotype were uncorrelated (while still following the same distribution in terms of bias), the robustness of a phenotype would simply be equal to its frequency [23]. The robustness of the most frequent phenotypes in the RNA, HP and Polyomino maps lies between 0.5 and 1, which means that the logarithmic increase in the robustness makes the phenotypes far more robust than they would be according to the linear increase of the null model. Across these GP maps, phenotypes are therefore significantly more robust that one would expect based on bias alone [23]. This also becomes clear if one considers the size of the largest neutral component of a given phenotype, and the number of its neutral components. In the null model, one sees the emergence of a giant component around a phenotype frequency of 0.37, which corresponds to the percolation threshold found by Schuster [19]. By comparison, the observations for the RNA, HP and Polyomino GP maps show that the sizes of the largest neutral components scale logarithmically with phenotype frequency, making them much larger than the null model expectation below the percolation threshold. Furthermore, at the percolation threshold the number of distinct neutral components in the null model, which grows linearly with phenotype frequency up to that point, drops rapidly, as expected. This again contrasts with the three observed GP maps, in which most phenotypes are divided into less than a hundred (and in many cases less than 10) components [23], rather that the many thousand components observed in the null model below the percolation threshold. These findings again show that genotypes are far more correlated than one would expect based solely on GP map bias. The consistency of the logarithmic relationship between robustness and frequency, across three entirely different GP maps and many orders of magnitude of phenotype frequency, strongly suggests that it is a fundamental property of GP maps. This logarithmic scaling has been discussed in the literature [53,54], and an analytical derivation of this relationship has been proposed, based on the scaling of the neutral space size as a result of the number of paired and unpaired bases in the sequence. We will return to this later on, in §9.

Figure 3.

Open in a new tab

9. Sequence constraints and genotype–phenotype map properties

The structural properties discussed above are remarkably consistent across different GP maps. This suggests that they are in some way fundamental to biological systems. A possible explanation for the universality of these properties is offered in [27], based on the fact that biological sequences are often divided into evolutionarily constrained and (relatively) unconstrained parts. In DNA, this distinction comes in the form of exons versus introns, and genes versus intergenic sequences. In RNA, this distinction can be made in terms of base pairs (which are more constrained) and loops (less constrained) [54]. In Polyominoes, one can distinguish interactive faces and neutral (or effectively neutral) faces. Furthermore, the level of constraint of a given sequence element can be changed by mutating the sequence. In DNA, a mutation can alter the position of the start and stop codons that determine the boundaries between coding and non-coding DNA. In RNA (or Polyominoes), a mutation can dissolve or establish a base pair (or an interaction between tiles), reducing or increasing the level of constraint at two sequence positions. The approach in [27] uses a very simple GP map, in which the binary genotypes are subdivided into a single ‘coding’ and a single ‘non-coding’ part by a simple ‘stop codon’ sequence. Any point mutation of the coding part of the sequence (including the stop codon) changes the phenotype, whereas the non-coding part is entirely free to mutate. This GP map can be used to study the effect of constrained and unconstrained sequences on the GP map structure. It is named the ‘Fibonacci GP map’ because it is closely related to the Fibonacci code in computer science, and because many of its analytical properties can be described in terms of Fibonacci numbers. The Fibonacci GP map exhibits all of the structural GP map properties discussed above—redundancy, bias, the logarithmic scaling of phenotype robustness as a function of phenotype frequency, the negative correlation of genotypic robustness and genotypic evolvability, the positive correlation of phenotypic robustness and phenotypic evolvability, and shape space covering. This implies that the structure of GP maps is heavily influenced by the fact that unconstrained subsequences of the genotypes result in phenotypes that occupy whole subspaces of genotype space, allowing them to form extensive neutral networks that combine robustness and evolvability. The logarithmic scaling relationship between phenotypic robustness and frequency is a particularly strong indicator that unconstrained genotypic sequences play an important role in determining the structure of GP maps. This aspect of the Fibonacci GP map is closely related to the analytical derivation of the same logarithmic relationship in RNA secondary structure [54], which is specifically based on the constraints imposed on RNA sequences by base-pairing interactions. The Fibonacci GP map also shows that the possibility of mutations that change the level of sequence constraint leads to much greater evolvability, and is likely to be one of the reasons why phenotypic robustness and evolvability are positively correlated. It is in fact the analytical term that accounts for mutations of the stop codon that turns the correlation between phenotypic robustness and evolvability from a negative one into a positive one. This is because mutations of the stop codon shift the boundary between the constrained and unconstrained parts of the sequence, enabling access to a much wider variety of different phenotypes.

10. Evolutionary implications of genotype–phenotype map structure

While selection remains the driver of phenotypic change, it has become increasingly clear how much the road map—in other words, the GP map—matters. Characteristics such as the skewed distribution of neutral network sizes and their shape in genotype space can strongly determine evolutionary outcomes. On an abstract level, this insight is already present in Iwasa's concept of free fitness [55]. But in terms of concrete GP maps evidence has accumulated only more recently, for RNA secondary structure [56] and regulatory network GP maps [57], that the neutral network size strongly determines how likely a phenotype is to arise in evolution. In the context of RNA structures, it was also shown more recently that fitter phenotypes remain undiscovered if their neutral network is too small [58]. The most striking result in this respect is that the functional RNA secondary structures that have arisen naturally in the course of biological evolution closely mirror the distribution one would expect based on the neutral network size [24]. This provides strong evidence that, rather than being the result of strong selective pressures and evolutionary optimization, functional non-coding RNAs represent the most accessible solutions to a given evolutionary challenge. More generally speaking this result implies that the accessibility of phenotypes may be just as important a determinant of evolutionary outcomes as the pressures of natural selection.

Data accessibility

This article has no additional data.

Competing interests

I declare I have no competing interests.

Funding

The author was supported by the Royal Society and the Gatsby Foundation.