Lineage Tracking for Probing Heritable Phenotypes at Single-Cell Resolution

Determining the phenotype and genotype of single cells is central to understand microbial evolution. DNA sequencing technologies allow the detection of mutants at high resolution, but similar approaches for phenotypic analyses are still lacking. We show that a drop-based millifluidic system enables the detection of heritable phenotypic changes in evolving bacterial populations. At time intervals, cells were sampled and individually compartmentalized in 100 nL drops. Growth through 15 generations was monitored using a fluorescent protein reporter. Amplification of heritable changes–via growth–over multiple generations yields phenotypically distinct clusters reflecting variation relevant for evolution. To demonstrate the utility of this approach, we follow the evolution of Escherichia coli populations during 30 days of starvation. Phenotypic diversity was observed to rapidly increase upon starvation with the emergence of heritable phenotypes. Mutations corresponding to each phenotypic class were identified by DNA sequencing. This scalable lineage-tracking technology opens the door to large-scale phenotyping methods with special utility for microbiology and microbial population biology.


Introduction
Populations of microbes are increasingly used as experimental tools to study evolutionary processes [1,2]. To monitor evolutionary change competition experiments between neutrallymarked evolved strains and ancestral strains are often used [2,3]. These assays provide a measure of the fitness improvement of derived clones or populations, but do not reveal phenotypic or genetic changes. In recent years, technological advances in DNA sequencing have made it possible to track evolution at the level of individual mutations [4][5][6]. However, connecting mutations to their phenotypic effects remains a significant challenge, often requiring years of painstaking genetic and physiological analyses [7,8]. Desirable are methodologies that allow detection of subtle phenotypic differences at the same level of resolution as provided by deep sequencing. To date, a handful of studies have exploited fortuitous differences in colony morphology among variant types and adaptive phenotypes after plating samples on agar media [9][10][11]. Detection of mutant types based on differences in colony morphology requires that adaptive mutants affect morphology and reach sufficient frequency to be detectable by plate culture. If it were possible to track heritable phenotypes in the descendants of many individual cells, then subtle heritable phenotypic differences among cells could be discerned within an evolving population. Ideally, these different cells would be non-destructively sorted and coupled to DNA sequencing to detect and to understand genetic and phenotypic evolution.
Here, we use a drop-based millifluidic system (Fig 1) that allows the compartmentalization of individual cells in~100 nL drops, which function as independent microreactors [12][13][14]. While the growth of cells in emulsion droplets is not new [15,16] the present system allows the descendants of hundreds of independently isolated cells to be tracked for~15 generations. This allows heritable phenotypic changes in individual cells that are undetectable at the single cell level to be exponentially amplified and quantified. To demonstrate the utility of the technology, we monitored the evolution of a single E. coli clone carrying a yellow fluorescent protein (YFP) reporter [17,18] during 30 days of starvation. During prolonged stationary phase, populations are expected to evolve and diversify. Indeed, previous analysis of such populations was instrumental in the discovery of Growth Advantage in Stationary Phase (GASP) mutants [19][20][21][22]. Moreover, in the absence of shaking, evolving populations may become heterogeneous due to alterations in the selective environment wrought by metabolic changes in the founding The resulting populations were maintained in starvation for 0, 1, 3, 6, 15 or 30 days. Populations were then analysed by growing colonies from single cells in drops of 100 nL fresh LB medium. Individual drops were collected in microplate wells for further analysis with complementary methods. (B) Millifluidic droplet analyser principles. Droplets are prepared in a cross junction (see S1 Movie) and stored in a PTFE tube (500μm inner diameter). By controlling the pressure at each side the sequence of drops can be moved back and forth in front of the measurement point (see S2 Movie).
population and spatial heterogeneity [9,23]. In general, depletion of the nutrients by one dominant clone may open opportunities for new mutants to arise that take advantage of changes in ecological opportunity [22,24].

