Evolution of a minimal cell – Nature

Strains and growth conditions

We maintained synthetic M. mycoides JCVI-syn1.0 and synthetic M. mycoides JCVI-syn3B in SP4 medium with KnockOut Serum Replacement (Gibco) substituted for fetal bovine serum (Supplementary Table 3). Cultures of these non-motile bacteria were grown in a dark, static growth chamber at 37 °C. The non-minimal JCVI-syn1.0 strain has been described in detail previously⁴⁵. The minimal JVCI-syn3B is identical to the strain synthesized in previous studies³ with the following exceptions: JVCI-syn3B possesses a second rRNA operon copy, lacks a gene (MMSYN1_0531) encoding an efflux protein, and has 19 genes that were added back into the minimal genome to render the cell easier to use^4,33 (Supplementary Table 4). The strain also contains a landing pad system (cre recombinase and loxP) facilitating genetic manipulation. For competition experiments used to quantify relative fitness, we used a JCVI-syn1.0 strain that expresses mCherry, which enabled us to distinguish it in mixed culture from other strains using flow cytometry and also factor out any costs associated with production of the fluorescent protein (see below).

Mutation accumulation experiment

Overview

Mutation accumulation (MA) experiments are designed to reduce the influence of natural selection through repeated bottlenecks of evolving populations¹⁹. When used with microbial populations, this is typically achieved by transferring single colonies, which have undergone single-cell bottlenecks. Before initiating MA experiments, we acclimatized JCVI-syn1.0 and JCVI-syn3B to laboratory conditions by maintaining populations in SP4 liquid medium. We took a clone of each acclimated strain to begin the MA experiment and propagated replicate lineages (n = 87 and n = 57 for JCVI-syn1.0 and JCVI-syn3B, respectively) for 20 to 36 weekly transfers.

Number of generations

To compare rates of mutation across replicates, we normalized all rates as per-generation values. To calculate the number of generations per transfer in the MA, we grew cells on SP4 agar for 1 week and diluted a sample of seventh day colonies into 1 ml of phosphate-buffered saline (pH 7.4). Cells were fixed with 20 μl of 25% glutaraldehyde and stained with 2× SYBR Green, and then counted with a NovoCyte flow cytometer (ACEA Biosciences). We used the dilutions to calculate the number of cells in the original colony, from which we inferred the number of generations (log₂[N], where N is the number of cells in the undiluted colony) that must have occurred to reach a colony of that size⁴⁶, assuming each colony is formed by a single progenitor cell. As the growth rate and other fitness components can decrease during an MA experiment⁴⁷, we also measured the number of cells per colony during and at the end of the MA, averaging across timepoints to estimate the total number of generations. We then used the number of generations per transfer to estimate the effective population size (N_e) using the harmonic mean method⁴⁷. Specifically, N_e was approximated as the harmonic mean of the series (2⁰, 2¹, 2², …, 2^f), where f is equal to the number of generations per transfer inferred from the previous step.

Whole-genome sequencing and sequence analysis

We performed DNA extractions from evolved MA cell lines using the DNeasy UltraClean Microbial Kit (Qiagen) according to the manufacturer’s instructions, with the additional step of adding 50 μl of 50 mg ml⁻¹ lysozyme to improve cell lysis. Genomic DNA was sequenced using Illumina MiSeq sequencing to a depth of at least 35× coverage. Library preparation and DNA sequencing were conducted by the Indiana University Bloomington Center for Genomics and Bioinformatics. Whole-genome sequencing reads were quality controlled using cutadapt⁴⁸ to trim low-quality base pairs and remove residual adapter sequences. We used breseq with the default parameters^49,50 to call mutations using the trimmed reads. We only considered fixed mutations for the MA cell lines. We checked for mutations that had arisen in experimental ancestor strains before evolution. Ancestral mutations were removed from the analysis of all evolved MA lines derived from that strain using gdtools^49,50. We used the sequencing data to check for contamination or cross-contamination in the evolved cell lines.

Statistical analyses

To compare the mutation rate and spectrum between strains, we used two-sample t-tests for numerical response variables and two-sample χ² tests with continuity correction for comparing proportions. For comparing proportions to theoretical expectations within a strain, we used one-sample χ² tests with continuity correction.

Adaptive evolution

Overview

In contrast to the mutation accumulation experiments, we conducted experiments that allowed bacteria to achieve large population sizes to increase the efficacy of natural selection. This involved serial passaging of cells in liquid cultures with limited bottlenecking at each transfer. For example, in our experiment, the minimum population size was 2× 10⁷–4 × 10⁷ for both JCVI-syn1.0 and JCVI-syn3B. We passaged replicate 3 ml liquid cultures of each strain (n = 4 per strain) in 13 mm glass test tubes by 1% (v/v) serial transfer each day for 300 days in a dark, static incubator held at 37 °C. We calculated the number of generations per day as the log₂ of the dilution factor, that is, log₂[101], the number of binary fissions needed to regenerate the original population size after the 1% (v/v) transfer⁵¹. Thus, we estimate that the M. mycoides strains were maintained for 1,997 generations, which, based on other experiments, is long enough for the majority of adaptation to occur^51,52.

