Continuous evolution of user-defined genes at 1 million times the genomic mutation rate

Abstract

When nature evolves a gene over eons at scale, it produces a diversity of homologous sequences with patterns of conservation and change that contain rich structural, functional, and historical information about the gene. However, natural gene diversity accumulates slowly and likely excludes large regions of functional sequence space, limiting the information that is encoded and extractable. We introduce upgraded orthogonal DNA replication (OrthoRep) systems that radically accelerate the evolution of chosen genes under selection in yeast. When applied to a maladapted biosynthetic enzyme, we obtained collections of extensively diverged sequences with patterns that revealed structural and environmental constraints shaping the enzyme's activity. Our upgraded OrthoRep systems should support the discovery of factors influencing gene evolution, uncover previously unknown regions of fitness landscapes, and find broad applications in biomolecular engineering.

Fig. 1.. Engineering orthogonal DNA polymerases for increased mutation rates.

(A) Architecture of the OrthoRep system. A DNA polymerase (TP-DNAP1) that exclusively replicates a specific cytoplasmically localized plasmid via protein primed replication at a high error rate enables in vivo targeted mutagenesis without mutagenizing genomic DNA. (B) Schematic for a directed evolution approach to engineer TP-DNAP1’s mutation rates and mutation spectrum incorporating both a direct selection for rare transversion mutations as well as high accuracy mutation rate measurement using a mutation accumulation and high throughput sequencing (HTS) assay. (C) Mutations identified in TP-DNAP1 variants presented in this study. Color and single letter amino acid code are used for only the first TP-DNAP1 in which a mutation is identified. Color only is shown for all other instances of a mutation. (D-F) Mutation rate measurements via mutation accumulation for a series of TP-DNAP1 directed evolution intermediates showing either overall mutation rates as boxplots (D), mean mutation rates for individual substitution types as heatmaps (E), or mutation rates for individual substitution types for three individual TP-DNAP1 variants as boxplots (F). Points are representative of individual biological replicates, each representing 3-4 timepoints with >50 sequences each. Boxplots and heatmaps are representative of n=4 biological replicates. Box plot central line, boxes, and whiskers represent the median, interquartile range, and minimum / maximum values, respectively.

Figure 2.. Massively parallel continuous diversification and evolution of TrpB.

(A) Schematic for OrthoRep-driven evolution of the tryptophan synthase β-subunit from Thermotoga maritima (TrpB) for standalone function in yeast. TrpB was integrated onto the p1 plasmid in a yeast strain lacking the native yeast tryptophan synthase gene (TRP5). 96 independent cultures of the resulting strain were passaged mostly under selective pressure for Trp production using exogenously supplied indole over ~540 generations. DNA from fifteen timepoints throughout the evolution campaign was harvested and sequenced using HTS TrpB illustration generated using Illustrate (59). (B) Selection pressure for TrpB function is applied by lowering or eliminating exogenously supplied Trp and lowering exogenously supplied indole over time. The schedule of selection pressure imposed throughout extensive evolution is plotted. Timepoints at which cultures were harvested and sequenced using HTS are indicated. Selection periods were characterized as ‘no selection’, ‘mostly positive selection’, and ‘mostly purifying selection’ based on the Trp and indole amounts supplied and the progress of evolution. See Methods and table S2 for a description of media conditions and selection pressure derivation. (C-E) Distribution of amino acid Hamming distances for both pairwise comparisons and comparisons to the wild-type sequence, at the first and last timepoint (D and E respectively) and as pairwise cumulative distributions for all timepoints (C).

Figure 3.. Revealed effects of selective constraints.

(A) AlphaFold structure of Thermotoga maritima TrpB with different regions colored according to their structural role. First shell, β-β interface, and β-α interface residues are designated as such if they are within 5 Å of the substrate and cofactor (Trp and PLP), the other β subunit, or the other α subunit in the αββα heterotetramer holoenzyme, respectively. Alignment to Pyrococcus furiosus TrpB crystal structures (PDB codes 5E0K and 5DW3) were used to determine distances from substrate and cofactor, α-subunit, and β-subunit. Mean solvent accessible surface area (SASA) was used to categorize all remaining residues as either surface (SASA≥0.2) or buried (SASA<0.2). (B) Heatmap of mutations among OrthoRep-evolved TrpB sequences applied to the AlphaFold structure. (C) Distributions of mutations among OrthoRep-evolved TrpB sequences within the 6 structural regions compared to a null model composed of a simulated dataset of sequences mutagenized in silico according to the mutation rates and preferences measured for TP-DNAP1 BadBoy2 until the number of synonymous mutations in the simulated sequences and real sequences were equivalent. (D) Heatmap of mutation frequency for all mutations among the 20 most frequently mutated positions in the timepoint corresponding to either 50 or 70 generations. Frequencies for simulated sequences were subtracted to account for bias due to BadBoy2’s mutation preference and wt TrpB sequence content. (E) Violin plots of isoelectric points for all OrthoRep-evolved and simulated TrpB sequences, split by timepoint. Points and black bars denote the means and interquartile range for all sequences within each timepoint. Isoelectric points of the wild type T. maritima TrpB, an N-terminally-truncated Trp5 homologous to TrpB (Trp5-ΔN), the native T. maritima TrpA-TrpB complex, and the TrpA-TrpB holoenzyme ortholog from Saccharomyces cerevisiae (Trp5) are shown for comparison. Note that the TrpB sequence used as a starting point for evolution includes a 6xHis tag that contributes an increase in pI of 0.2, and that the native TrpA-TrpB complex pI shown for reference does not include this 6xHis tag.

Figure 4.. Lineage barcodes reveal covarying residues.

(A) Schematic of computational processing used to reduce phylogenetic contingency of residue covariation. (B) Residue covariation plot for the most frequently covarying residues among all timepoints in the TrpB evolution dataset for 93 lineages downsampled to 100 sequences per lineage. (C-D) Heatmaps of most frequently mutated residues among sequences containing mutations A20V or A20T, downsampled to 100 sequences and chosen from specific timepoints. Heatmaps of all positions for early (generations 70 and 90) and late (generations 480 and 540) timepoints are overlaid onto a TrpB AlphaFold structure with the 20 most frequently mutated positions shown as spheres (C) or shown as mutation-specific heatmaps of the most frequently mutated 20 positions in late timepoints (D).

Figure 5.. Pooled measurement and TransceptEVE prediction of TrpB variant fitness.

(A) Schematic of pooled TrpB fitness assay using HTS. (B) Spot plating growth assay of control sequences included in the pooled fitness assay. (C) Hexbin plot of replicate concordance among pairs of replicates under growth conditions with Trp (no selection), without Trp and with 400 uM indole (weak selection) or without Trp and with 25 uM indole (strong selection) for highly functional sequences (enrichment score > −5) (D) Distribution of mean enrichment scores (average of n=2 biological replicates) among the three selection conditions. (E and F) Hexbin plot of TranceptEVE score vs. measured mean enrichment score with strong selection for either all enrichment scores (E) or enrichment scores for highly active sequences (F). The percentage of all sequences that fall in the upper or lower quartile of score predictions and are classified as either high or low function (enrichment score greater than or less than −5, respectively) are shown in the respective sections of the plot in (E). (G) Hexbin plot of TransceptEVE score vs. number of nonsynonymous substitutions for all sequences with a measured strong selection mean enrichment score. r, Pearson correlation.