Tuning evolvability via plasmid copy number and regulatory architecture

Domain: biology
Published: 2025-12-31
Editorial outcome: minor_revision

Abstract
Genetic modules are often designed and implemented with inspiration from engineering disciplines. Although this approach can be successful because of the similarities underpinning physical and biochemical systems, it neglects a key factor that affects the performance of living organisms: evolution. Thus, it is crucial to incorporate the impact of inevitable mutations into the design and analysis of genetic modules. Combining computational modeling and in vivo mutagenesis experiments in Escherichia coli, we characterize how the interplay of gene dosage via plasmid copy number (PCN) and regulatory architecture affect the phenotypic mutation rate. For example, while greater PCN facilitates the emergence of gain-of-function mutations, it instead curbs the spread of loss-of-function mutations. We further reveal that mutations in the coding region are often masked at the phenotypic level, unlike those occurring in the regulatory region which become more prominent as PCN increases, both when the regulator is expressed constitutively and when it is self-repressed. Together, our results shed light on evolutionary organizing principles and aid the rational design of both evolutionarily stable and highly evolvable biocircuits.

Introduction
The design of synthetic gene circuits is often inspired by engineering principles1,2. For instance, modules that realize sophisticated analog3 and digital functions4 can be implemented by borrowing concepts from electrical engineering5–10, and complex systems can be realized leveraging computer-aided automated design tools11,12 in a wide array of contexts ranging from biomanufacturing to living therapeutics13–18.

Unfortunately, the analogy is imperfect, as modules display context-dependent behavior due to a variety of factors, such as competition for shared cellular resources19,20 and the physical layout of circuits21. Representing an even more fundamental difference, while electrical circuits are static elements, genetic components are instead subject to constant changes due to mutations and evolutionary pressure, eventually leading to design failures2. Therefore, it is imperative to integrate evolutionary considerations into the design and analysis of synthetic gene circuits22–25 to fine-tune the balance between (i) evolvability to optimize or innovate novel traits26–30, and (ii) evolutionary stability to preserve the phenotype in the presence of genetic mutations31–34.

Synthetic gene circuits are typically deployed on plasmids, offering a convenient platform to modulate their gene dosage via their copy number. These self-replicating extrachromosomal DNA vectors are present in a wide range of bacteria, and provide an avenue for accelerated evolution35. Plasmids are most often shared among bacteria via horizontal gene transfer36 to aid the spread of beneficial evolutionary traits, such as resistance to antibiotics, degradation of toxic chemicals, the ability to synthesize essential amino acids and to undergo anoxygenic photosynthesis37–40. Unlike the chromosome, plasmids are generally present at multiple copies in each host, further accelerating evolution through gene amplification41: upon mutation, high copy plasmids are likely to be heterozygous, yielding enhanced genetic diversity and plasticity to overcome evolutionary trade-offs42.

Numerous in vivo mutation systems have been developed in recent years to improve directed evolution by providing continuous genetic diversification29,43–46 to closely mimic the dynamics in population genetics, where mutations occurs stochastically and across generations. Compared with in vitro mutations, in vivo evolution yields (i) mutations that are directly linked to host fitness; (ii) parallel streams of mutations due to the constantly replicating population; and (iii) potentially wider and longer mutation paths43. Thus, instead of utilizing pre-diversified genetic variants, we adopt an in vivo system to generate mutations using EvolvR: an enhanced Cas9 nickase fused with a modified error-prone DNA polymerase which is directed by a single guide RNA (sgRNA) to generate mutations at a constant rate at user-defined locations46.

Harnessing EvolvR, we implemented a collection of engineered circuits to study the emergence of plasmid-borne mutations. We show that the interplay between gene dosage and regulatory motifs shapes both genotypic and phenotypic heterogeneity, highlighting the value of integrating synthetic gene circuits and in vivo mutagenesis to study evolutionary dynamics. By developing a computational model to simulate mutation, plasmid replication and partitioning over subsequent generations, we show that the number of cells with at least one mutant plasmid increases with PCN, which we then verify experimentally. Following this, we demonstrate that while for gain-of-function mutations the rate at which phenotypic mutations occur depends on the detection limit, for loss-of-function mutations this rate strictly decreases with PCN. Our work further reveals that mutations in the coding region tend to be masked at the phenotypic level, unlike those in the regulatory region which become more prominent as PCN increases. Our results demonstrate that circuit-level features such as gene dosage and regulatory motifs that shape gene expression impact the emergence and spread of mutants at the phenotypic level despite identical mutation dynamics at the DNA level. Thus, our work provides explicit design guidelines to modulate the dominance of mutants within a population, and to incorporate evolutionary considerations into the design of complex synthetic gene circuits.

Results
Phenotypic mutation rate is expected to depend on gene dosage
To investigate how the phenotypic mutation rate is impacted by gene dosage, we adopted a modified version of the standard Wright-Fisher model47,48 to perform stochastic simulations. This framework is commonly used in population genetics, although typically for haploid and diploid organisms. Thus, we needed to slightly modify it to account for the fact that PCN generally exceeds two, corresponding to a multiploid population in the Wright-Fisher model.

Our computational framework (Fig. 1a) has discrete and non-overlapping generations with a constant population size N, and each member is replaced in every generation as follows. Let x_wt and x_m denote the number of wild-type and mutated plasmids in a given cell, respectively, where P is the PCN. Similarly, let N_wt and N_m = N - N_wt denote the number of wild-type and mutant (with at least one mutated plasmid) bacteria, respectively. We generate the number of new mutations Δx within each cell as follows: each wild-type plasmid becomes mutated with probability μ, modeled as a Poisson process. Next, the numbers of wild-type and mutated plasmids for each cell are updated (STEP 1 in Fig. 1a) as x'_wt = x_wt - Δx and x'_m = x_m + Δx, respectively. Subsequently, the number of new plasmids P' after replication is generated according to a second Poisson distribution, with its mean equal to the average PCN (see Methods), and this pool is split according to a binomial distribution with parameters P' and x'_m / (x'_m + x'_wt), that is, newly replicated plasmids are randomly assigned to be either wild-type or mutated according to the respective proportions within the cell (STEP 2 in Fig. 1a). Following this, both daughter cells receive either (P + P')/2 plasmids if P + P' is even, otherwise one of them ends up with one more than the other to ensure that there is no segregational loss, where a hypergeometric distribution49 determines how mutated and wild-type plasmids are distributed (STEP 3 in Fig. 1a). Finally, half of the population is randomly selected to initialize the next generation (STEP 4 in Fig. 1a), while the rest is discarded to keep the population size constant across generations.

Relying on parameter values typical in E. coli (see Methods), data in Fig. 1 reveal the following. First, the size of the mutant subpopulation increases with PCN (Fig. 1b), and while the average number of mutated plasmids also increases with PCN within this subpopulation (Fig. 1c), their fraction instead decreases (Fig. 1d). Greater gene dosage thus results in more bacteria with at least one mutated plasmid, and they will harbor more mutated variants in absolute numbers, but not in relative terms. These trends are insensitive to changes in parameter values (Supplementary Fig. 1), and while mutations that arise in early generations tend to dominate the mutant population for low PCN, increasing gene dosage contributes to greater genotypic diversity with mutations emerging across multiple generations (Supplementary Fig. 2).

In what follows, we consider the population-level effects of plasmid mutations, emphasizing how circuit-level properties (e.g., PCN, type and location of the mutation, regulatory architecture) modulate the phenotypic mutation rate within a population. While not the focus of our study, it is important to note that the co-occurrence of multiple plasmid variants within a single cell can result in sophisticated dynamics due to the interplay of evolutionary forces acting at two scales: in addition to competition between cells, host fitness is also impacted by replication interference between plasmids harbored within the same host. Studying this latter, intracellular competition, is possible via novel experimental tools such as engineered dimer plasmids50 or via Cas9-based lineage recording51 to quantify the fundamental trade-off between the two levels of selection, but falls outside the scope of our work.

Phenotypic mutation rate depends on gene dosage
To experimentally characterize whether and how gene dosage affects the phenotypic mutation rate, we first considered a simple genetic circuit comprising two plasmids (Fig. 2a). The first one with a p15A origin of replication (ori) contains the CRISPR-guided error-prone DNA polymerase EvolvR46 and sgRNA, both expressed from the arabinose-inducible P_BAD promoter. EvolvR is designed to edit within a window extending to around 60 nucleotides in the 3’ direction from the Cas9 nick site46. The second plasmid is equipped with a pSC101 variant ori, selected from a set with PCN ranging from 3 to 25 per chromosome52, and it harbors constitutively expressed sfGFP with a defect, rendering it non-fluorescent (sfGFP_d gene) by changing residue 88 from methionine to a stop codon. Relying on the expression of araC in the genome (Supplementary Fig. 3), synthesis of EvolvR and the sgRNA is induced via the addition of arabinose to target the stop codon region, thus correcting the defect in the sfGFP_d gene to regain fluorescence.

After co-transformation with both plasmids, cells were induced by arabinose and sfGFP expression was monitored using flow cytometry (Fig. 2b). This revealed that a small fraction of events (less than 1%) showed varying levels of high sfGFP expression, corresponding to bacteria containing mutated plasmids (with sfGFP instead of sfGFP_d). As expected, EvolvR activity is largely confined to the target editing window, correcting the premature stop codon and restoring function resulting in the dominance of a small subset of genotypes (Supplementary Fig. 5b). In general, the phenotypic mutation rate increases with PCN (Fig. 2c) as the potential targets increase, positively contributing to the mutational supply, eventually leading to a greater mutant subpopulation. Unfortunately, PCN also increases the total concentration of the target protein (sfGFP and sfGFP_d combined) for the circuit featured in Fig. 2a. Since we are interested in how the impact of PCN on the mutational supply alone affects the size of the mutant subpopulation, we next integrated a control module to ensure that the expression of the target protein remains unchanged as we explored how the phenotypic mutation rate depends on the gene dosage both in case of gain-of-function and loss-of-function mutations.