Measurements of fitness

First, we measured fitness as µ_max by conducting growth curves on cells that were isolated at different timepoints during the adaptive evolution experiment (Supplementary Fig. 5). Cryopreserved cells were thawed on ice before preculturing at 37 °C for 24–72 h in 3 ml of SP4 medium in a 13 mm test tube. Before initiating the experiment, we adjusted the start times of precultures to help ensure that cultures from different evolution timepoints were at the same stage of growth. Approximately 6 × 10⁵ cells from turbid precultures were then inoculated into replicate wells of a 96-well plate containing 200 µl of SP4 medium. Separately, each population was incubated in a 96-well plate for 24 h in a BioTek Synergy H1 microplate reader that recorded the absorbance every 15 min at 415 nm. This wavelength is close to a spectral peak for phenol red, a pH indicator that is a component of SP4 medium (Supplementary Table 3). Previous studies have shown that phenol red can be used as proxy for metabolism and growth⁵³ because bacteria like M. mycoides produce organic acids as a byproduct of carbohydrate metabolism⁴ (Supplementary Fig. 5). With the resulting data, we used maximum likelihood to estimate growth-curve parameters using a modified Gompertz equation⁵⁴:

where L is the lag time (h), A is the carrying capacity or yield (optical density at 415 nm), µ_max is the maximum growth rate (day⁻¹) and b₀ is the intercept (Supplementary Fig. 6 and Supplementary Table 5).

Second, we measured relative fitness by competing ancestral and evolved strains against a M. mycoides JCVI-syn1.0 reference strain labelled with mCherry (syn1.0::mCh)²⁶. Cryopreserved cells were used to make precultures in a similar manner to those in the growth curve experiment. Each strain was grown in liquid medium to log phase, and then the labelled and unlabelled strains were simultaneously diluted into a mixed culture in fresh medium. We immediately sampled the axenic cultures or the mixed culture (t₀), fixed the cells with 20 μl of cold 25% glutaraldehyde, incubated them at 4 °C for 20 min and then stained the samples with 2× SYBR Green. After 24 h of growth (t_f), the mixed culture was sampled and processed again in an identical manner. For samples in the adaptive evolution experiment, we quantified the abundance of each strain using a an LSR II flow cytometer (BD Biosciences) at Indiana University’s Flow Cytometry Core Facility. For measuring the relative fitness of engineered ftsZ mutants, we used the NovoCyte flow cytometer (ACEA Biosciences). While measurements were being made, we vortexed the samples every minute to prevent multiple cells from clumping together and being scored as single events. The purity was assessed during every run using negative controls and axenic controls. We detected 1,800–2,700 events per second and abundances on the order of 1 × 10⁸ cells per ml. With the resulting data, we differentiated cells on the basis of the expression of mCherry. Using NovoExpress, FACSDiva and FCS Express software, we established gates on pure cultures of the non-mCherry-expressing experimental strains and the syn1.0::mCh reference strain (Supplementary Figs. 7 and 8). For the experimental strains, boundaries were established by gating axenic mCherry-negative cells that were positive for only SYBR Green fluorescence. For the reference strain, boundaries were established by gating axenic syn1.0::mCh cells that were positive for SYBR Green and mCherry (Supplementary Fig. 9). In the competition assays used to quantify relative fitness, we applied the axenically established gates to samples that contained a mixture of the reference strain and experimental strain. We obtained the proportion of false-negative mCherry cells by applying the mCherry-negative gate to axenic mCherry-expressing cells; this proportion was then used as a correction factor in mixed populations. Last, we calculated relative fitness as the change in the relative abundance of the strain of interest during the 24 h period of competitive growth versus syn1.0::mCh. Specifically, the relative fitness versus the mCherry reference strain W_C is

where N₀ represents the initial abundance of the experimental strain, N_f the abundance of the experimental strain after 24 h, and N_C0 and N_Cf are initial and final abundances of the reference strain (syn1.0::mCh), respectively²⁶. We normalized fitness values to be relative to the original M. mycoides JCVI-syn1.0 ancestor strain. In other words, we represent the fitness (W) as \(\frac{{W}_{C}}{{W}_{{\rm{J}}{\rm{C}}{\rm{V}}{\rm{I}}-{\rm{s}}{\rm{y}}{\rm{n}}1.0}}\), where W_{JCVI -syn1.0} is the value of W_C for M. mycoides JCVI-syn1.0.