Results
We demonstrate the utility of the millifluidic colony analyser to monitor heritable phenotypic changes at high resolution by studying changes in populations of E. coli propagated under conditions of starvation over the course of 30 days. The experiment involved growth of six replicate populations under static culture conditions, and three replicate populations under shaken conditions, for a total of 30 days. At set intervals (days 0, 1, 3, 6, 15 and 30) independent microcosms were destructively harvested and cells were analysed in order to detect heritable phenotypic diversity (Fig 1A).
The millifluidic colony analyser [12,13] is based on droplet-based microfluidics [14,25,26]. A detailed description of the technology is available in previous publications [12,13]. Two distinct phases are embedded in a continuous flow of fluorinated oil in millimetre-scale tubing (inner diameter of 500 μm). Each nutrient drop reservoir (~100 nL) is separated within the "train" from its neighbour by a hydrocarbon drop of the same size. The train of drops flows continuously ( Fig 1B, S1 and S2 Movies) ensuring that the contents of each drop are continually mixed. This also ensures that a thin layer of fluorinated oil always lubricates the interface between the drops and the surface of the tube. Before injection into the device, bacteria are diluted with fresh medium so that, assuming a Poisson distribution of cells in drops [27,28], the average number of cells per drop, λ, is~0.5 and most droplets contain either no cells, or just a single cell. For a given sample, the founding inoculum size is readily deduced from the fraction of droplets showing growth (see Materials and Methods). About 1,000 droplets were maintained in a single tube, 10 meters long, and approximately one droplet in two displayed growth (~500 inoculated droplets per experiment). Lineage tracking is based on fluorescence measurement. The E. coli clone used to initiate the evolution experiment carries a chromosomally-integrated yellow fluorescent protein (YFP) [17,18]. YFP fluorescence is measured in each drop every 10 minutes (S2 Movie). For each single cell the fluorescence signal is measured during the course of growth. Expression of YFP is driven from a constitutive lac promoter (see Materials and Methods) [17,18], but is nonetheless responsive to changes in the metabolic status of the cell, for example, those wrought by mutation. This provides a sensitive read-out of phenotype [29]. By tracking the descendants of individual cells, minute heritable changes in growth, yield, and YFP expression stand to be amplified, thus allowing discrimination between phenotypes that are otherwise indistinguishable at the single cell level.
At successive time intervals, a sample of cells from static and shaken microcosms was collected and analysed with the millifluidic machine (Fig 2; S1 and S2 Figs). Analysis of~500 occupied drops inoculated with the ancestral E. coli strain (on entry to stationary phase) showed that the increase in fluorescence acquired over 15 generations of growth (from 1 to 10 5 bacteria), is identical for all drops (Fig 2). Indeed, the coefficient of variation (CV) in the millifuidic system (5%) is greatly reduced compared to the CV using flow cytometry (34%) (S3 Fig).
In the first selection experiment based on populations maintained in static broth microcosms (Fig 2), after 1 day, a single phenotype was observed. Fluorescence profiles were homogeneous. At day 3, two distinct profiles were detected: one with the same final fluorescence intensity as the ancestor (blue lines), and a new type with a lower final fluorescence intensity (red lines). At day 6, three profiles were detected: two with the same fluorescence intensity as observed at day 3 (blue and red lines), with one displaying an intermediate intensity (green lines). At day 15, two main classes remained evident (blue and red lines) with the green type no longer detected, however the derived population (red lines) represents 98% of the population. At day 30, two new profiles were detected (purple and yellow lines) with a single droplet harbouring populations manifesting phenotypes of the two previously dominant classes from day 15 (blue and red lines).
Next we asked whether the different fluorescence profiles have a mutational cause. The millifluidic device allows the recovery of individual droplets associated to different fluorescence profiles (S4 Fig). Individual drops corresponding to four different classes (numbered 1 to 4 in Fig 2) from day 6 to 30, were plated on LB-agar. They produced distinct colonies with varied morphology (S5 Fig) that were distinctly different to that of the ancestral strain. Types isolated at day 6 and day 15 of static starvation (phenotypes 1 and 2) had a mucoid colony morphology [30], whereas colonies isolated at day 30, corresponding to phenotypes 3 and 4, showed ancestral-like morphology, but the diameter of each colony was respectively smaller, and larger, than that of the ancestor. To test whether these changes were heritable, cells from droplets expressing altered phenotypes were propagated in LB broth and the phenotype rechecked both in drops and via agar-plate culture. In every case the phenotypes mirrored those of the parental types.
In order to determine the genetic basis of the different profile types, we performed whole genome sequencing on four clones from each of the four phenotypically distinct derived populations (Table 1). Mutations were found within clones representative of each profile. The four clones of phenotype 1 with a mucoid colony had identical genotypes, containing just a single point mutation in yrfF (G99V). Little is known about yrfF but deletion is known to induce mucoidy [31]. All four clones of phenotype 2, which also produce a mucoid colony on agar plates, also contain a mutation in yrfF (Y98D), however the position of this mutation is different from that in phenotype 1. The difference reflects the fact that each test tube is an independent biological replicate, but the gene-level parallelism points to evidence of selection and thus an adaptive effect of these mutations [8,32]. In addition, all replicates of phenotype 2 carry a mutation in wzc (D642G), a further known contributor to the mucoid phenotype [33]. For phenotype 3, all four genotypes contain an insertion of IS5 between csrA (regulator of glycogen synthesis and catabolism) and alaS (alanine synthesis). IS5 movement is known from prolonged stationary phase experiments [34]. In contrast to the first three phenotypes, where all the clones shared mutations, for phenotype 4, one genotype had a unique mutation in barA (D472N) and the other three clones each had a different mutation in mglC. No mutations were observed in lrp or rpoS, which are characteristic of Growth Advantage in Stationary Phase (GASP) mutants [35]. Nevertheless, the mutations in barA and csrA are related to global regulation and are reminiscent of rpoS expression modifications. barA has been reported to regulate rpoS expression [36], while crsA play a role in the regulation of barA [37,38] and represses genes expressed in stationary phase [39].
The individual drops corresponding to phenotypes 3 and 4 from the static microcosms (detected after 30 days of starvation), were also analysed by flow cytometry (Fig 3A) after recovery into wells of a microtitre plate. The difference in mean fluorescence between phenotypes 3 and 4 was almost identical when measured in the millifluidic system or using flow cytometry (Fig 3B), with phenotype 3 being 1.4-fold more fluorescent than phenotype 4. However, the CV of the fluorescent signal obtained using cytometry is large (Fig 3A), reflecting the impact of noise at the single cell level [17,40,41]. The high CV using flow cytometry makes discrimination of different phenotypic variants in the static microcosm (with mixed phenotypes) impossible, whereas they are clearly distinguishable using the millifluidic system.   In Fig 4, the final yield fluorescence intensity levels measured in drops for a total of 30 independent starvation experiments are plotted as a function of time spent under static or shaken conditions. This provides a view of the diversification dynamics in the microcosms. For each sample, one to four discrete phenotypic classes are always present (Fig 2, S2 Fig). Moreover, the normalized fluorescence intensity and the frequency of the different phenotypic classes strongly fluctuate from one sample to another, even within a single time point (S2 Fig). Diversity is particularly large after 30 days, where fluorescence intensities vary more than 10 fold (Fig 4A). As a general feature, the number of phenotypic profiles (classes), their intensities, and the number of members in each class are found to strongly fluctuate from one experiment to another, confirming the existence of a range of possible evolutionary solutions and trajectories [2,3,21,22]. Comparison with the results obtained for unstructured (shaken) microcosms ( Fig  4B) confirms the importance of spatial structure for the process of diversification [9,42]. When using the Shannon index [43] (see Methods) to quantify diversity between treatments (Fig 4C), early-stage diversification is higher under static relative to shaken conditions, but the final diversity is similar after 30 days in both conditions (Fig 4, S1 Fig). Thus, diversity under starvation builds up in about six days and then is dynamically maintained, as already suggested [21].