Gene regulation impacts the phenotypic mutation rate
To ensure that the expression level of the reporter protein remains unaffected by variations in PCN, we integrated an incoherent feed-forward loop (IFFL) control circuit52 to regulate the synthesis of the target protein (Fig. 3a). Instead of constitutively expressing sfGFP (Fig. 2a), its synthesis is driven by the promoter P_UPsp2, repressed by a transcription activator-like effector (TALEsp2). Changes in PCN affect both TALEsp2 and sfGFP expression, but the impact on the latter can be canceled by fine-tuning the repression of P_UPsp2 via adjusting the TALEsp2 synthesis, yielding constant sfGFP concentration irrespective of the PCN (Supplementary Fig. 9). Integrating the IFFL-based control module therefore decouples the expression of the target gene from PCN variations, allowing the effects of PCN on mutation to be studied independently.

For a given PCN, the concentration of sfGFP within a cell is expected to increase linearly with the number of mutated plasmids (Fig. 3b). Additionally, the strength of each individual P_UPsp2 promoter (i.e., the amount of sfGFP expressed, either defective or functional) is now inversely proportional to the PCN due to repression by TALEsp2, as the total amount is independent of the PCN because of the control circuit. Thus, for a given number of mutated plasmids, the amount of functional sfGFP decreases with PCN (Fig. 3b). As the integration of the control module does not impact the mutation events and the simulation results presented in Fig. 1, we expect that increasing gene dosage results in (i) greater mutant subpopulation (Fig. 1b), but (ii) smaller fraction of mutated plasmids within each mutated cell (Fig. 1d). Therefore, two competing factors emerge that shape the population-level behavior: more mutant bacteria are present as PCN increases, but at the same time, P_UPsp2 repression also becomes more prominent, decreasing sfGFP production in mutated cells (Fig. 3b). Assuming that the detection limit is low and only non-defective sfGFP produces fluorescence, the latter factor is irrelevant, and we would expect the mutant ratio to increase with PCN (gray in Fig. 3c). However, for any non-zero threshold, the initial increase in the mutant ratio may be followed by a subsequent decrease, as typically there are only a handful of mutated plasmids per mutated cell and the amount of sfGFP they each produce decreases with PCN due to the control module. This is precisely what emerges both in stochastic simulations (Fig. 3c) and in experimental data (Fig. 3d).

To further test the model prediction when the detection limit is low, we next replaced the sfGFP_d gene in the circuit outlined in Fig. 3a with a defective ampicillin resistance marker (ampR_d), as even low expression of an antibiotic resistance gene is expected to ensure survival in selective media. The mutated marker was obtained by changing two consecutive residues (104 methionine and 105 serine) with stop codons to minimize possible ribosomal read-through, which could lead to unintended survival in selective media. Expression of this gene remains regulated by the IFFL-based control module, and the sgRNA was redesigned to target the defect-causing stop codons (Fig. 3e). After induction with arabinose to activate EvolvR targeting these stop codons, cultures were plated on carbenicillin selective agar, then survival colonies were counted to estimate the rate of gain-of-function mutations. Fluctuation analysis reveals that the phenotypic mutation rate sharply increases with PCN (Fig. 3d), confirming the model prediction for the low detection threshold case (gray in Fig. 3c).

Phenotypic loss-of-function mutation rate decreases with gene dosage
To reveal how the phenotypic mutation rate changes with gene dosage for loss-of-function mutations, we next included a LacI-repressible sfGFP cassette using the P_Tac-lacO promoter on the plasmid with a p15A ori carrying EvolvR (Fig. 4a). Together with the control module regulating its expression, the lacI gene is harbored on the other plasmid with a pSC101 variant ori. Arabinose induction activates EvolvR, and by targeting the DNA-binding domain of lacI with the sgRNA yields its mutated variant LacI_d. This protein either misfolds or loses its binding affinity to P_Tac-lacO, thus interfering with its inhibitory function on sfGFP expression (Fig. 4b). The total concentration of the regulator (LacI and LacI_d together) modulating the expression of the fluorescent reporter is independent of PCN due to the control module.

For a given PCN, LacI_d is expected to increase linearly with the number of mutated plasmids (as sfGFP did in Fig. 3b), with decreasing sensitivity as the gene dosage increases due to the control module. Since the fraction of mutated plasmids per mutated cell decreases with PCN (Fig. 1d), more LacI will remain functional (Fig. 4c), mitigating the mutation impact and yielding a downward pressure on sfGFP expression (Fig. 4c). Thus, while the number of mutated bacteria increases with PCN, they each express negligible amounts of sfGFP as the non-mutated lacI genes can compensate, hence the population size above the detection threshold is expected to decrease with PCN (Fig. 4c), confirmed experimentally in Fig. 4d.

Mutations in the coding region are masked at the phenotypic level, whereas those in the regulatory region become more prevalent as PCN increases
Autoregulated motifs play a fundamental role in modulating gene expression and shaping evolutionary adaptation10,53–58. Given their wide-spread nature in both natural systems and synthetic biology applications, we next focused on the phenotypic mutation rate in the presence of this central motif. Additionally, since the behavior of a gene circuit may be compromised due to mutations that occur in either the coding sequence of a regulator or its cognate binding site59,60, we designed sgRNAs that target either of these sites.

In particular, we first implemented sfGFP expression from a promoter that is repressed by LacI, which is synthesized constitutively (Fig. 5a). Both of these genes are harbored on the same plasmid equipped with an ori selected from the set of pSC101 variants to explore the impact of changes in PCN. Next, we placed sfGFP and lacI in a polycistronic configuration within the same operon (Fig. 5b) to implement negative feedback. In both of these cases, the EvolvR-associated sgRNA targets either the DNA binding domain of the lacI gene or its cognate promoter.

When LacI is expressed constitutively from the pSC101 variants (Fig. 5a), its concentration increases with PCN, whereas sfGFP synthesis is governed by two opposing forces. While greater PCN results in an upward pressure due to the increased gene dosage of the sfGFP cassette, its activity is also down regulated due to the elevated prevalence of the repressor LacI. While mutating either the lacI gene or its cognate promoter results in elevated sfGFP expression and thus in the emergence of a mutant subpopulation, the way this is achieved is fundamentally different, yielding considerable discrepancies in the phenotypic mutation rate and how it depends on PCN. In particular, mutating the lacI gene results in protein misfolding or impaired binding. However, other plasmids with intact lacI can easily compensate, hence mitigating or altogether eliminating the overall impact of the mutation. Furthermore, this compensation effect is expected to increase with PCN as LacI becomes overabundant, thus saturating the available binding sites and sharply decreasing the average sfGFP per mutant bacteria. This in turn makes it difficult for mutations to manifest at the phenotypic level when the coding region is targeted (Fig. 5a), as the number of wild-type variants that are ready to compensate increases with PCN (Fig. 1). Conversely, mutations in the binding site that prevent its repression via LacI will not be masked by intact copies of the promoter: even a single copy with compromised binding affinity to LacI will express sfGFP independent of the other copies, thus manifesting the mutation at the phenotypic level. Greater PCN increases the mutational supply, which in turn yields increased phenotypic mutation rate when the promoter is targeted (Fig. 5a).

Note that the distribution of mutated cells and mutated plasmids within each cell are independent of the genetic layout, i.e., whether LacI is expressed constitutively (Fig. 5a) or repressing its own synthesis (Fig. 5b). As a result, the observations uncovered in Fig. 1 hold for both circuits in Fig. 5: the number of mutated hosts increases with PCN, but the fraction of mutated plasmids within each cell decreases. Therefore, the phenotypic mutation rate is expected to follow the same trend when lacI is self-repressed (Fig. 5b) as in the constitutive case (Fig. 5a).

To confirm these predictions experimentally, we implemented both circuits outlined in Fig. 5 on various pSC101 variant plasmids (each with a different PCN), and co-transformed them with a p15A plasmid harboring EvolvR alongside with the sgRNA targeting either the lacI coding region or its cognate promoter (Fig. 6a). Additionally, we integrated a TetR-based repression module into the two-plasmid system to repress both LacI and sfGFP expression during EvolvR activation. We included this additional layer of control as Cas9-based systems are known to suffer from off-target effects, especially when binding sites becomes saturated due to the high relative prevalence of Cas9 with respect to its target61–63. This likely occurs when the sgRNA targets the promoter, as LacI binding to lacO blocks EvolvR from accessing the promoter64, diverting EvolvR to unintended sites and resulting in off-target mutations (Supplementary Fig. 6). By modifying the promoters of the circuits in Fig. 5 to include additional repression via TetR, the transcription of sfGFP and lacI are repressed upon arabinose induction both when lacI is expressed constitutively (Fig. 6b) and when it is self-repressed (Fig. 6c). Data in Fig. 6bc confirm the theoretical predictions underpinned by the simulations in Fig. 5: mutations in the coding region are masked at the phenotypic level, whereas those in the regulatory region become more prevalent as PCN increases.