Whole-genome sequencing and sequence analysis

DNA extraction, sequencing and bioinformatics were performed according to the same methods as for the mutation accumulation experiment with a few exceptions. Specifically, each replicate population was sequenced to a depth of at least 100× coverage, and polymorphic mutations were included in our analyses. As an indicator of selective pressure, we used the Jukes–Cantor method⁵⁵ to compute the per-site d_N/d_S value on the basis of the number of nonsynonymous and synonymous SNMs within each of the evolved replicate populations normalized by the total nonsynonymous and synonymous target sizes. We counted the number of synonymous and nonsynonymous AT to CG, AT to GC, AT to TA, CG to GC, CG to TA and CG to AT sites using the gdtools module of breseq, which is a computational pipeline that identifies mutations from short-read DNA resequencing studies⁵⁰. We next combined that information with the empirical mutation spectrum from the MA experiment to account for the differing probabilities of each of the six SNM types, and thereby calculate the total expected number of SNMs at nonsynonymous and synonymous sites⁵⁶. The observed numbers of synonymous and nonsynonymous substitutions were obtained directly from breseq outputs. Synonymous and nonsynonymous polymorphisms were included in the observed count with probability equal to their allele frequency in mapped reads. We added a pseudocount of 1 synonymous substitution for all calculations⁵⁷ because two of the populations had 0 synonymous substitutions.

To identify mutations possibly contributing to adaptation, we looked for genes that had mutations across two or more replicate populations for each genotype. Mutations in the same gene, arising and increasing in frequency in independent lineages, suggests that that mutation’s rise could be driven by positive selection⁵⁸. To test this hypothesis, we statistically assessed whether multiply-mutated genes (that is, genes mutated in >1 replicate evolved population) had acquired more mutations than would be expected by chance under the assumption that the mutations were neutral⁵⁸. To do this, we recorded all of the polymorphic and fixed mutations that were called within genes. Synonymous mutations were excluded. We then used Python⁵⁹ to simulate the placement of these mutations at random across all genes. The probability of any given gene receiving any given mutation was relativized to the gene’s length and GC content using the known mutation rates of G:C nucleotides and A:T nucleotides from the mutation-accumulation experiment. We repeated this random placement of mutations 100,000 times. In each simulation, we counted the number of mutations received by each gene, with each fixed mutation increasing the count by 1 and each polymorphism increasing the count by an amount equal to its allele frequency. For each multiply-mutated gene from the real adaptation experiment, we calculated the proportion of the 100,000 simulations in which the gene received at least as many mutations as were truly observed and called this proportion the P value. We then used the Benjamini–Hochberg method^60,61 to generate corrected P values (P_adj) to account for multiple tests with the false-discovery rate set to be α = 0.05 (Extended Data Table 2). As a negative control, we repeated the simulations using only synonymous mutations. This process returned two false-positive significant genes, which was small compared with the 52 significant signatures detected among nonsynonymous mutations, although we also acknowledge that synonymous gene analysis had less power due to the smaller number of synonymous mutations.

Generation of ftsZ E315* mutant cells

This process required mutating the bacterial genomes while they were yeast centromeric plasmids (YCPs) followed by genome transplantation of the mutated genomes. The YCPs were mutated using rounds of CRISPR–Cas9 and yeast homologous recombination that is a modification of a method used previously to mutate M. mycoides strains⁶².

In the first CRISPR–Cas9 step, the molecule to be mutated was cleaved and the donor DNA comprising sequences from the two flanking genes was recombined with the cut JCVI-syn1.0 or JCVI-syn3B YCP, removing parts of genes of the flanking genes and all of the target gene. The donor DNA had 40 bp overlaps to both genes flanking the target gene and had a 22 bp Mycoplasma gallisepticum 161 CRISPR–Cas9 target sequence with a protospacer adjacent motif (PAM) (5′-GTATAAATACATCCAGGAGTGG-3′) that had no homology elsewhere in JCVI-syn1.0 or JCVI-syn3B. The M. gallisepticum sequence put a new PAM in the genome that was used in the second round of CRISPR–Cas9.

The second round of CRISPR–CAS9 cut the JCVI-syn1.0 or JCVI-syn3B YCP at the new M. gallisepticum PAM. The cut YCP was then recircularized using a donor DNA containing the desired point mutation. The mutagenized regions of the YCPs were PCR amplified and the mutation was confirmed by Sanger sequencing. Correctly mutagenized JCVI-syn1.0 or JCVI-syn3B YCPs were then transplanted into Mycoplasma capricolum recipient cells as reported previously^{3,59,60,63,64}. The mutagenized regions of the transplants were PCR-amplified and sequenced to confirm the presence of the desired mutations.

Microscopy and image analysis