Discussion
In contrast to approaches based on flow-cytometry [40,41], time-lapse microscopy [17], or plating [44], our millifluidic system enables detection, quantitative analysis and isolation of clones. By tracking a signal reflecting the physiological state of the descendants of individual cells over about 15 generations, i.e. from a single cell to about 10 5 cells, our methodology identifies heritable changes. Amplification across 15 generations overcomes gene expression noise occurring at each generation. According to the law of large numbers the contribution of this noise to the variance of our measurement scales as σ 2 /N, thus signal becomes more apparent with increasing generation number. In essence, the millifluidic system operates by a principle similar to used in agar plate culture as pioneered by Robert Koch in 1881 [45]. However, cells in colonies on plates experience different conditions depending on the proximity of colonies to one another, and depending upon the location of cells within colonies [44,46,47]. In contrast, in the millifluidic system, all individuals express their phenotype in strictly equivalent wellmixed environments. The phenotypic resolving-power of our method relies on the equivalence in growth and measurement conditions and capacity to amplify via growth, phenotypic differences that are indistinguishable at the single cell level.
Application of our methodology relied on the fluorescence signal derived from expression of yfp [17]. Expression of this gene is essentially constitutive except for the regulation by cAMP Receptor Protein (CRP). cAMP is a second messenger whose concentration depends on the growth medium, being in particular, sensitive to glucose levels (cAMP concentration is low when glucose is the carbon source). cAMP is also involved in regulation of numerous metabolic pathways [48] and therefore regulation of yfp by CRP stands as a reporter of the physiological state of the cell, and hence sensitive to many phenotypic changes. Thus our phenotypic trait of interest is the expression of yfp integrated over 15 generations. The detection of discrete classes in the starved static microcosm confirms the usefulness of the yfp reporter (Fig 2). The classes identified are strongly reminiscent of the consequences of evolution underpinned by genetic (mutational) change. Indeed, the emergence of clones with a fitness advantages associated with a heritable change would account for the occurrence of such classes. Although YFP fluorescence is not a direct measurement of fitness, phenotypic changes are nevertheless associated with measurable modifications of YFP fluorescence, allowing YFP fluorescence to be used to directly probe the dynamics of phenotypic diversity. By focussing on phenotypic variation across many generations of single cell lineages, our approach distinguishes variation in gene expression at the single cell level (noise) from heritable changes, such as those caused by mutations that are relevant for longer-term evolution. Evolutionary outcomes during long-term stationary phase have been shown to be sensitive to environmental conditions [49].
Here, we have demonstrated that a drop-based millifluidic system can be used to detect heritable phenotypic changes in evolving bacterial populations by tracking, in parallel, the lineage of descendants of single cells. By amplifying and averaging heritable changes over multiple generations variants, subtle phenotypic differences can be resolved in a way that is not possible using single-cell analysis techniques such as flow-cytometry. This approach offers a technological response to the growing need for reliable high-throughput phenotyping methods [50,51] to keep pace with current theoretical developments and advances in genomics [52]. Furthermore, since individual drops can be sorted into microtitre plates, and analyzed by barcoded whole genome next-generation sequencing, it is feasible to map genotype to phenotype. While the proof of principle work described here was based on the analysis of~500 single cells from each population, the technology is readily expanded to encompass the capture and analysis of the lineage of descendants of hundreds of thousands of cells.
The drop-based culturing strategy has utility in basic microbiology well beyond the simple application recounted here. It offers diverse applications in microbial physiology, ecology and evolution. For example, limiting dilutions of microbes from numerous sources can be readily partitioned among 1000s of drops (with the composition of drops being determined by the investigator), with growth from single cells being observed in real time, and independent populations then being available for genome sequencing, further phenotypic characterisation, and archiving. In addition, capacity to manipulate drops opens the door to a range of applications in experimental evolution and beyond.