Methods
Parameter values of the modified Wright-Fisher model
The mutation rate μ depends primarily on the concentration of the EvolvR fusion protein and the sgRNA, as well as on the kinetics of DNA binding. Considering prior estimates46, we chose μ = 2.5 ⋅ 10−6 mutations per nucleotide per generation. The constant population size was selected to be N = 105 in all simulations. In our experimental setup, cells are grown at 30°C for 22 h with arabinose induction, which corresponds to approximately 15 generations given that E. coli doubling time is about 90 min at this temperature89. Parameter values of the ordinary differential equations governing the deterministic dynamics of each circuit are provided in Supplementary Sections 4, 5.

Stochastic plasmid replication
At each generation, new plasmids are produced in each cell following a Poisson distribution with the mean equal to the PCN of the corresponding pSC101 variant, as has been previously developed for R1 plasmids90 since both systems use RepA-mediated plasmid replication to control PCN. While the specifics of the replication mechanisms differ in R1 and pSC101 plasmids (co-regulation of RepA synthesis via CopA and CopB in the former, and negative autoregulation of RepA expression in the latter), they are both underpinned by the same feature: replication commences when sufficient RepA initiator protein binds to the iterons within the origin of replication, and the process is regulated via the repression of RepA to stabilize PCN.

Bacterial strains
Clonings were performed using E. coli NEB Stable (C3040H). Experiments in Figs. 2, 3 were carried out using DH5α from NEB (C2978H), while those in Figs. 4 and 6 were performed using DH5αΔlacI::ampR. To generate DH5αΔlacI::ampR, we replaced the genomic lacI gene in DH5α with an ampR cassette using homologous recombination via the λ-Red system91, following the protocol described in ref. 88. The pRED-Cas9 plasmid carrying the λ-Red was a gift from Tao Chen (Addgene plasmid #71541)92. The genome sequence was confirmed via whole genome next generation sequencing (NGS) (Novogene Singapore). Detailed genotypes for the strains used in the study are specified in Supplementary Table 2.

Media and growth conditions
LB media was prepared using premix powder from Sigma Aldrich (L3522). Antibiotics were sourced from Fisher Scientific: kanamycin sulfate (10031553), chloramphenicol (10324980) and carbenicillin (10396833). Final working concentrations are 50 μg/mL for kanamycin, 25 μg/mL for chloramphenicol, and 100 μg/mL for carbenicillin. L-arabinose was obtained from Sigma Aldrich (A3256), and its final working concentration is 2 mM for experiments in Figs. 2–4 and 200 μM for those in Fig. 6. Anhydrotetracycline (aTc) was obtained from IBA Lifesciences (2-0401-002) with a final working concentration of 100 ng/mL.

Cloning and plasmid assembly strategies
The plasmid harboring the EvolvR gene is constructed based on pEvolvR-enCas9-PolI3M-TBD46 (Addgene plasmid #113077). The IFFL control circuit (Addgene plasmid #109254) together with the set of pSC101 variants (Addgene plasmids #109240–#109246) are adopted from ref. 52. All plasmids are built using standard molecular cloning procedures93, including PCR amplification, restriction digestion, gel electrophoresis, ligation, and transformation. Primers and gene blocks were synthesized by Integrated DNA Technologies. All constructed plasmids in the study were confirmed by either Nanopore sequencing (Plasmidsaurus) or NGS (Sangon Biotech). A comprehensive list of all plasmids and DNA parts can be found in Supplementary Tables 3, 4, together with the plasmid maps in Supplementary Fig. 26.

Transformation procedure
Heat shock transformation was carried out using chemically treated NEB Stable competent cells for cloning, and E. coli DH5α or DH5αΔlacI::ampR for experiments. To optimize transformation efficiency, co-transformations were performed in two steps. Following different transformation with the p15A plasmid, cells were grown the next day, then made chemically competent to transform with the pSC101 variant plasmid.

Design of sgRNA sequences
Four sgRNAs were used in the study, targeting (i) sfGFP_d in Fig. 2a and Fig. 3a; (ii) ampR_d in Fig. 3e; (iii) the DNA binding domain of lacI in Fig. 4a and Fig. 6bc; and (iv) the P_tetO-lacO promoter in Fig. 6bc. All sgRNA sequences were examined using ViennaRNA94 for folding. To minimize undesired secondary structures, the sgRNA targeting ampR_d and lacO contain a scaffold with a hairpin lock to stabilize the structure95. The other sgRNAs contain the unmodified Streptococcus pyogenes scaffold. The sgRNA sequences and the corresponding targeting regions are annotated in Supplementary Table 3. All four sgRNA targeting regions have also been checked in CRISPRscan96, showing sufficient efficiencies (Supplementary Table 6).

General protocol for flow cytometry experiments
Co-transformants of EvolvR and pSC101 variant plasmids are grown overnight on LB agar plates at 37°C with appropriate antibiotics and 0.2% glucose (to suppress the P_BAD promoter). Colonies are picked on day 1 for overnight growth at 37°C in LB liquid media supplemented with antibiotics and 0.2% glucose. On day 2, each overnight culture is diluted approximately 10,000-fold into 1 mL LB with antibiotics and 2 mM arabinose, except for experiments in Fig. 6, where 200 μM arabinose is used for induction. To ensure consistent growth, the culture is first diluted approximately 30-fold to reach an OD of 0.15, and then diluted another 333-fold. Next, each 1 mL diluted culture is distributed into 8 wells of a 96-well plate, with 105 μL per well. The plate is inserted into a microplate reader (Biotek Cytation 5M) to incubate with shaking at 30°C for 22 h (induced culture). Meanwhile, each overnight culture from day 1 is diluted 10,000-fold into LB media with antibiotics and 0.2% glucose, and grown at 37°C for 5 h (pre-induced culture). On day 3, induced cultures are diluted 5000-fold into LB media with antibiotics and 0.2% glucose, then grown at 37°C for 6 h or longer to reach the exponential phase. For experiments in Fig. 6, 100 ng/mL aTc is supplemented in the media during the post-induction regrowth stage to relieve the repression from residual TetR, as well as the pre-induction stage for consistent comparison. Both pre-induced and induced cultures were analyzed with a flow cytometer (Thermo Fisher Attune NxT) and the FlowCal Python package97 with at least 100,000 events for each sample. For details on gating and identifying mutant events, see Supplementary Figs. 7, 8.

Fluctuation assay protocol for restoring antibiotic resistance
We use the following protocol adapted from ref. 98 for the experiment in Fig. 3d to restore the function of ampR_d. On day 0, colonies of co-transformed plasmids are picked to grow in LB media with appropriate antibiotics (kanamycin and chloramphenicol) containing 0.2% glucose at 37°C. On day 1, the overnight culture is diluted 10,000-fold into LB media with antibiotics (kanamycin and chloramphenicol) and 2 mM arabinose. The diluted culture is transferred to a 96 well-plate by distributing it over wells, each containing 105 μL volume, and incubated with shaking in a microplate reader (Biotek Cytation 5M) at 30°C for 22 h. On day 2, induced cultures are plated on selective agar (containing carbenicillin) and non-selective agar (containing kanamycin and chloramphenicol). For selective agar plating, 15 μL induced culture from each well is mixed with 150 μL LB containing carbenicillin, then plated on 12-well plates each containing 1.5 mL agar (with carbenicillin). For non-selective agar plating, the induced culture is diluted 1,000,000-fold, of which 500 μL is plated onto an agar Petri dish (15 mL agar). On day 3, colony forming units (CFU) are counted and used as the input to calculate the mutation probability utilizing the R package flan99. The experiment is repeated three times on different days.

Sample preparation for genome and plasmid sequencing
To prepare genome samples, co-transformants of EvolvR and pSC101 variant plasmids are picked to prepare overnight cultures at 37°C, which are subsequently induced at 30°C with 2 mM arabinose for 22 h. Next, induced cultures are either sorted for high sfGFP mutants (for the sgRNA targeting sfGFP_d, lacI, and lacO) using a cell sorter (BD FACSAria III), or grown on carbenicillin containing agar plates to isolate mutated colonies (for the sgRNA targeting ampR_d). Mutated cells are subsequently grown overnight at 37°C and collected for genomic DNA using the QIAamp DNA kit (56304). In particular, the regrown culture is first concentrated to form a pellet and stored at −80°C for 1 h. Then 2 μL protease K (NEB P8107S) is added to the pellet culture. The mixture is incubated at 37°C for 2 h, followed by incubation at 56°C for 1 h, then heat inactivated at 80°C for 5 min. The mixture is cleaned following the steps suggested by the QIAamp DNA kit using the provided buffers and spin columns. The purified genome samples are verified on an agarose gel and shipped on dry ice for NGS genome sequencing (Novogene Singapore). Variant calling is performed on genome sequencing results, using a minimal quality score of 20 and a minimal read depth of 10 for variant detection, as recommended in ref. 100. Results are summarized in Supplementary Table 1.

To prepare plasmid samples for mutation analysis, we induce cultures for 22 h at 30°C with 2 mM arabinose, sort the cultures for high sfGFP expressing mutants using a cell sorter (BD FACSAria III), and collect at least 1000 events for each culture. Sorted cells are regrown on agar plates to isolate individual sfGFP-expressing colonies for plasmid preparation. At least 15 colonies are collected for each sgRNA (targeting sfGFP_d, lacI, and lacO), and whole plasmids are subsequently analyzed using Nanopore sequencing (Plasmidsaurus).

Statistical analysis
Centers in all figures represent mean values, except for the boxplot figures where the centers represent medians. Error bars correspond to the 95% confidence intervals throughout both the manuscript and the Supplementary Information, except for the boxplot figures. The lower and upper bounds of the box represent the first quartile (Q1) and third quartile (Q3), respectively. With IQR = Q3–Q1 denoting the interquartile range, the whiskers extend from the box to the minimum and maximum non-outlier values such that the lower and upper thresholds are defined as Q1−1.5IQR and Q3+1.5IQR, respectively.

Software tools
For data analysis and visualization, Python (version 3.10.14) was used in the Jupyter Lab interface (version 4.2.2). FlowCal python library (version 1.3.0) was used to analyze flow cytometry data. R (version 4.1.2) and the flan package (version 0.9) were used for fluctuation analysis. Adobe Illustrator (version 28.5) was used for creating figures. Microsoft Word for Mac (version 16.62) was used to create the tables presented in the Supplementary Information. Overleaf (version v2) online LaTeX editor was used to prepare and compile the manuscript.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_2": "This paper aims to study the effect of plasmid copy number and various regulatory motifs affect the rate of mutation of a phenotypically silenced fluorescent protein using a two plasmid system with the CRISPR-based EvolvR mutagenesis system. The authors have developed a reasonable strategy for interrogating the phenotypic mutation rate, but there are several questions that arise from their experimental design.\nFirst, it's not clear how the authors are tracking the genotype of individual plasmids within each cell. As the authors point out, the copy number of a plasmid can cause there to be a handful to dozens of the plasmid in a cell. As the plasmid replicates, there is intrinsic mutational variation or heterogeneity in the plasmid sequence from plasmid to plasmid within a cell (see the work of Michael Baym for example). Its not clear how the authors take this heterogeneity into account but it seems intrinsic to the goal. For example, variations in the promoter sequence or the terminator region, or even regions outside of the reporter gene that change the rate of transcription initiation could affect the measured outcomes beyond the action of just the EvolvR system. In Section 7 of the supplement, this is partially addressed by bulk NGS experiments on what seems to be the plasmids (also, what about the genome?... see my next comment below), but this is an bulk, averaged measurement of a culture and does not constitute genotyping of each plasmid. This could be clarified.\nSecond, it's not clear whether EvolvR finalizes its mutation in a single plasmid at the end of a gene of interest or if it continues onwards as a processive polymerase to replicate and mutate other elements of the plasmid. This probably is stated somewhere in the paper, but I could not find any discussion on the point, which seems to be a critical premise in the analysis. How do we conceptually model the EvolvR's mutational action - just confined to the gene of interest or potentially spilling over into the other elements of the plasmid? This seems to be a critical point.\nThird, there will be basal mutation rates in the plasmid and the genome that confound the measured mutation rates from sfGFP expression. These effects will potentially be exacerbated by off-target binding of the EvolvR system, as many Cas-based systems have been shown to have non-specific or off-target binding to other genetic elements on the genome. None of the experiments/data in the figures of the manuscript quantify the genetic drift of the genome or attempt to control for this possibility. To assume it is entirely absent seems a bit too strong of assumption for the goals of the study.\nLastly, it seems like having two plasmids in the same cell with a CRISPR system on one would impose a potential metabolic burden on the cell, especially after activation with a metabolically relevant compound such as arabinose. I did not see strong evidence of the orthogonality of these two-plasmid systems to the fitness and metabolism of the cell. Some basic control experiments measuring fitness or impact would improve the rigor of the paper's argument.\nIn general, it seems like a thoughtful series of arguments in the introduction or perhaps the results section are missing to justify the rigor of the experimental approach presented. The results are interesting and merit publication, but there are some minor yet basic, conceptual revisions here that need to be addressed prior to publication.\nMinor comments:\nIt is not clear if the content in Figures 3 and 4 are wholly based off of stochastic simulation anaylsis or experiment. The delineation between experiment and simulation could be clearer in Figure captions.", "reviewer_3": "This manuscript tackles a novel question in synthetic biology on the impact of plasmid mutations on the behavior of circuits and regulatory motifs, and how this relationship is itself impacted by circuit topology and plasmid copy number. This study has the potential to help significantly improve the robustness of synthetic circuits to mutations, as well as improve performance in directed evolution applications. I find that the authors could better introduce the novelty of their study in the introduction and its potential impact. The authors could also better introduce the motivation behind each of the topologies they test throughout the study. The manuscript is otherwise well written and the results are presented clearly.\nMy main concern is that the authors use a statistics-based method to identify the phenotypic mutation rate in their cultures (Supplementary Fig. 2). This method seems sound and it is at the center of most results in this manuscript, but it needs to be described more thoroughly both in the Main text and in the Supplementary text. I assume they are using the pre-induced cultures that are grown in parallel, mentioned in the Methods section (l. 360-362), as reference for their t-tests but it is not explicitly said anywhere in the text. However the pre-induced cultures surprisingly do not follow the same protocol as the induced cells, as they grow under different temperatures and durations, which might bias the t-tests. The flow cytometry data collected for these pre-induced reference populations should also be shown in Supplementary Fig. 2 and 6-12.\nMinor comments:\nLine 84: \"After their duplication (STEP 2 in Fig. 1a), plasmids from...\" the duplication step appears to be deterministic, where the number of plasmids exactly doubles, and all plasmids are duplicated. This is not consistent with actual plasmid replication mechanisms, did the authors consider making the duplication probabilistic in their model?\nOn Figures 2-6 the authors should describe what the error bars represent and give the number of replicates (8?). For all phenotypic mutation plots, the authors should also consider testing for statistical significance between, at least between PCN=3 and PCN=39/79.\nFigure 3a & 3e: \"p15A harboring EvolvR and sfGFP\" -> sfGFP is not on that plasmid.\nLines 237-244: This part is speculative and the same issue should also happen in Fig. 5b. The authors should either test the binding affinity of their mutated promoters, or rephrase this section.\nIn Supplementary Fig. 5, the authors should specify which plasmid copy number variant they used. This Figure should also be referenced earlier in the Main text to support the validity of their approach."}, "author_response": "Response to Referees\nWe thank the Editor and the Reviewers for their constructive feedback that enabled us to enhance the rigorandclarityofourmanuscript. Weundertooksubstantialexperimental,computational,andeditorial revisionstoaddressthehelpfulcommentsthatwehavereceived(highlightedinblue). Inparticular:\n1. Expandedexperimentalvalidation: Weperformednewwhole-plasmidsequencingofindividual mutant colonies to address plasmid-level heterogeneity, quantified background mutation rates us- ing genome sequencing, and measured metabolic burden via growth rate assays. These additions directly addressed concerns about mutation specificity, off-target effects, and potential confound- ingfactors.\n2. Clarified EvolvR mutational scope: We explicitly defined EvolvR’s editing window, provided experimentalevidenceonmutationconfinement, andimplementedaTetR-basedrepressionmod- ulewithreducedinductionlevelstomitigateoff-targetactivity.\n3. Enhancedcomputationalanalyses: Weintegratedstochasticplasmidreplicationintothemodel, modeledplasmidgenotypedistributionswithinsinglecells,andprovideddetailedjustificationsfor parameterchoices. Thesechangesimprovedbiologicalrealismwithoutalteringkeyconclusions.\n4. Strengthenedmethodologicaltransparency: Weexpandeddescriptionsofthephenotypicmuta- tionrateestimation,includedallpre-inductionflowcytometrycontrols,clarifiedstatisticaltesting procedures, andupdatedfigurecaptionstoexplicitlydistinguishbetweenexperimentalandsimu- lateddatasets.\n5. Improved contextual framing and rationale: The Introduction and Discussion were expanded tobetterarticulatethenovelty,impact,anddesignmotivationforspecificcircuittopologies(inco- herentfeedforwardloopandnegativeautoregulation),inadditiontolinkingourresultstobroader syntheticbiologyanddirectedevolutionprinciples.\n6. Refined data presentation: Figures were reorganized to separate simulations from experiments, replicate numbers and error bar definitions were added, and statistical comparisons across condi- tionswereincluded.\nCollectively,webelievethatthesechangesaddressallReviewerconcerns,reinforcethescientificrigorof ourapproach,andclarifythebroadersignificanceofthefindings. Inadditiontotherevisedmanuscript, we also uploaded a copy indicating all the changes (generated using the latexdiff package: added textisblue,discardedtextisred).\nThispaperaimstostudytheeffectofplasmidcopynumberandvariousregulatorymotifsaffecttherateof mutationofaphenotypicallysilencedfluorescentproteinusingatwoplasmidsystemwiththeCRISPR- based EvolvR mutagenesis system. The authors have developed a reasonable strategy for interrogating thephenotypicmutationrate,butthereareseveralquestionsthatarisefromtheirexperimentaldesign.\nWe thank the Reviewer for recognizing the merit of our strategy and for providing us with constructive and helpful comments regarding our experimental design. In the revised manuscript, we have substan- tially expanded both experimental and computational analyses to address these concerns. Specifically, wenow:\n1. Address plasmid heterogeneity: We performed additional simulations (Supplementary Fig. 2) to model the distribution of genotypes across plasmid copies within single cells, and conducted whole-plasmidsequencingfromisolatedmutantcoloniesreplacingbulkNGSmeasurements(Sup- plementaryFigs.5–6).\n2. Clarify EvolvR’s mutational scope: We explicitly state that EvolvR edits within ∼60 nt of the Cas9 nick site, and present sequencing data showing that mutations are predominantly on-target, withonlyoccasionalspilloverintoadjacentplasmidelements(SupplementaryFigs.5–6).\n3. Quantifybackgroundmutationrates: Usingliterature-derivedratesandgenomesequencingof induced samples, we show that EvolvR activity causes virtually no genomic mutations under our experimentalconditions(SupplementarySection2.3,SupplementaryTable1).\n4. Assess metabolic burden: We now include growth rate measurements (Supplementary Fig. 4) showingthatarabinoseinductionintheabsenceofEvolvRdoesnotaffectgrowth,whereasEvolvR expressiondoesimposeameasurableburden.\nFinally, we reorganized the Results and Discussion sections to improve the logical flow and rigor, and updated figure captions to clearly distinguish between experimental and simulation data. We believe theseadditionsresolvetheconceptualuncertaintiesidentifiedandstrengthenthefoundationofourstudy.\nFirst, it’s not clear how the authors are tracking the genotype of individual plasmids within each cell. As the authors point out, the copy number of a plasmid can cause there to be a handful to dozens of the plasmid in a cell. As the plasmid replicates, there is intrinsic mutational variation or heterogeneity in the plasmid sequence from plasmid to plasmid within a cell (see the work of Michael Baym for example). Its not clear how the authors take this heterogeneity into account but it seems intrinsic to the goal. Forexample,variationsinthepromotersequenceortheterminatorregion,orevenregionsoutside of the reporter gene that change the rate of transcription initiation could affect the measured outcomes beyondthe actionofjust theEvolvR system. InSection7 ofthesupplement, thisis partiallyaddressed by bulk NGS experiments on what seems to be the plasmids (also, what about the genome?...see my next comment below), but this is an bulk, averaged measurement of a culture and does not constitute genotypingofeachplasmid. Thiscouldbeclarified.\nFrom a computational perspective (Supplementary Section 1), to estimate the distribution of genetic variations at the plasmid level, we first note that EvolvR mutates the target at a rate of 2.5 · 10−6 per nucleotidepergeneration,suggestingthatmutationsaresparseacrossgenerations. Thisalignswellwith the scenario underpinning fluctuation analysis, in which a few early-generation mutations expand into descendants. Consequently, many mutant plasmids are expected to share similar genotypes from a few mutatedancestorsthatariseearlyon. Toinvestigatetheissueofintrinsicplasmidheterogeneity,wecate- gorizedmutatedplasmidsbasedonthegenerationinwhicheachmutationfirstoccurred. Weassumethat eachplasmidmutatesatmostonce,astheEvolvR-associatedsgRNAlosesbindingaffinitytoamutated target. Thus,eachmutatedplasmidisassignedtoauniquegeneration. Thedistributionofthesegenerations(thateventuallygivesrisetothedistinctmutationsweobserveinthefinalsimulatedgeneration)is presentedinSupplementaryFig.2. Aswenowstateinthemanuscript(lines101–103),thesedatareveal that “while mutations that arise in early generations tend to dominate the mutant population for low PCN,increasinggenedosagecontributestogreatergenotypicdiversitywithmutationsemergingacross multiplegenerations(SupplementaryFig.2).”\nFromanexperimentalperspective(SupplementarySection2.2),sincenext-generationsequencing(NGS) measurements characterize mutation genotypes that emerge only in bulk at the population level, we isolated individual colonies exhibiting mutant phenotypes and prepared the plasmids for sequencing. In particular, cells were grown after arabinose induction, and following regrowth, the ones showing GFP expression were sorted and regrown on agar plates (see Methods). Thus, plasmids in the resulting individual colonies are expected to be genetically homogeneous, hence sequencing them provides fine- grained genetic information that is not available via bulk NGS. As it is now included in the manuscript (lines109–110): “EvolvRisdesignedtoeditwithinawindowextendingtoaround60nucleotidesinthe 3’directionfromtheCas9nicksite(47),”togetherwithadditionaldetailsinSupplementarySection2.2 to acknowledge that mutations can potentially occur beyond this range with decreasing efficiency and thatthesizeofthiswindowdependsontheprocessivityoftheerror-proneDNApolymeraseinEvolvR (PolI3M), and can be further modulated via the insertion of a thioredoxin-binding domain (TBD) into thepolymerase(PolI3M-TBD).Mutationsthatoccurwithintheintendedregionareconsideredon-target, otherwise they are off-target. On-target ratios are calculated by dividing the on-target mutation counts with respect to the total number of mutations. As we now state in the manuscript (lines 120–122): “EvolvRactivityislargelyconfinedtothetargeteditingwindow,correctingtheprematurestopcodonand restoringfunctionresultinginthedominanceofasmallsubsetofgenotypes(SupplementaryFig.5b).”\nFinally,wesurveyedtherelevantliteratureassuggestedbytheReviewerandacknowledgeintheDiscus- sion(lines299–302)that“theroleandimpactofintrinsicplasmidheterogeneitycanbefurtherstudied using advanced experimental techniques such as induced plasmid monomerization (90) or via Cas9- based lineage recording (91) to genetically encode plasmid replication events.” This line of inquiry, however,liesoutsideofthescopeofourmanuscript.\nToclarifytheexpectedbehaviorofthemutagenesisplatform, wenowexplicitlystateinthemanuscript (lines109–110)that“EvolvRisdesignedtoeditwithinawindowextendingtoaround60nucleotidesin the3’directionfromtheCas9nicksite(47),”alongsidewithfurtherdetailsinSupplementarySection2.2. Additionally, to examine whether mutations are confined to this window, we sequenced whole plasmid DNAsfromindividualcoloniesthatshowmutantphenotypes(asdetailedinthepreviouscomment). For thesgRNAtargetingGFP (Fig.2aandFig.3a)wenowincludethat(lines120–122)“EvolvRactivityis d\nlargelyconfinedtothetargeteditingwindow,correctingtheprematurestopcodonandrestoringfunction resultinginthedominanceofasmallsubsetofgenotypes(SupplementaryFig.5b).” Complementingthis, wealsoacknowledge(Fig.4caption)that“Whileoff-targetmutationsbecomemoreprevalentincaseof LacIthanwhentargetingGFP ,themajorityoftheseunintendedmutationsarestillrestrictedtothelacI d\ngene,givingrisetothesamemutationphenotype(SupplementaryFig.5).” Finally,weobservedecreased EvolvR specificity for the sgRNA recognizing the lacO binding sites (Supplementary Fig. 6), which is notsurprising,andisnowstatedinthemanuscript(lines233–237): “asCas9-basedsystemsareknown to suffer from off-target effects, especially when binding sites become saturated due to the high relative prevalence of Cas9 with respect to its target (65–67). This likely occurs when the sgRNA targets the promoter, as LacI binding to lacO blocks EvolvR from accessing the promoter (68), diverting EvolvR to unintended sites and resulting in off-target mutations (Supplementary Fig. 6).” To address this issue (lines231–232),“weintegratedaTetR-basedrepressionmoduleintothetwo-plasmidsystemtorepress both LacI and GFP expression during EvolvR activation” for the constructs in Fig. 6. Additionally, as we now include in the caption of Fig. 6: “Arabinose induction is reduced from 2 mM to 200 µM to decreasetheactivityofEvolvRandthechanceofitsbindingtooff-targetsites. Subsequenttoarabinose induction, anhydrotetracycline (aTc) at a final concentration of 100 ng/mL is supplemented to relieve the effects of residual TetR.” As noted in the caption of Fig. 6c, “On-target mutations are found at a few specific positions that are functionally important, disrupting the lacO binding sites, resulting in the dominanceofasmallsubsetofgenotypes(SupplementaryFig.6a),”alongsidewiththeacknowledgment that“evenwiththeinclusionoftheTetRrepressionmoduleandthereducedEvolvRexpression,off-target mutationspersistforpSC101variantswithPCNbelow9.” Therefore,thecorrespondingdataareomitted inFig.6c,andinsteadincludedinSupplementaryFigs.17–18.\nThe data in Supplementary Figs. 5–6 alongside with the extensive characterization during the original developmentofEvolvRillustrate(aswenowstateintheDiscussion,lines286–292)that“EvolvRactivity ismostlyconfinedtothegeneofinterest(thoughnotnecessarilytotheeditingwindow),causingvirtually nogenomicmutations,onlyoccasionallyspillingoverintootherelementsoftheplasmid. Ourdataalso highlightthatbesidesintroducinganadditionalcontrollayer(suchasourTetR-basedrepressionmodule) orreducingEvolvRexpression,increasingtheabundanceofthetargetsiteviaPCNcontrolcanofferan effectivealternativestrategytocombatthisundesirablecharacteristicofCas9-basedtools.”\nComplementingtheworkoutlinedintheprevioustwocommentsaddressingplasmidheterogeneityand EvolvR’s targeting accuracy around the editing window, we included additional computational analysis and experimental verification to estimate the basal mutation rates in the plasmid and the genome that couldconfoundthemeasuredmutationrates.\nFromacomputationalperspective,consideringpreviouslyreportedmutationrates,weestimatelessthan a single nucleotide mutation in the genome throughout our experiments. In particular, as we detail in Supplementary Section 2.3, E. coli has a basal spontaneous mutation rate of approximately 10−10 mutations per nucleotide per generation. The Cas9 variant in EvolvR (enCas9-PolI3M-TBD) is engi- neeredsothatithasreducednon-specificDNAbinding,resultinginon-targetandglobalmutationrates of approximately 10−5 and 10−8 mutations per nucleotide per generation, respectively. We utilized a P -induciblesystemto tightlyregulateEvolvRandsgRNA expression, hencetheapproximateback- BAD\nground mutation rate can be estimated as follows. Prior to induction, cultures are grown for around 24 hours (i.e., 32 generations at 37◦C) on agar (12 hours) and liquid media (12 hours) with 0.2% glucose (Fig. 2b) to suppress basal P expression through catabolite repression, ensuring a low spontaneous BAD\nmutation rate during this stage. Subsequently, we induce EvolvR expression using arabinose for 22 hours at 30◦C, i.e., for approximately 15 generations. Considering the typical E. coli genome size, this is expected to generate approximately 5·106 ·(10−10 ·32+10−8 ·15) = 0.766 mutations, i.e., less than a single mutation throughout the experiment (and even fewer for the plasmids considering their maximal size and PCN: 13·12·103 ·(10−10 ·32+10−8 ·15) = 0.024 for p15A with 13 copies and 100·8·103·(10−10·32+10−8·15) = 0.123forpSC101withnomorethan100copies).\nTovalidatethisestimate,wecollectedmutatedsamplesafterarabinoseinductionandperformedgenome sequencing to identify potential mutations (targeting GFP , ampR , lacI, and lacO). As the sgRNA tar- d d\ngetingthelacIgeneharboredontheplasmidalsorecognizestheoneintheDH5αgenome,wegenerated theDH5α∆lacI::ampRstrainbyknockingoutthegenomiclacIgene. Followingthis,weperformedvari- antcallingonthegenomesequencesobtainedusingnextgenerationsequencing(seeMethods). Among all tested samples, only a single mutation was detected in the genome (Supplementary Table 1), when thesgRNAtargetsGFP ,yieldingasinglenucleotidepolymorphism(AtoT)inthePqqLgeneencoding d\nfor a zinc protease (turning the 122nd amino acid from Val to Glu; no nearby sequence resembles the sgRNA target, but several NGG motifs are present). These data confirm the low background mutation levelsuggestedbyourcomputationalestimate.\nThisisnowexplicitlystatedintheDiscussion(lines286–288): “EvolvRactivityismostlyconfinedtothe geneofinterest(thoughnotnecessarilytotheeditingwindow),causingvirtuallynogenomicmutations, onlyoccasionallyspillingoverintootherelementsoftheplasmid.”\nTo address this shortcoming, we performed additional microplate reader experiments to quantify the impact of arabinose induction, as well as the metabolic burden that each of the plasmids impose on the host. GrowthratedatainSupplementaryFig.4 highlightthat(i)inductionviaarabinoseintheabsence ofEvolvRonthep15AplasmidorharboringthetargetplasmidexpressingGFP fromapSC101plasmid d\ndonotimpactcellgrowth;whereas(ii)theexpressionofEvolvRfromap15Aplasmiddoes. Thisisnow alsoacknowledgedinthecaptionofFig.2b.\nIngeneral,itseemslikeathoughtfulseriesofargumentsintheintroductionorperhapstheresultssection are missing to justify the rigor of the experimental approach presented. The results are interesting and meritpublication,buttherearesomeminoryetbasic,conceptualrevisionsherethatneedtobeaddressed priortopublication.\nWe appreciate the Reviewer’s constructive feedback and helpful recommendations. As outlined above, we have taken several steps to justify the rigor of our experimental approach: additional data collected during the revision process have been integrated into the manuscript and Supplementary Information to reinforce key claims, which range from growth rate measurements to sequencing analysis. These additions,complementingthechangesthatwemadetothemanuscriptitselftoimprovethelogicalpro- gressionfromhypothesistoexperimentaldesign,andfromdatatoconclusions,renderourfindingsmore transparent to the reader. We believe that together these revisions address the Reviewer’s concern and strengthentheoverallrigorofourstudy.\nMinor comments: It is not clear if the content in Figures 3 and 4 are wholly based off of stochastic simulationanaylsisorexperiment. Thedelineationbetweenexperimentandsimulationcouldbeclearer inFigurecaptions.\nToimproveclarityandthedelineationbetweenexperimentsandsimulations,eachfigurecaptionnowex- plicitlystateswhethertheresultsareobtainedinsimulationsorexperiments,inadditiontothecomplete reorganizationofFigs.5–6toseparatethetwotypesofdata.\nThismanuscripttacklesanovelquestioninsyntheticbiologyontheimpactofplasmidmutationsonthe behaviorofcircuitsandregulatorymotifs,andhowthisrelationshipisitselfimpactedbycircuittopology and plasmid copy number. This study has the potential to help significantly improve the robustness of syntheticcircuitstomutations,aswellasimproveperformanceindirectedevolutionapplications.\nWe thank the Reviewer for their constructive and encouraging feedback, which recognized the novelty and potential impact of our study on understanding how plasmid copy number and regulatory architec- ture influence mutation dynamics in synthetic circuits. These insightful comments have significantly strengthenedthemanuscriptand helpeduspresentourfindings with greaterclarityandrigor. Inpartic- ular, as we detail in our point-by-point response, we made the following changes to address the helpful pointsraisedbytheReviewer:\n1. Expanded Introduction and Discussion: We provided further context to articulate the novelty andpotentialimpactofthestudy,emphasizinghowourresultscaninformsyntheticbiologydesign anddirectedevolutionstrategies, aswellastobetterconnecttheexperimentalandcomputational resultstobroaderdesignprinciplesforbuildingmutation-resilientcircuitsandtuningevolvability.\n2. Motivationforcircuittopologies: Weincludedtherationaleforintegratingtheincoherentfeed- forwardloopandnegativeautoregulationmotifsintoourconstructs,explainingtheirrelevancefor modulatinggeneexpression,evolutionaryadaptation,andmitigatingmetabolicburden.\n3. Clarifiedphenotypicmutationratemethodology: Weincludeddetailedexplanationsdescribing thestatisticalframework,controls,anddatasourcesforfluctuationanalysis. Wealsoincorporated pre-inductionflowcytometrydataforcompletenessandtransparency.\n4. Additional experimental and computational validation: We conducted further sequencing ex- perimentstocharacterizetheoff-targetimpactofEvolvR,integratedanadditionalcontrolmodule tomitigatetheseeffects,andincorporatedstochasticplasmidreplicationintothesimulations.\n5. Improved data presentation: We updated figure captions to clearly distinguish simulation from experimental datasets, added the number of replicates and error bar definitions, performed pair- wise statistical tests between PCN conditions, and explicitly stated the PCN throughout both the manuscriptandtheSupplementaryInformation.\nWe sincerely thank the Reviewer for their thoughtful and detailed feedback, which has directly con- tributedtoimprovingboththescientificdepthandclarityofourmanuscript. Theirsuggestionsprompted ustoexpandourcontextualframing,refineourexperimentalandcomputationalmethodologies,andim- provethetransparencyofourdatapresentation. Furthermore,theencouragementtoprovidedeepermo- tivationforourdesignchoicesandtovalidateourfindingsthroughadditionalexperimentshasstrength- ened the robustness of our conclusions. The revisions inspired by these comments not only address the specific points raised but also broaden the scope and potential impact of our study within the fields of syntheticbiologyanddirectedevolution.\nToaddressthisshortcoming,werevisedboththeIntroductionandtheDiscussiontobetterintroducethe noveltyofourstudyanditspotentialimpact.\nWe first incorporated the following in the Introduction (lines 53–56): “Harnessing EvolvR, we imple- mentedacollectionofengineeredcircuitstostudytheemergenceofplasmid-bornemutations. Weshow thattheinterplaybetweengenedosageandregulatorymotifsshapesbothgenotypicandphenotypichet- erogeneity,highlightingthevalueofintegratingsyntheticgenecircuitsandinvivomutagenesistostudy evolutionary dynamics.” Additionally, we also included in the Discussion (lines 251–256) that “it is imperative to understand how circuit-level features link genetic mutations to the spread of mutants in the presence of inevitable evolutionary forces to uncover the organizing principles of living cells, and to devise design guidelines that ensure the predictable behavior of engineered biocircuits. This is es- peciallycrucialwhenevolvinggeneticcircuitsatanintermediatescalebetweenfocusingonoptimizing singlegeneticcomponentsandfacilitatinggenome-wideadaptationsinseriallypropagatedcellpopula- tions(76).”\nTheseadditionsclarifythenoveltyofourworkandbettercontextualizeandhighlightthepotentialimpact of our results, complementing the closing paragraph of the Discussion (lines 305–312): “Our work provides explicit design guidelines by revealing how regulatory architecture and gene dosage via PCN togethershapetheemergenceofmutantswithinapopulation,despitethesimplicityofourcomputational model (e.g., mutations can be neutral or even increase protein stability and enhance binding (92, 93)) and the imperfect nature of the mutation platform (e.g., off-target effects of EvolvR (47)). Therefore, the results presented here shed light on evolutionary driving forces, and can be leveraged to guide the designofevolutionarilystablesyntheticgenecircuitsortodynamicallyprogramgeneticdiversification byleveragingincreasedevolvability(88,89),forinstance,viaflexiblePCNcontrol(81,94).”\nFinally, we provide further motivation and detail how the the regulatory motifs used in our study (i.e., IFFLandnegativeautoregulation)impactcircuitperformanceandevolvabilitythroughoutthemanuscript:\n• lines191–192: “Autoregulatedmotifsplayafundamentalroleinmodulatinggeneexpressionand shapingevolutionaryadaptation(10,54–59).”\n• lines 263–267: “Integrating an additional control layer (like the IFFL-based regulatory module or negative autoregulation) can relax such bioenergetic constraints (57), allowing the circuits to persistonhighcopyplasmids. Thisinturncandirectlyimpacttheirevolvabilityamongotherben- efits(58,59),likelycontributingtotheubiquitousnatureofthesemotifsintranscriptionnetworks (82,83).”\n• lines 276–280: “Additionally, inclusion of a control layer can modulate the mapping from geno- typicmutationstophenotypicmutationratedependingonthedetectionthreshold,whichcouldbe crucialwhenestablishinggeneticdominance(85). Thatis,byregulatinggeneexpression,wecan adjustthevisibilityandsignificanceofinevitablemutations,effectivelyalteringthedominanceof genetictraits(85).”\nWethankthereviewerforpointingoutthisshortcoming. Toclarifythesesourcesofconfusion,wehave madethefollowingthreechangesthroughoutthemanuscriptandtheSupplementaryInformation.\nFirst, we provide all pre-induction flow cytometer data alongside post-induction measurements in Sup- plementaryFigs.19–25. Growthofbothculturesinvolvedilutionfollowedbyregrowthat37◦Cpriorto FACSanalysis. Themajordifferencebetweenpre-andpost-inductionworkflowsisthatthelatterrequires longer regrowth, caused mainly because EvolvR proteins accumulate during induction and continue to exist in cells, and the corresponding metabolic burden negatively impacts cell growth (Supplementary Fig.4). Tomitigatethisstress,weadopted30◦Cduringtheinductionstageinsteadof37◦C(asindicated in Fig. 2b and explicitly stated in the figure caption). Additionally, the regrowth time may also vary depending on several other factors, including final cell density after induction, manual dilution, batch media quality and the uneven temperature distribution in shaker incubators. We therefore allow a flexi- bleregrowthwindowthatspansapproximately6hours,rangingfrom5to7hours. Thisapproachensures that culturesreach comparablegrowth states inexponential phasebefore subsequent FACSanalysis for mutant detection (which is also included in the caption of Fig. 2). The protocol is included in Fig. 2b for easy comparison of the two phases, and it is explicitly stated in the figure caption that the mutant subpopulationisidentifiedusingflowcytometerdatacomparingthepre-andpost-inducedsamples.\nSecond, SupplementaryFig.8hasalsobeenrefinedforenhancedclarity, includingthefollowinginthe caption: “Cultures were screened prior to arabinose induction as a control to exclude spontaneously mutatedindividualsandtoprovideareferenceformutantdetection.”\nThird, we revised the description of how mutants are detected in Supplementary Section 3. Key modi- ficationsaretheinclusionofthedetailsthatunderpinthestatistics-basedmethodforidentifyingmutant cells,aswellasfurtherstepsonthecalculationofthephenotypicmutationrate. Inparticular,welever- aged the flow cytometer data to quantify mutant cells, which were then used to calculate phenotypic mutation rate using fluctuation analysis. The mutation rate represents the probability of mutation per generation,estimatedfromthedistributionoftheobservedmutantsatthefinaltimepoint.\nMinorcomments:\nLine 84: “After their duplication (STEP 2 in Fig. 1a), plasmids from...” the duplication step appears to be deterministic, where the number of plasmids exactly doubles, and all plasmids are duplicated. This is not consistent with actual plasmid replication mechanisms, did the authors consider making the duplicationprobabilisticintheirmodel?\nin a given cell, respectively, where P is the PCN. Similarly, let N and N = N − N denote the m wt m\nnumberofmutant(withatleastonemutatedplasmid)andwild-typebacteria(thosewithoutanymutated variants), respectively. We generatethe numberof newmutations∆xwithineach cellas follows: each wild-typeplasmidbecomesmutatedwithprobabilityµ,modeledasaPoissonprocess. Next,thenumber of wild-type and mutated plasmids for each cell are updated (STEP 1 in Fig. 1a) as x′ = x −∆x wt wt\nand x′ = x +∆x, respectively. Subsequently, the number of new plasmids P′ after replication is m m\ngenerated according to a second Poisson distribution, with its mean equal to the average PCN (see Methods),andthispoolissplitaccordingtoabinomialdistributionwithparametersP′ andx′ /x′ , m wt\nthat is, newly replicated plasmids are randomly assigned to be either wild-type or mutated according to the respective proportions within the cell (STEP 2 in Fig. 1a). Following this, both daughter cells receive either (P +P′)/2 plasmids if P +P′ is even, otherwise one of them ends up with one more thantheothertoensurethatthereisnosegregationalloss, whereahypergeometricdistribution(50) determines how mutated and wild-type plasmids are distributed (STEP 3 in Fig. 1a). Finally, half of thepopulationisrandomlyselectedtoinitializethenextgeneration(STEP4inFig.1a),whiletherestis discardedtokeepthepopulationsizeconstantacrossgenerations.”\nAdditionally,tojustifythechoiceofthePoissondistributiongoverningthestochasticnatureofplasmid replication, we have included the following in the Methods (under “Stochastic plasmid replication”): “At each generation, new plasmids are produced in each cell following a Poisson distribution with the mean equal to the PCN of the corresponding pSC101 variant, as has been previously developed for R1 plasmids (96) since both systems use RepA-mediated plasmid replication to control PCN. While the specifics of the replication mechanisms differ in R1 and pSC101 plasmids (co-regulation of RepA synthesis via CopA and CopB in the former, and negative autoregulation of RepA expression in the latter), they are both underpinned by the same feature: replication commences when sufficient RepA initiatorproteinbindstotheiteronswithintheoriginofreplication,andtheprocessisregulatedviathe repressionofRepAtostabilizePCN.”\nFinally, we updated Fig. 1a to reflect the above changes in the mathematical model, as well as the simulationdatathroughoutboththemanuscriptandtheSupplementaryInformation. Importantly,while incorporating stochastic plasmid replication into our computational framework captures an important aspectoftheunderlyingcellularprocesses,itdoesnotimpacttheoveralltrendsandourconclusions.\nOn Figures 2–6 the authors should describe what the error bars represent and give the number of repli- cates (8?). For all phenotypic mutation plots, the authors should also consider testing for statistical significancebetween,atleastbetweenPCN=3andPCN=39/79.\nTo address these shortcomings, we made the following changes. First, to clarify how the phenotypic mutation rate is calculated and the meaning of the error bars, we updated Fig. 2c. In particular, we in- 10\nclude data on the detected mutant ratio for each PCN, and state in the figure caption that “The mutant ratioisobtainedinflowcytometerexperiments(SupplementaryFig.8)using8replicates,whichisthen leveragedusingfluctuationanalysis(52–53)toderivethephenotypicmutationrate(errorbarsrepresent thestandarddeviationoftheestimatedmutationprobability).” Additionally,inSupplementarySection3 we include a detailed explanation of applying fluctuation analysis to estimate the phenotypic mutation rate from mutant detection experiments using flow cytometer data. Finally, we assessed the statistical significanceoftheobservedvariationinthephenotypicmutationrateacrossPCNconditionsusingpair- wise t-tests. The resulting p-values are presented as heatmaps in Supplementary Fig. 11–18, which are referencedinthecorrespondingfigurecaptionsofthemanuscript.\nFigure3a&3e: “p15AharboringEvolvRandsfGFP”→sfGFPisnotonthatplasmid.\nWethankthereviewerforpointingoutthiserror,thefigurehasbeenupdatedaccordingly.\nLines 237–244: This part is speculative and the same issue should also happen in Fig. 5b. The authors shouldeithertestthebindingaffinityoftheirmutatedpromoters,orrephrasethissection.\nMotivated by this helpful comment, we performed additional experiments. We found that while muta- tions are mostly confined to the editing window when targeting GFP or lacI (Supplementary Fig. 5), d\nthe sgRNA targeting the lacO promoter tends to cause off-target mutations including in the lacI gene, resultinginphenotypessimilartolacOmutants(SupplementaryFig.6).\nAs we detail in Supplementary Section 2.2: “Cas9-based systems are known to suffer from off-target effects,especiallywhenbindingsitesbecomesaturatedduetothehighrelativeprevalenceofCas9with respect to its target (4). Therefore, it is not surprising that we observe decreased EvolvR specificity for thesgRNArecognizingthelacObindingsites,especiallywhenthePCNharboringthetargetsequenceis low(SupplementaryFig.6a). Tounderstandwhythishappens,unlikewhentargetingGFP orlacI,note d\nthatLacIbindingtolacOblocksEvolvRfromaccessingthepromoter(5),divertingittounintendedsites and resulting in off-target mutations, for instance, in the lacI gene (Supplementary Fig. 6bc). Inducers suchasIPTGcanrelieveLacIfrombindingtolacO,however,thisleadstoconstitutiveGFPexpression, which would accumulate over generations and interfere with downstream mutant detection using flow cytometry. Thus,weinsteadintegratedaTetRrepressionmoduleintothetwo-plasmidsystemtorepress both LacI and GFP expression during EvolvR activation (Fig. 6a). By blocking LacI production this way, the lacO binding sites become accessible for EvolvR to initiate on-target mutations. Additionally, we lowered the expression level of EvolvR by using 200 µM arabinose instead of 2 mM to reduce the chance of EvolvR binding to off-target sites. Just as before, on-target mutations are found at a few specificpositionsthatarefunctionallyimportant,inthiscase,todisruptthelacObindingsites,resulting inthedominanceofasmallsubsetofgenotypes(SupplementaryFig.6a).”\nThisisnowalsobrieflyincludedinthemanuscript(lines231–240): “Additionally,weintegratedaTetR- based repression module into the two-plasmid system to repress both LacI and GFP expression during EvolvR activation. We included this additional layer of control as Cas9-based systems are known to suffer from off-target effects, especially when binding sites become saturated due to the high relative prevalence of Cas9 with respect to its target (65–67). This likely occurs when the sgRNA targets the promoter,asLacIbindingtolacOblocksEvolvRfromaccessingthepromoter(68),divertingEvolvRto 11\nunintendedsitesandresultinginoff-targetmutations(SupplementaryFig.6). Bymodifyingthepromoters ofthecircuitsinFig.5toincludeadditionalrepressionviaTetR,thetranscriptionofGFPandLacIare repressed upon arabinose induction both when LacI is expressed constitutively (Fig. 6a) and when it is self-repressed (Fig. 6b).” Additionally, in the caption of Fig. 6 we also state that “Arabinose induction is reduced from 2 mM to 200 µM to decrease the activity of EvolvR and the chance of its binding to off-targetsites.”\nAsexpectedbasedonthesimulationdatapresentedinFig.5,thehighlighteddiscrepancybetweencon- stitutive and self-repressed LacI expression (previously Figs. 5b and 6b, now Figs. 6b and 6c) is now resolved when relying on the redesigned genetic circuits in Fig. 6a. In particular, as we state in the manuscript (lines 240–243): “Data in Fig. 6bc confirm the theoretical predictions underpinned by the simulationsinFig.5: mutationsinthecodingregionaremaskedatthephenotypiclevel,whereasthose intheregulatoryregionbecomemoreprevalentasPCNincreases.”\nInSupplementaryFig. 5,theauthorsshouldspecifywhichplasmidcopynumbervarianttheyused. This FigureshouldalsobereferencedearlierintheMaintexttosupportthevalidityoftheirapproach.\nIn the revised manuscript, we now state the copy number of each plasmid variant when presenting the sequencingdatathatweleveragedtocharacterizeplasmidheterogeneityandtoquantifythebackground mutationrate,followingthehelpfulrecommendationsofReviewer#2.\nFirst, we isolated individual colonies exhibiting mutant phenotypes and prepared the plasmids for se- quencing due to the limited resolution of NGS measurements characterizing mutation genotypes that emergeonlyinbulkatthepopulationlevel. Accordingly,SupplementaryFig.5isnowreplacedbySup- plementary Figs. 5–6. To obtain these data, we sorted the arabinose induced cultures, regrew them on agarplates,andisolatedthecoloniesshowingmutantphenotypes. Plasmidswerethencollectedfromin- dividualcoloniesforwholeplasmidsequencing(seeMethods). ThecaptionsofSupplementaryFigs.5–6 explicitlystatethePCNofthevariantsweused.\nSecond, to validate our computational estimate of the background mutation rate (Supplementary Sec- tion 2.3),wecollectedmutatedsamplesafterarabinoseinductionforgenomesequencingtoidentifypo- tentialmutations(targetingGFP ,ampR ,lacI,andlacO).Followingthis,weperformedvariantcalling d d\nonthegenomesequencesobtainedusingNGS(seeMethods). Theresultsoftheanalysisaresummarized inSupplementaryTable1,alongsidewiththePCNofthevariantsweused.\nTakentogether,thesequencingdatasupportthevalidityofourapproach,whichisnowhighlightedearlier whenintroducingtheinvivomutagenesisplatform(lines120–122): “EvolvRactivityislargelyconfined to the target editing window, correcting the premature stop codon and restoring function resulting in thedominanceofasmallsubsetofgenotypes(SupplementaryFig.5b),”aswellasemphasizedagainin the Discussion (lines 286–288): “EvolvR activity is mostly confined to the gene of interest (though not necessarily to the editing window), causing virtually no genomic mutations, only occasionally spilling overintootherelementsoftheplasmid.”\n12"}, {"round": 2, "reviewer_comments": {"reviewer_2": "I appreciate the authors thoughtful responses, their additional experiments, and their dutiful follow-up on the questions I raised regarding intracellular plasmid variability. The authors have addressed the majority of my concerns. My primary concern is still regarding the framing of this paper. Figure 1 introduces the theme of this paper as intracellular plasmid variability and presents beautiful simulations to quantify the statistics of accumulation of mutation. As I see it (and I may have missed something), no other figure in the paper continues on this vein, but rather speaks towards general phenotypic outcomes of mutation using the EVOLVR system cleverly developed by the authors.\nThe supplementary figure 5 shows the outcomes of whole plasmid sequencing experiments on individual colonies, which are presumed to be isogenic variants in the traditional sense. However, the issue of intracellular plasmid variability as introduced in Figure 1 is still unaddressed, in this regard. A single cell can have multiple mutants of a plasmid and whole plasmid sequencing, NGS methods do not quantify this intrinsic genetic variation but rather present a measurement of a bulk sample of hundreds of thousands of cells.\nI agree with the authors that development of methodologies to measure such intrinsic variability would be outside the scope of this work, but I think Figure 1 is quite confusing, in that it sets the expectation that the focus of this work will be on that theme. The subsequent figures are fine and tell a strong story, but the framing of the paper just needs some adjustment. No additional experiments are required, as I can tell.\nOne minor comment: In Supplementary Figure 5 (and elsewhere throughout the paper), standard error is used to plot the error bar. This is misleading, as the standard error represents the standard deviation normalized by an additional sampling factor n — (sqrt(n) to be precise). I recommend the authors use 95% confidence intervals instead, as standard deviation is often too generous of a measure while standard error is too conservative. Using a confidence interval is generally statistically more rigorous, as well.", "reviewer_3": "The authors have thoroughly and convincingly addressed my comments and concerns. I recommend this study for publication in Nature Communications."}, "author_response": "Response to Referees\nWe thank the Editor and the Reviewers for their constructive feedback that enabled us (i) to enhance the clarity of our manuscript by refining its framing, and (ii) to improve its rigor via updated statistical analysis.\nReviewer #2\nIappreciatetheauthorsthoughtfulresponses,theiradditionalexperiments,andtheirdutifulfollow-upon thequestionsIraisedregardingintracellularplasmidvariability. Theauthorshaveaddressedthemajority ofmyconcerns. Myprimaryconcernisstillregardingtheframingofthispaper. Figure1introducesthe themeofthispaperasintracellularplasmidvariabilityandpresentsbeautifulsimulationstoquantifythe statistics of accumulation of mutation. As I see it (and I may have missed something), no other figure inthe papercontinues onthis vein, but ratherspeaks towardsgeneral phenotypicoutcomesof mutation usingtheEVOLVRsystemcleverlydevelopedbytheauthors.\nThesupplementaryfigure5showstheoutcomesofwholeplasmidsequencingexperimentsonindividual colonies, which are presumed to be isogenic variants in the traditional sense. However, the issue of intracellularplasmidvariabilityasintroducedinFigure1isstillunaddressed,inthisregard. Asinglecell canhavemultiplemutantsofaplasmidandwholeplasmidsequencing,NGSmethodsdonotquantifythis intrinsic genetic variation but rather present a measurement of a bulk sample of hundreds of thousands ofcells.\nI agree with the authors that development of methodologies to measure such intrinsic variability would be outside the scope of this work, but I think Figure 1 is quite confusing, in that it sets the expectation thatthefocusofthisworkwillbeonthattheme. Thesubsequentfiguresarefineandtellastrongstory, but the framing of the paper just needs some adjustment. No additional experiments are required, as I cantell.\nToclarifytheframingofthepaper,weincludedthefollowingrightafterpresentingthesimulationdatain Fig.1: “Inwhatfollows,weconsiderthepopulation-leveleffectsofplasmidmutations,emphasizinghow circuit-levelproperties(e.g., PCN,typeandlocationofthemutation, regulatoryarchitecture)modulate the phenotypic mutation rate within a population. While not the focus of our study, it is important to note that the co-occurrence of multiple plasmid variants within a single cell can result in sophisticated 1\ndynamics due to the interplay of evolutionary forces acting at two scales: in addition to competition betweencells,hostfitnessisalsoimpactedbyreplicationinterferencebetweenplasmidsharboredwithin the same host. Studying this latter, intracellular competition, is possible via novel experimental tools suchasengineereddimerplasmids(51)orviaCas9-basedlineagerecording(52)toquantifythefunda- mental trade-off between the two levels of selection, but falls outside the scope of our work.” This way we (i) emphasize that the present work focuses on population-level phenotypic effects of plasmid mu- tations, which are first revealed computationally and then subsequently confirmed experimentally; and (ii) acknowledge that high PCN can introduce competition both within and across cells, but addressing this lies beyond the scope of our work. We believe that this revision eliminates the source of confusion identifiedbytheReviewerandprovidesacleararticulationofthemanuscript’scentralfocus.\nOne minor comment: In Supplementary Figure 5 (and elsewhere throughout the paper), standard error is used to plot the error bar. This is misleading, as the standard error represents the standard deviation normalized by an additional sampling factor n — (sqrt(n) to be precise). I recommend the authors use 95%confidenceintervalsinstead,asstandarddeviationisoftentoogenerousofameasurewhilestandard erroristooconservative. Usingaconfidenceintervalisgenerallystatisticallymorerigorous,aswell.\nFollowingthishelpfulrecommendation,weupdatedthestatisticalanalysissothaterrorbarscorrespond to the 95% confidence intervals throughout both the manuscript and the Supplementary Information, except for the boxplot figures where the lower and upper bounds of the box represent the first quartile (Q1) and third quartile (Q3), respectively. With IQR = Q3-Q1 denoting the interquartile range, the whiskers extend from the box to the minimum and maximum non-outlier values such that the lower and upper thresholds are defined as Q1-1.5IQR and Q3+1.5IQR, respectively. This is now explicitly mentionedintheMethodsaswellasinthefigurecaptions.\nReviewer #3\nTheauthorshavethoroughlyandconvincinglyaddressedmycommentsandconcerns. Irecommendthis studyforpublicationinNatureCommunications.\nWethanktheReviewerfortheirconstructivefeedbackandsuggestionsthathelpedusnotonlyimprove theresultspresentedintheoriginalsubmission, butalsosparkednovellinesofresearchinquirysignifi- cantlyexpandingthescopeofourworkthroughoutthereviewprocess.\n2"}]

Decision letter