{'fraction_non_alphanumeric': 0.0308300395256917, 'fraction_numerical': 0.020796594709638187, 'mean_word_length': 4.90214140447422, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pairwise and higher order genetic interactions during the evolution of a tRNA

Summary A central question in genetics and evolution is the extent to which mutations have outcomes that change depending on the genetic context in which they occur1–3. Pairwise interactions between mutations have been systematically mapped within4–18 and between19 genes, and contribute substantially to phenotypic variation amongst individuals20. However, the extent to which genetic interactions themselves are stable or dynamic across genotypes is unclear21,22. Here we quantify >45,000 genetic interactions between the same 87 pairs of mutations across >500 closely related genotypes of a yeast tRNA. Strikingly, all pairs of mutations interacted in at least 9% of genetic backgrounds and all pairs switched from interacting positively to interacting negatively in different genotypes (FDR<0.1). Higher order interactions are also abundant and dynamic across genotypes. The epistasis in this molecule means that all individual mutations switch from detrimental to beneficial in even closely-related genotypes. As a consequence, accurate genetic prediction requires mutation effects to be measured across different genetic backgrounds and the use of higher order epistatic terms.

A central question in genetics and evolution is the extent to which the outcomes of mutations change depending on the genetic context in which they occur 1-3 . Pairwise interactions between mutations have been systematically mapped within 4-18 and between 19 genes, and have been shown to contribute substantially to phenotypic variation among individuals 20 . However, the extent to which genetic interactions themselves are stable or dynamic across genotypes is unclear 21,22 . Here we quantify more than 45,000 genetic interactions between the same 87 pairs of mutations across more than 500 closely related genotypes of a yeast tRNA. Notably, all pairs of mutations interacted in at least 9% of genetic backgrounds and all pairs switched from interacting positively to interacting negatively in different genotypes (false discovery rate < 0.1). Higher-order interactions are also abundant and dynamic across genotypes. The epistasis in this tRNA means that all individual mutations switch from detrimental to beneficial, even in closely related genotypes. As a consequence, accurate genetic prediction requires mutation effects to be measured across different genetic backgrounds and the use of higher-order epistatic terms.
To directly test the extent to which the effects of mutations and the interactions between mutations are stable or change depending upon the genotype in which they occur, we designed an experiment in which mutation effects and interactions are quantified across a large number of closely related genetic backgrounds. As a model system, we used the singlecopy arginineCCU tRNA (tRNAArg(CCU)) gene that is conditionally required for the growth of budding yeast (Extended Data Fig. 1a) and for which pairwise interactions have been previ ously mapped in one genetic background 15 . The small size of the gene allowed us to design a library that covered all 5,184 (2 6 × 3 4 ) genotypes containing the 14 nucleotide substitutions observed in ten positions in postwholegenomeduplication yeast species 26 (Fig. 1a, b). Each genotype therefore varies from zero to a maximum of ten nucleo tides divergence from the Saccharomyces cerevisiae tRNA sequence (Extended Data Fig. 1b). After transformation of the library into S. cerevisiae, we performed six selection experiments in parallel to quantify the relative fitness of each of the 5,184 variants under restrictive conditions (high temperature and 1 M NaCl) (Fig. 1c). The fitness of each genotype was quantified as the change in its abundance in each culture between the beginning and end of the competition period determined using deep sequencing with a hierarchical error model and normalized in log scale to the fitness of the S. cerevisiae genotype (henceforth 'fitness'). After filtering, we obtained fitness measurements for 4,176 variants (Supplementary Table 1) that corre lated well across replicates (Fig. 1d). The median fitness declines as the number of mutations increases but there are still many combinations of mutations with high fitness amongst genotypes that are far from the reference genotype (Fig. 1e).
We first examined the fitness consequences of single mutations and how these change across different genetic backgrounds (Fig. 2a). In the S. cerevisiae genotype, six of the 14 individual mutations were detri mental (Fig. 2b). However, when the same 14 mutations were made in the tRNA genotypes of the other six extant species (these alternative 'wildtype' tRNAs have fitness very close to the S. cerevisiae tRNA when expressed in S. cerevisiae, Supplementary Table 2), their effects changed substantially (Fig. 2b). For example, the mutation C66A had no effect in the S. cerevisiae background but became detrimental in the Candida glabrata tRNA, which only differs by two substitutions (paired ttest, q = 0.006, n = 6). Indeed, 11 out of 14 mutations had effects that changed across these seven tRNAs from different species (Extended Data Fig. 2a, false discovery rate (FDR) < 0.1).
We next compared the effects of the single mutations across the com plete set of genetic backgrounds in the library. In total, we tested each mutation in a median of 1,449 genetic backgrounds (minimum = 1,088, maximum = 1,993, Extended Data Fig. 1c, d). Notably, we found that every mutation was both detrimental and beneficial in a substantial number of genetic backgrounds (Fig. 2b, c, median number of back grounds in which the less frequent sign was observed = 6.4%; mini mum = 3.4%; maximum = 11.9% across all 14 mutations, FDR < 0.1, n = 21,450). Restricting the analyses to background genotypes with high or intermediate fitness, to genotypes with high input read counts, or to genotypes with few mutations did not change this conclusion (Extended Data Fig. 2b). Thus, all mutations have effects that switch from beneficial to detrimental in closely related genotypes.
To investigate the interactions between mutations that underlie these changes in mutation effects, we first quantified pairwise genetic inter actions between the 14 mutations, which is a total of 87 pairs in any genotype. We define epistasis as the difference between the fitness of each double mutant and the sum of the fitness of the two corresponding individual mutations. Consistent with previous results 15 , in the S. cerevisiae genotype, many pairs of mutations (40.2%, 35 out of 87) had combined fitness effects that were more detrimental than expected (negative epistasis) and only a few had effects that were less detrimental than expected (positive epistasis, 5.7%, 5 out of 87, FDR < 0.1, Fig. 3a). However, these interactions changed when they were tested in tRNAs from different species (Fig. 3b, c, Extended Data Fig. 3), with 83 out of the 87 interactions differing across the species (n = 1,000 paired ttests, FDR < 0.1, Extended Data Fig. 4).
We next analysed how the 87 interactions changed across all the genetic backgrounds in the library. Each interaction was quantified in a median of 506 genetic backgrounds (minimum = 240, max imum = 946, Extended Data Fig. 1d). Notably, all 87 interactions switched from positive to negative in a substantial proportion of the Letter reSeArCH genetic backgrounds (Fig. 3a). Restricting our analyses to genetic backgrounds with high or intermediate fitness, to combinations with high expected fitness or to genotypes with high input read counts did not change this conclusion (Extended Data Fig. 5b). Across all genetic backgrounds, positive and negative interactions were similarly preva lent (11.4% and 10.3% for positive and negative epistasis respectively, FDR < 0.1, n = 47,649).
Changes in base pairing only partially explained changes in the sign and magnitude of the effect of single mutations (Extended Data Fig. 6). The four pairs of mutations that restore Watson-Crick base pairs (WCBPs) were amongst the most robust positive interactions (Fig. 3e). However, even these combinations interacted negatively in a large frac tion of backgrounds (5.9-8.4%). This is consistent with the presence of nonWCBP nucleotides in these positions in the tRNAs from other    Letter reSeArCH species 27 (Extended Data Fig. 5c). Double mutants in the same RNA strand of the acceptor stem were enriched for negative epistasis (odds ratio (OR) = 1.23, Fisher's exact test P = 2.15 × 10 −6 , Extended Data Fig. 5d, e) and the restoration of a WCBP was also more likely to result in a negative interaction when the stem harboured multiple additional mutations in a single strand (Extended Data Fig. 5f). This suggests that other mechanisms, for example stacking interactions, are also impor tant determinants of tRNA function. We next tested whether pairwise interactions changed in back grounds containing each additional single mutation (Fig. 4a, Extended Data Fig. 7a). Notably, when averaging across genetic backgrounds, a total of 138 out of 316 possible thirdorder interactions were found (Extended Data Fig. 7b, FDR < 0.1), meaning that 76 out of 87 pairwise interactions were altered by the presence of a single additional mutation in the background (Fig. 4b). All 14 individual mutations altered at least eight pairwise interactions (median = 16.5, maximum = 24, Fig. 4c). As with secondorder interactions, thirdorder interactions were enriched amongst proximal mutations and mutations found on the same strand (Extended Data Fig. 7c, d).
However, as for pairwise interactions, all thirdorder interactions (316 out of 316) also switched from positive to negative across different genetic backgrounds, indicating the presence of even higherorder epistasis (Fig. 4d). 260 out of 316 thirdorder interactions changed in the presence of a fourth mutation (FDR < 0.1, n = 740). Indeed, interactions can be detected in this dataset up to the eighth order (Extended Data Fig. 7b, a total of 763 backgroundaveraged epistatic interactions from 3,961 possible interactions tested from order one to eight, FDR < 0.1). Consistent with the behaviour of the lowerorder interactions, the signs of many higherorder interactions also switch from positive to negative as the genetic background changes (Fig. 4d, 1,981 out of 3,691 interactions in the total dataset interact both posi tively and negatively in different genetic backgrounds at FDR < 0.1).
Finally, we evaluated the extent to which epistasis affected our ability to predict phenotypes from genotypes. We quantified the accuracy of genetic prediction in the 76 complete diallelic sublandscapes of eight mutations using models restricted to a single genetic background as a reference or models that averaged epistatic terms across multiple backgrounds (see Methods section 'Genetic prediction'). Although Letter reSeArCH individual mutation effects quantified in a single genetic background provide very poor prediction (Fig. 4e), the average effect of each mutation across all genotypes within a sublandscape improves the prediction (Fig. 4e, percentage of variance explained, PVE = 58% on heldout data, tenfold crossvalidation). Including a limited number of significant interaction terms further improves the prediction (Fig. 4f, Extended Data Fig. 8a, PVE = 64%). The best models evaluated by crossvalidation contain first and second order coefficients, but also higherorder interactions (Fig. 4g) that progressively reduce the pre diction error (Fig. 4h). However, these models contain a relatively small number of coefficients (20 out of 256 coefficients on average across sublandscapes, Extended Data Fig. 8b), suggesting that although pair wise and higherorder epistasis is important, reasonably sparse models can provide good genetic predictions when coefficients are measured across different genetic backgrounds.
Taken together, our results show that even single steps in sequence space substantially change the effects of both individual mutations and how these mutations combine to alter fitness. By a range of metrics, the combinatorially complete tRNA fitness sublandscapes are most sim ilar to rugged theoretical fitness landscapes 28 that constrain evolution (Extended Data Fig. 9). Indeed, the abundance of sign epistasis (Fig. 3d) limits the number of accessible evolutionary paths 29 , for example, paths between the genotypes of extant species (Fig. 4i, j, Extended Data Fig. 10). These results add to a growing body of evidence 2 that evolution is highly contingent at the molecular level. As a consequence, models that use coefficients averaged across different genetic backgrounds and that incorporate higherorder epistatic terms provide more accurate genetic prediction.

Online content
Any Methods, including any statements of data availability and Nature Research reporting summaries, along with any additional references and Source Data files, are available in the online version of the paper at https://doi.org/10.1038/s41586 01801707.  | Averaging coefficients across genetic backgrounds and using higher order epistatic terms is important for genetic prediction. a, Changes in the distribution of pairwise epistasis when the genetic backgrounds contain or do not contain the indicated mutation (left) and the distribution of the corresponding thirdorder epistasis values (right). b, Distribution of pairwise interactions that are altered by a third mutation. c, Distribution of single mutations that are involved in a third order interaction. d, Proportion of genetic backgrounds in which each combination of mutations from third to eighth order interact positively (orange) or negatively (green) at a FDR < 0.1. e, Agreement between observed and predicted fitness values of all eighth order complete sub landscapes (n = 19,456 genotypes, 76 sublandscape with 256 genotypes each) when using up to first order epistatic coefficients, relative to a single background genotype (left) or averaged across backgrounds (right, tenfold crossvalidation). f, Agreement between observed and predicted fitness values for all complete eighth order sublandscapes for the best models incorporating epistatic coefficients according to the rank of their significance and evaluated by crossvalidation (an average of 20 out of 256 epistatic coefficients per model). g, Mean orders of the most significant epistatic coefficients for the models used in f (bottom, relative to the possible number of coefficients per order; top, absolute counts). Error bars are 95% confidence intervals. h, Mean rootmeansquare error (RMSE) across the 76 eighth order sublandscapes when cumulatively adding the most significant coefficients determined by crossvalidation (inset, colour indicates the median order of the coefficient added across the 76 sublandscapes) or all significant coefficients from the same order (main). Error bars are 95% confidence intervals. i, Example of shortest paths between two extant species (top) and the accessible proportion (bottom). j, Average frequency of accessible paths between species.

MEthodS
Library design. tRNAs orthologous to S. cerevisiae tRNAArg(CCU) (encoded by HSX1) were collected from the Genomic tRNA Database 30 or extracted from the genome of each species using BLAST 31 ('blastall' 2.2.25). The sequences were aligned with Clustal Omega 32 . Across the 12 species closest to S. cerevisiae, only the six species shown in Fig. 1a had substitutions in the gene, with a total of 14 substitutions in ten positions. Allowing all of these substitutions to cooccur results in a total library size of 5,184 (2 6 × 3 4 ) possible mutation combinations. Plasmid library construction. An oligonucleotide of 115 nucleotides containing 72 nucleotides of tRNA flanked by 21 and 22 nucleotides of the yeast endogenous promoter and terminator was synthesized by IBA Lifesciences. At ten of the 72 positions of the tRNA, two or three different nucleotides were mixed in equal proportions during synthesis. For example, position one can be G or A, but position two can be T, G or C.
The oligonucleotide was amplified using PCR for ten cycles (Q5 Hot Start HighFidelity DNA Polymerase, NEB), purified using an Egel electrophoresis system (EGel SizeSelect Agarose Gel 2%) with column purification (MinElute PCR Purification Kit, Qiagen). Subsequently, the purified oligonucleotide was cloned into a version of the yeast centromeric plasmid pRS413 (HIS3 marker) 33 that contained the HSX1 gene flanked by 218 bp of upstream and 202 bp of downstream genomic sequences (pJD001). pJD001 was linearized from the HSX1 flanking regions (excluding the HSX1 sequence) using PCR (Q5 Hot Start HighFidelity DNA Polymerase, NEB) and then purified using gel extraction (QIAquick Gel Extraction Kit, Qiagen). The library of oligonucleotides was cloned into 400 μg of linearized pJD001 substituting the wildtype HSX1 gene using a Gibson reac tion (prepared in house) at 50 °C for 12 h with a ratio 5:1 of insert:vector. After dialysing the reaction with 0.025 μm VSWP membrane filters (Merck Millipore) for 1.5 h, the product was concentrated 4× using speedvac. Six microlitres of the concentrated reaction was transformed into 100 μl of electrocompetent Escherichia coli (NEB 10beta Electrocompetent E. coli, NEB) according to the manufacturer's protocol. Cells were allowed to recover in SOC (NEB 10beta/ Stable Outgrowth Medium) for 30 min and later transferred to 150 ml of LB medium with ampicillin 4× overnight. The total number of transformants was estimated to be ~9.59 × 10 6 . Given the complexity of the library, each variant was therefore represented ~1,849 times on average. 50 ml of E. coli saturated culture was harvested to extract the plasmid library using plasmid midi prep (QIAfilter Plasmid Midi Kit, Qiagen). Selection experiment. Yeast strain and conditional growth defect in different environmental conditions. The HSX1 deletion strain was obtained by replacing the HSX1 gene with a nourseothricin resistance cassette in the haploid laboratory strain BY4742 (MATα his3Δ1 leu2Δ0 lys2Δ0 ura3Δ0 HSX1::natMX4) and later confirmed using colony PCR. The deletion of the single copy tRNAArg(CCU) (HSX1) in yeast was previously reported to lead to a conditional growth defect when the temperature is raised from 30 °C to 37 °C 15 . We found that a similar growth defect is observed if the growth medium contains high salt concentrations (1 M NaCl), and that a combination of high temperature and high salt gives an even stronger defect (Extended Data Fig. 1a). Synthetic complete medium lacking histidine (SCHIS) 1 M NaCl at 37 °C was therefore used as the selective condition for the library selection experiment. Large-scale yeast transformation. The highefficiency yeasttransformation pro tocol was derived from a previously described method 7 . Two precultures of the tRNA deletion strain were grown independently in 25 ml standard YPDA at 30 °C overnight. The next morning, the cultures were diluted into 175 ml of fresh YPDA to OD 600 nm = 0.3. The two cultures were incubated at 30 °C for 4 h (~2-3 gener ations). After the growth period, the cells were harvested and centrifuged for 5 min at 3,000g, washed in sterile water and later in SORB (100 mM LiOAc, 10 mM Tris pH 8.0, 1 mM EDTA, 1 M sorbitol). The cells were resuspended in 8.6 ml of SORB and incubated at room temperature for 30 min. After incubation, 175 μl of 10 mg ml −1 boiled salmon sperm DNA (Agilent Genomics) was added to each tube of cells, as well as 3.5 μg of plasmid library. After 10 min of gentle shaking at room temperature, 35 ml of Plate Mixture (100 mM LiOAc, 10 mM TrisHCl pH 8, 1 mM EDTA/NaOH, pH 8, 40% PEG3350) was added to the cells and incubated at room temperature for 30 more min. 3.5 ml of DMSO was added to each tube and the cells were then heat shocked at 42 °C for 20 min (inverting tubes from time to time to ensure homogenous heat transfer). After heat shock, each independent tube of cells was centrifuged and resuspended in 350 ml of YPD + 0.5M Sorbitol and allowed to recover for 1 h at 30 °C. The cells were then centrifuged, washed twice with SCHIS medium and resuspended in 350 ml SCHIS. The two independent transforma tions were grown at 30 °C for ~60 h until saturation. For the two independent trans formations, 1.5 × 10 6 and 1.1 × 10 6 transformants were obtained, which ensured that each variant of the library was on average represented ~250 times 34 . Competition assay. The competition experiment had two different phases. In phase one, the environment had minimal selection on the tRNA functionality (SCHIS at 30 °C), allowing the pool of variants to be amplified and the cells to enter the exponential growth phase (input library) 34 . In the second stage, the medium was changed to a condition (SCHIS 1 M NaCl medium at 37 °C) in which non functional tRNA variants would lead to a severe growth defect phenotype (output library). The assay was performed immediately after yeast transformation to avoid recovering cells from frozen glycerol stocks. Once the two independently trans formed cultures reached saturation (~60 h after plasmid transformation), they were inoculated at an OD 600 nm of 0.08 in 500 ml of SCHIS medium and grown for four generations at 30 °C (~11 h). When exponential phase was reached after four generations of growth, the cells were harvested and washed with selection medium (warm SCHIS NaCl 1 M) and then inoculated in 500 ml of selection medium at an OD 600 nm of 0.015. The remainder of the cells was harvested and stored at −20 °C for later DNA extraction of the input libraries. Each independent input library was divided into three different output libraries (six replicates in total). Cells were grown in selective conditions for ~6.5 generations (~26.5 h). This number of generations was chosen so that the average read coverage in the input would be of ~150 reads per variant and that null alleles, which grow ~0.18 generations every 3 h, would be detected in the output after sequencing. After 6.5 generations, the cells were harvested and the cell pellets stored at −20 °C for later DNA extraction of the output libraries. DNA extraction and quantification. Cell pellets (eight tubes, two inputs and six outputs) were resuspended in 1.5 ml extraction buffer (2% TritonX, 1% SDS, 100 mM NaCl, 10 mM TrisHCl pH 8, 1 mM EDTA pH 8), frozen using an dry ice-ethanol bath and incubated at 62 °C in a water bath twice. Subsequently, 1.5 ml of phenolchloro isoamyl alcohol (25:24:1 ratio, equilibrated in 10 mM TrisHCl, 1 mM EDTA, pH8) was added, together with 1.5 g of glass beads and the samples were vortexed for 10 min. Samples were centrifuged at room temperature for 30 min at 3,200g and the aqueous phase was transferred into new tubes. The same step was repeated twice. 0.15 ml of NaOAc 3 M and 3.3 ml of cold ethanol 100% were added to the aqueous phase. The mix was incubated at −20 °C for 30 min and then centrifuged for 30 min at full speed at 4 °C to precipitate the DNA. The ethanol was removed and the DNA pellet allowed to dry overnight at room tem perature. DNA pellets were resuspended in 900 μl TE 1× and treated with RNaseA (10 mg ml −1 , Thermo Scientific) for 30 min at 37 °C. To desalt and concentrate the DNA solutions, a QIAEX II Gel Extraction Kit was used (75 μl of QIAEX II beads suspension). The samples were washed three times with PE buffer and eluted twice in 375 μl of 10 mM Tris·Cl buffer, pH 8.5. Sequencing library preparation. The plasmid concentration in each total DNA sample was quantified in triplicate by realtime quantitative PCR, using prim ers that had homology to the origin of replication region of the pJD001 plasmid backbone (Supplementary Table 3). On average, we obtained ~3.5 × 10 6 plasmid molecules per μl of DNA sample.
A twostep PCR using high fidelity Q5 Hot Start HighFidelity DNA Polymerase (NEB) was used to amplify the input and output libraries for sequencing. For each sample, ~150 million plasmid molecules were amplified for ten cycles using prim ers with overhang homology to Illumina sequencing adapters (Supplementary Table 3). The samples were then treated with ExoSAP (Affymetrix) and cleaned using bead purification with a QIAEX II kit (10 μl pf QIAEX II beads suspen sion). The whole eluates, corresponding to the entire first PCR reactions, were used for the second PCR reactions (15 cycles), in which the rest of the Illumina adaptor was added as overhangs on the primers, in addition to samplespecific indexes. The DNA concentration of each individual second PCR was quantified by fluorometric quantitation (QuantiT PicoGreen dsDNA Assay Kit) and pooled together at an equimolar ratio. Finally, the pooled sequencing library was gel puri fied (QIAEX II Gel Extraction Kit) and subjected to 125 bp pairedend sequencing on an Illumina HiSeq 2500v5 sequencer at the EMBL Genomics Core Facility (Heidelberg, Germany). From sequencing reads to fitness values. The sequencing reads of each sam ple (two inputs and six outputs) were processed and filtered independently. Each sequencing read covered the entire tRNA. The 5′ and 3′ constant regions of the read (primers annealing sites) were removed with the 'cutadapt' software 35 . The forward and reverse reads were merged using PEAR 36 and sequences that were either not assembled owing to low quality or unexpected length were discarded. Unique gen otypes were called and quantified with custom Python scripts. Genotypes with less than nine input reads in any input replicate, unexpected nucleotide substitutions (sequencing or PCR errors) or zero reads in the outputs were discarded. After filtering, we ended up with a total of 4,176 sequence genotypes quantified in all inputs and outputs.
To obtain accurate fitness and error estimates for each variant we took into account the hierarchical structure of the replicates 37 as well as sampling error owing to the low number of read counts 38 . Input and output frequencies for each genotype from each of the independent competition experiments were first calculated and then these were combined into a single output measurement for each input repli cate. in which g is the genotype (from 1 to l, with l being the total number of genotypes after filtering), i is the number of input replicates (1 or 2) and j is the number of output replicates per input replicate (1 to 3). These formulae assume that each read derives from an individual cell, so that by multiplying the frequency of reads in the output with the final (OD out ) and initial culture density (OD in ) we can estimate the number of cells for a particular genotype at the beginning (f in ) and end (f out ) of the competition experiment.
Each input and output frequency is associated to a Poisson variance given the number of read counts of each genotype and the total read count 38 : We calculated a single output frequency score for each input replicate using a weighted average in which the weight of each score ( f out gij ) is the inverse of the variance of the genotype (σ out The output frequency errors of each replicate were then combined to yield an overall output frequency error: The number of generations n gi was then calculated as the log2 ratio of the nor malized input and output frequencies: with an associated error of: The number of generations in each input replicate (n g1 and n g2 ) was combined using a weighted average as before to obtain a single growth measurement and an error for each genotype: In logspace, if a particular genotype grew faster or slower than the wild type, the ln(fitness) value would be >0 or <0, respectively. Single mutation effects, pairwise genetic interactions and higher order epistasis. On a logscale, the fitness effect of a mutation A on a genetic background X was calculated as the relative fitness gain of the variant AX respect to X: This fitness effect of a mutation can also be referred to as the first order epistatic term (ε 1 ) 39 .
A pairwise epistatic interaction between two mutations was defined as the dif ference between the observed fitness of the double mutant AB and the expected fitness obtained by the addition of the two single mutant fitness values (A and B). The fitness effects of the mutations A, B and AB can be calculated on each genetic background X by subtracting the fitness of X itself from the fitness of AX, BX and ABX, as described above. Pairwise epistasis (or secondorder epistasis ε 2 ) is then the change in the effect of each single mutation in the presence of the second mutation: This same analysis can be expanded to higher order terms 22,39 . For example, a thirdorder interaction (ε 3 ) is the degree to which secondorder epistasis is differ ent when a third mutation is present in the background: Higher order terms follow the same principle, so we can calculate any n th order term using the formula 39 : ( 1) (1) in which ω n are all fitness terms of order n in a specific genetic background. It is important to note that an epistatic term of any order n can only be calculated if the genotype space is complete (that is, that the fitness of all genotypes from order 0 to n were quantified in the experiment). In our dataset, higherorder epistasis was quantified up to order eight (76 cases in this dataset), which was the highest order in which the fitness of a combinatoriallycomplete set of genotypes could be quantified after data filtering (Extended Data Fig. 1d).
To quantify how many epistatic terms were significantly positive or negative across all the backgrounds in which they were tested, a onesample ttest was per formed (using the epistatic term and its respective propagated error). The FDR was adjusted across all the tests performed (a total of 203,240 tests for all interactions of all orders across all backgrounds) using the Benjamini-Hochberg method 40 . Controlling for background fitness, sequence divergence and the number of input sequencing reads. Across all the data, there was a weak correlation between the fitness of the genetic background and both the fitness effect of the single muta tions and pairwise epistasis (Extended Data Figs. 2c, 5a). We therefore repeated all of the analyses on the subset of the genetic backgrounds with fitness close to the wildtype S. cerevisiae (−0.15 < fitness < 0.15, n = 1,479 library genotypes) and also on genetic backgrounds with moderate fitness decreases (−0.3 < fitness < −0.15, n = 1,577). We also repeated all of the analyses on the genetic backgrounds that were closest to the S. cerevisiae sequence (one to four mutations away, n = 1,040) or excluding all variants with a mean input frequency of less than 100 reads (n = 1,315). With each of these filters we excluded approximately two thirds of the original number of variants in the library. Classifying pairwise epistasis. Significant pairwise interactions in the dataset (n = 10,330 out of 47,649 tested) were classified into three categories: magnitude, sign, and reciprocal sign epistasis 41 . Pairwise epistasis was thus classified as follows. When the fitness effect of both single mutants differs in magnitude but not in sign in the presence of the other mutation, the epistatic interaction was classified as Letter reSeArCH magnitude epistasis. For sign epistasis, the sign of one of the individual fitness has effects on changes in the presence of a second mutation. Finally, if the sign of effect changes for both individual mutations, the interaction was classified as reciprocal sign epistasis. The way a single mutation effect changes in the presence of another mutation can be inferred if the fitness effect and sign of the single mutations (A and B) and the fitness of the double mutant (AB) are known. For instance, if the two single mutations A and B have significantly beneficial (positive) effects and the double mutant has higher fitness than both single mutants, then none of the single mutations are changing sign, so this interaction would be classified as magnitude. However, if the double mutant has a fitness value lower than both single mutations, then this interaction would be classified as reciprocal sign (both single mutations are changing sign in the presence of the other). Otherwise, this interaction will be classified as sign (fitness of the double is lower than only one of the singles).
The sign of each of the single mutants in the dataset (n = 21,450) was assigned after performing a onesample ttest (Benjamini-Hochberg FDR controlled across all tested interactions of all orders from one to eight, n = 203,240 as described in the Methods section 'Single mutation effects, pairwise genetic interactions and higher order epistasis'). Single mutants with q ≥ 0.1 were assigned as neutral (or notsignificant) and the rest as positive (beneficial) or negative (deleterious) when the fitness effect of the mutation was more or less than 0 respectively. Exceptional interactions between two mutations in which both single mutations had a neutral category (no significant fitness effect at FDR < 0.1) were classified as magnitude epistasis (either positive or negative). When only one of the single mutations had a neutral category they were then classified as sign or magnitude epistasis depending on whether the other single mutation changed sign or not. Whenever both single mutations had either positive or negative categories, epistasis was classified as explained above. Background-averaged epistatic interactions. We quantified the backgroundav eraged epistatic interaction of a particular mutation combination (ranging from order one to eight) by averaging all epistatic coefficients of that same combination of mutations across all backgrounds in which it was found. To assess the signifi cance of the average epistatic coefficient, the errors of all individual fitness terms were propagated and a onesample ttest was performed. The P value was adjusted for all tests performed from order one to eight (a total of 3,691 tests) using the Benjamini-Hochberg FDR method 40 .
After identifying those mutations that interacted significantly when averaging across backgrounds (at FDR < 0.1), we counted the number of times the interac tions between two mutations changed owing to another mutation in the back ground, or calculated the number of times a single mutation was able to change a pairwise interaction (Fig. 4b, c). Genetic prediction. As described in the section 'Single mutation effects, pairwise genetic interactions and higher order epistasis' , epistatic terms were calculated as linear combinations of the fitness values of genotypes of different orders. This system of linear combination can be represented in a matrix form, which allows the epistatic coefficients to be calculated from fitness values, and fitness values back from epistasis 39 .
In a complete n loci diallelic genotype space, in which each locus can harbour two different nucleotides, epistatic terms can be calculated as follows: in which ω corresponds to a vector with the fitness values of the 2 n genotypes from order 0 to n, ε is a vector with all the corresponding epistatic terms and G is a matrix that defines the linear mapping between ω and ε for all orders. G can be recursively constructed as follows: In this case, epistatic terms are calculated relative to a single background (0 th order genotype or 'wt'). However, within a complete landscape, epistatic terms can be calculated across many different backgrounds. For instance, in a diallelic landscape of three loci, the same single mutation effect (epistasis term of order one) can be measured four times from four different backgrounds. To obtain epistatic coefficients averaged amongst backgrounds we can use a similar version of the previous equation:= Fitness values can be obtained by a linear combination of epistatic coefficients using the inverse mapping, for either relative or backgroundaveraged epistatic coefficients:̄ω For an overview and extended definitions, we refer the reader to the previously published description 39 . Cross-validation. To detect model overfitting, we used a tenfold crossvalidation approach in which the backgroundaveraged epistatic coefficients were quantified using 90% of the genotypes (training set) within each of the 76 eightloci tRNA sublandscapes with the remaining 10% used for evaluation (test set). With 10% of genotypes missing, computation of seventh or eighth order coefficients is no longer possible. Coefficients of other orders were averaged across backgrounds for which all intermediate genotypes were available. To assess the significance of each epistatic coefficient, the estimates of fitness errors were propagated accordingly and the tstatistic for a one sample ttest was calculated. Within each of the ten training sets for each complete sublandscape, the coefficients were ranked by their absolute tstatistic and cumulatively used to predict fitness of the heldout test set genotypes (least significant coefficients were iteratively set to zero before predicting fitness values) using the inverse of the Walsh-Hadamard transform as described above (using a weighting matrix V in which the weights correspond to the number of backgrounds each coefficients had been averaged across). The best predictive model for each of the ten training sets of each sublandscape was selected as the model that gave the lowest prediction error on the corresponding test set (Extended Data Fig. 8).
The accuracy of all the above predictions was quantified using root mean square error (RMSE): = n RMSE SS res in which SS res is the residual sum of squares and n is the total number of predicted genotypes. To calculate the percentage of variance explained (PVE) we used the formula: res total in which SS total is the total sum of squares. Comparisons to theoretical fitness landscapes. We used three different landscape statistics (γ statistic 28 , roughnesstoslope ratio 42 and the proportion of epistasis types 42 ) to compare the tRNA fitness landscape to theoretical landscapes. To esti mate the robustness of these measurements, all the statistics were calculated for all possible diallelic (two possible nucleotide substitutions per position) complete tRNA sublandscapes from three to eight loci that started from the wildtype S. Letter reSeArCH was modelled with a mix of 50% additive and 50% HoC and the K parameter of the NK model (each locus interacts with K loci) was set to K = n/2. These parameters were selected as they resulted in landscape statistics most similar to those of the tRNA sublandscapes (data not shown). γ statistic: correlation of fitness effects. The γ statistic was recently introduced 28 and extended by others 13 . γ quantifies the correlation of fitness effects of the same mutation in singlemutant neighbours. It measures how the effect of a focal muta tion is altered by another mutation at another locus in the background, averaged across the whole landscape. The statistic is bounded between −1 and 1. In a sce nario without epistasis (the effect of a mutation is completely independent of the background), γ = 1. The γ measure gives information on the amount of epistasis in a combinatoriallycomplete landscape, but does not discriminate between different landscape topographies (two landscapes that differ in structure can have the same γ value). As with γ, γ d (the decay of correlation of fitness effects with mutational distance) can be defined as the correlation of fitness effects of mutations between genotypes that are 1, 2, 3…, d mutations away. γ d gives extra information about the structure of the landscape, as it describes the cumulative epistatic effect of d mutations 13,28 . In a completely additive landscape, γ d is always 1 because the effect of a mutation is independent of the background genotype that is 1, 2, 3 or d muta tions away. However, in a maximally rugged fitness landscape (in which the effect of a mutation depends entirely on its genetic background) γ 1 is 0 and γ d is 0 for all values of d. The behaviour of γ d as a function of d varies for different theoretical landscape models 13,28 (Extended Data Fig. 9a). We calculated γ d values for all possible complete diallelic tRNA sublandscapes of three to eight mutations combinations that contained the S. cerevisiae genotype using the software MAGELLAN (eight being the maximum number of loci in which a complete genotype space is available in the dataset). We later compared the statistic to the values for the theoretical landscapes. As a measure of similarity, we calculated the Euclidean distance between the γ d of all tRNA sublandscapes and the γ d of the theoretical models (each tRNA landscape, n = 73,250, 142,000, 159,500, 100,750, 33,000 and 4,500 for tRNA landscapes from three to eight mutations respectively, was compared to the 250 simulations of each theoretical landscape). Other quantitative measures of landscape ruggedness. In addition to the γ statistic, for all complete tRNA and theoretical sublandscapes from three to eight loci, we also calculated the roughnesstoslope ratio (r/s ratio) and characterized the local pairwise epistatic interactions. The r/s ratio measures how well the landscape can be described by a linear model, which corresponds to the purely additive limit 42 . The roughness is given by the variance of the residuals from the linear model and the slope is given by the average of the absolute values of the linear coefficients. The higher the r/s, the higher the deviation from the linear model and the more epistasis is present (in a nonepistatic scenario, r/s = 0). To characterize the local interactions of each landscape we calculated the fraction of magnitude, sign or reciprocal sign pairwise epistasis within each landscape. We used the software MAGELLAN to calculate all the described statistics. Accessible paths between extant species. An accessible path between two geno types in the landscape was defined as a mutation trajectory in which none of the intermediate genotypes has significantly lower fitness than both the initial and final genotypes that they connect (ttest between all the intermediate genotypes against the origin and endpoint genotypes, n = 1-8 tests). A path that had at least one deleterious intermediate genotype (P < 0.05) was classified as inaccessible. We measured the number of accessible direct (shortest) paths between 20 pairwise comparisons of the extant genotypes in the landscape using the R package igraph. Statistical analyses. All statistical analyses were performed in R (v.3.3.3) and figures were made using the R package ggplot2. Lower and upper hinges of box plots correspond to the first and third quartiles (25 th and 75 th percentiles). The upper and lower whiskers extend from the hinge to the largest and lower value no further than 1.5 × IQR (interquartile range) respectively. Higher or lower points (outliers) are plotted individually (or not plotted in those cases were the box plot is plotted together with a violin plot). Notches give roughly 95% confidence interval for comparing the medians. Reporting summary. Further information on experimental design is available in the Nature Research Reporting Summary linked to this paper. Data availability. The complete dataset is available as Supplementary Table 1. Custom code used in this study is available from the authors upon request. Raw sequencing data has been submitted to GEO (accession number GSE99418).  27 ; sequences with indels were excluded). d, Number of positive, negative or not significant pairwise interactions at FDR < 0.1 within the acceptor stem of the tRNA (n = 23,237) when both mutations are found in the same helix strand or when each mutation is located in a different strand (n = 13,615). log2 odds ratio shown below together with twosided Fisher's exact test P values. e, Number of positive, negative and nonsignificant backgroundaveraged pairwise interactions between pairs of mutations in the acceptor stem that are found in the same RNA strand and between mutations that are in positions that base pair with each other. log2 odds ratio and twosided Fisher's exact test P values are shown below. f, Distribution of pairwise epistasis values of mutation pairs that restore a canonical WCBP depending on the location of their background mutations in the acceptor stem (P values from Welch's twosided ttest, n = 263 or n = 1,368 when more than one background mutations are in the same strand or not, respectively). The same result is obtained when epistasis values are corrected for the ln(fitness) of the background (residuals of a linear model using background ln(fitness) to predict epistasis, data not shown).

Letter reSeArCH
Extended Data Fig. 7 | Background-averaged third and higher-order interactions. a, The most significant backgroundaveraged thirdorder interactions (8 out of 74, FDR < 0.1, n = 3,691 tests for all interactions across all orders). The first three plots of each row show how the distribution of pairwise epistasis of two mutations across different genetic backgrounds (each double mutation can be found in a median of 506 different genetic backgrounds) changes in the presence or absence of a third mutation. The paired differences between pairwise interactions in those three cases correspond to third order epistatic coefficients. Distributions of thirdorder epistasis for the same three mutations are shown to the right. Horizontal lines correspond to the background averaged thirdorder epistatic term, coloured by sign (orange or green for positive or negative respectively). b, Number of significantly positive and negative backgroundaveraged epistatic interactions of order one to eight (at FDR < 0.1). c, Distribution of the absolute magnitude of averaged thirdorder interactions plotted against the mean nucleotide distance between the three mutations (n = 316 triple mutations). Welch's twosided ttest P values for differences between the groups are shown. Significant interactions (onesample twosided ttest at FDR < 0.1) are coloured in orange or green for positive or negative epistasis respectively. d, Top, Number of positive, negative or nonsignificant backgroundaveraged thirdorder interactions (FDR < 0.1) within the acceptor stem of the tRNA when both mutations are found in the same helix strand or not (n = 129). Bottom, the log2 odds ratios (when all three mutations are found in the same strand of the tRNA acceptor stem) of significantly positive interactions versus others (negative or not significant interactions) and significantly negative interactions versus other double mutants. P values reported from the twosided Fisher's exact test.
1 nature research | reporting summary

Statistical parameters
When statistical analyses are reported, confirm that the following items are present in the relevant location (e.g. figure legend, table legend, main text, or Methods section).

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement An indication of whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistics including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Software and code
Policy information about availability of computer code Data collection All the data analysed in this study was produced in this study. The sequencing reads of each sample (two inputs and six outputs) were processed and filtered independently. Each sequencing read covered the entire tRNA. The 5' and 3' constant regions of the read (primers annealing sites) were removed with the 'cutadapt' software. The forward and reverse reads were merged using 'PEAR' and sequences that were either not assembled due to low quality or unexpected length were discarded. Unique genotypes were called and quantified with custom python scripts. Genotypes with less than nine input reads in any input replicate, unexpected nucleotide substitutions (sequencing or PCR errors) or 0 reads in the outputs were discarded.

Data analysis
All data analysis were performed in R (version 3.3.3). We used the software 'MAGELLAN' to generate theoretical fitness landscapes and calculate the gamma statistic to compare the tRNA landscape to theoretical models.
For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information. Data