Materials and Methods
Strains E. coli MC41000-YFP [17], contains yfp at the galK locus under control of the lac promoter. The lac operon (including lacI) is deleted from the MC4100 background, thus expression is not regulated by the lac repressor and is essentially constitutive, but the yfp locus remains under CRP control.

Culture conditions
All liquid cultures were in 5 mL of LB medium (Sigma-Aldrich) in closed polypropylene conical tubes (50 mL, VWR) and incubated at 37°C. First, 6 E. coli MC41000-YFP colonies were isolated on LB agar and then grown in liquid culture in closed tubes on an orbital shaker (VWR) at 100 rpm for five hours. Each of the 6 pre-cultures was used to inoculate 2x6 independent tubes (dilution 20x). Six tubes were incubated on an orbital shaker (VWR) at 100 rpm and six tubes were left static. Within each group of six tubes, the tubes were left in stationary phase for 1, 3, 6, 10, 15 and 30 days. At the end of the chosen duration each tube was mixed for homogeneous sampling. The samples were then stored at 80°C in 20% glycerol (VWR).

Millifluidic droplet analyser
The analysis with the millifluidic lineage tracking machine was performed after dilution in fresh LB medium (Sigma-Aldrich) to obtain, an average number of cell per drop λ 0.5. Drops were incubated and measured at 37°C as described in [12,13]. The fluorescence was measured using a photomultiplier (Hamamatsu) and a YFP filter set (479/40nm excitation and 530/40nm emission, Thorlabs) and an LED source (490nm, Thorlabs). In order to compare experiments, the data were normalized according to a reference measurement on the ancestral population.

Viable cell density determination
Viable cell density measured at day 1, 3, 6, 10, 15 and 30 was estimated using the millifluidic analyser. When generating drops, the cell encapsulation probability follows the Poisson law. If the average number of viable cells per drop is λ, the probability to form a drop with no bacteria is exp(−λ). This probability is estimated according to the fraction of drops with no detectable growth, x empty : λ*− ln(x empty ). Knowing the drop volume, we can finally determine the viable cells density within the sample: cfu = λ.V drop .

Colony morphology and Cytometry
Isolated clones (S3 Fig) of new phenotypes were plated on LB agar (BD). Colony morphologies were compared after 15h of incubation at 37°C. In order to test the relationship between fluorescence and the biomass, different isolated phenotypes were analysed by cytometry. Comparison of measurements as shown in Fig 3 were obtained for populations of each phenotype grown separately under the same conditions: 37°C shaking at 100 rpm (cultures were grown to stationary phase). Samples were then analyzed separately after dilution in phosphate buffered saline using a Partec cytometer (CyFlow1 Space).

Genome sequencing
The ancestral and 16 isolated clones (4 clones per phenotype) were grown overnight in 5 ml of M9. DNA was purified from these cultures using a bacterial genomic DNA kit (Sigma-Aldrich). DNA library construction and sequencing were carried out at the NGS platform of the Institut Curie. The 16 clones were DNA barcoded and paired-end sequenced (2 x 150 bp reads) in a single run on an Illumina HiSeq 2500. Mutations in the genomes were identified using the computational pipeline BRESEQ [53]. Accession number on SRA database at NCBI is SRX1564688.

Shannon index
The Shannon index [43] is defined as ∑ p i .ln(p i ) where p are the observed occurrences of phenotypic classes, i.e. p i ¼ n i = N , where n i is the observed number of the i-th phenotype and N the total number of observed cells.