Scanning electron microscopy (SEM) was used to compare changes in the cell size of evolved populations. All of the populations were grown in the same batch of medium and under identical conditions in a single incubator. The start times of cultures were adjusted so that they reached stationary phase at the same time. We centrifuged stationary-phase cultures and resuspended the pellet in 1 ml of phosphate-buffered saline (pH 7.4). The resuspended cells were fixed by adding 20 μl of cold 25% glutaraldehyde and incubating at 4 °C for 20 min. For microscopy observation, fixed cells were concentrated 4× by centrifugation and resuspension. The centrifugation steps were performed at 25 °C for 4 min at 2,000g. SEM was performed at the Indiana University Bloomington Electron Microscopy Center. Fixed cells in PBS were pelleted and resuspended in 100 mM sodium cacodylate buffer (pH 7.2) with 2 mM calcium chloride and 2% sucrose. We coated 12-mm-diameter glass coverslips with 0.1% poly-l-lysine for 5 min, after which coverslips were washed with a few drops of double distilled water. Resuspended cells were added to the coverslip surface and allowed to adhere. After 5 min, the coverslips were washed twice with 100 mM sodium cacodylate buffer (pH 7.2) with 2  mM calcium chloride and 2% sucrose. Next, 300 µl of 2% osmium tetroxide in 100 mM sodium cacodylate buffer (pH 7.2) with 2% sucrose was added to the surface of the coverslips while on ice. After 30 min, the coverslips were washed with double-distilled water. The coverslips were placed into a CPD coverslip holder (Electron Microscopy Sciences, 70193-01). The samples were dehydrated in a graded ethanol series (30%, 50%, 70%, 90%, 95%) while on ice. At room temperature, the coverslips were rinsed three times with 100% ethanol. Each dehydration step lasted for 2 min. Critical-point drying was performed using the Tousimis Samdri 790 critical-point dryer. The dried coverslips were placed on aluminium SEM stubs and sputter-coated using the Safematic CCU-010 with SP-010 Sputter Head with 45 nm of gold/palladium (80%/20%), which is accurate in the Angstrom range. All of the samples were coated simultaneously to minimize variance among samples. We viewed the samples using the FEI Teneo scanning electron microscope at 2.0 kV, 25 pA probe current and 3.0 mm working distance. The T2 detector was used. We calibrated the measurements using line grating replicas (2,160 lines per mm) with 0.261 μm latex spheres (Electron Microscopy Sciences). We analysed the SEM image data using ImageJ⁶⁵. We used the straight and measure features combined with image scale metadata to measure the vertical diameters of imaged cells that met the following criteria: cells must be round; cells must not have apparent holes or punctures; cells must be completely within the field of view; cells must have an unambiguous perimeter; there must be no suggestion that a cell is currently or has recently undergone binary fission; cells must be ≥0.1 μm across. Each image was processed counterclockwise starting from east. The samples were processed in a randomized order.

Statistical analyses

For the growth-curve experiments, we used a generalized linear mixed model to test for the fixed effects of time (generation) and cell type (minimal versus non-minimal) on growth curve parameters (µ_max, lag time, yield) while fitting random intercepts for the replicate evolved populations (Supplementary Table 3). We used variance partition coefficients to estimate the contribution of the replicate populations (random effect) to the total variation explained in the models (Extended Data Fig. 1, Extended Data Table 1, Supplementary Figs. 1 and 2 and Supplementary Table 5). For the adaptative evolution experiment (Figs. 2 and 4), we tested hypotheses using a general linear model (GLM) after subtracting observations of each replicate-evolved population (generation 2,000) from its corresponding ancestor (generation 0). With the intercept term excluded, the GLM tests whether the evolutionary trajectory for each group is different from zero. With the intercept term included, the GLM tests whether the evolutionary trajectories are different among groups. We also used two-way ANOVA with Tukey’s honest significant difference test to test hypotheses about the effects of cell type (minimal versus non-minimal) and ftsZ E315* (wild type versus mutant) on relative fitness and cell size. When necessary, data were log₁₀-transformed to meet statistical assumptions.

We compared the composition of genes acquiring mutations among the evolved replicate populations by first constructing a gene-by-population matrix. Here, each row represented an evolved population and each column represented a gene that had acquired at least one mutation among all of the populations. Each cell of the matrix was filled with the sum value of mutations occurring in that gene in that population, where fixed mutations were valued at 1 and polymorphisms were valued equal to the allele frequency. Only essential genes, shared between JCVI-syn1.0 and JCVI-syn3B, were considered. We used PERMANOVA on the Bray–Curtis distances generated from the gene-by-population matrix to test for the significance of cell type (minimal versus non-minimal) on the composition of mutations using the adonis function in the R package vegan⁶⁶. For visualization, the Bray–Curtis distances were decomposed into two dimensions using principal coordinate analysis using the cmdscale function.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.