Design and Experimental Evaluation of a Minimal, Innocuous Watermarking Strategy to Distinguish Near-Identical DNA and RNA Sequences.

The construction of powerful cell factories requires intensive and extensive remodelling of microbial genomes. Considering the rapidly increasing number of these synthetic biology endeavors, there is an increasing need for DNA watermarking strategies that enable the discrimination between synthetic and native gene copies. While it is well documented that codon usage can affect translation, and most likely mRNA stability in eukaryotes, remarkably few quantitative studies explore the impact of watermarking on transcription, protein expression, and physiology in the popular model and industrial yeast Saccharomyces cerevisiae. The present study, using S. cerevisiae as eukaryotic paradigm, designed, implemented, and experimentally validated a systematic strategy to watermark DNA with minimal alteration of yeast physiology. The 13 genes encoding proteins involved in the major pathway for sugar utilization (i.e., glycolysis and alcoholic fermentation) were simultaneously watermarked in a yeast strain using the previously published pathway swapping strategy. Carefully swapping codons of these naturally codon optimized, highly expressed genes, did not affect yeast physiology and did not alter transcript abundance, protein abundance, and protein activity besides a mild effect on Gpm1. The markerQuant bioinformatics method could reliably discriminate native from watermarked genes and transcripts. Furthermore, presence of watermarks enabled selective CRISPR/Cas genome editing, specifically targeting the native gene copy while leaving the synthetic, watermarked variant intact. This study offers a validated strategy to simply watermark genes in S. cerevisiae.

A DNA watermark is a unique synthetic nucleotide sequence that enables the identification and traceability of its carrier when applying PCR amplification and sequencing techniques. Application of the watermarks in living organisms started recently with a purpose to protect R&D investments, to create an information storage source, or to enable traceability of pathogenic or endangered species.[1−3] The literature reports successful embedding and subsequent detection of the watermarks in DNA strands in vitro,[4] as well as in vivo using several model microorganisms (i.e., Bacillus subtilis, Escherichia coli, Saccharomyces cerevisiae, Mycoplasma mycoides, and Mycoplasma capricolum), plants, and viruses.[1,2,5−11] All these studies focused on a single locus for the watermark introduction, with a few notable exceptions. First, the Mycoplasma genome de novo synthesis in which four large watermarks (ca. 1 kb) were introduced to enable the differentiation between the natural and synthetic copies of the Mycoplasma genome.[5] Second, the Synthetic Yeast 2.0 (Sc2.0) project, where approximately 28 bp regions of each open reading frame were recoded to distinguish synthetic from native genes by PCR.[12] Lastly, the recoding of the E. coli genome, such that it uses 61 instead of 64 codons.[10] The successes of these projects reveal the potential of the watermarks for future development in synthetic biology, particularly during large-scale genome remodeling projects, where tagging the synthetic gene copies can enable the discrimination between synthetic and native homologues. For instance, Kuijpers and co-workers recently reported the pathway swapping strategy that enables to redesign large, native essential pathways.[13] Pathway swapping was demonstrated on the glycolytic and fermentation pathways of S. cerevisiae, involving 12 catalytic steps encoded by 26 genes. After a first genetic reduction leading to a minimal glycolysis set of 13 genes,[14] a second, synthetic set of these 13 genes was integrated in a single locus on chromosome IX. Subsequently, the native copies of these 13 glycolytic genes were removed from their original chromosomal loci, leading to SwYG, a yeast strain with a single locus, minimal glycolytic pathway. However, the presence of two identical gene copies for all glycolytic genes during the strain construction process led to complications. First, in this intermediate strain carrying both native, scattered, and synthetic, colocalized glycolytic genes, removal of the native gene copies without harming the synthetic, identical copies integrated on chromosome IX was challenging. Second, expression of the native and synthetic genes could not be measured and compared. Both problems can easily be addressed by embedding watermarks in the synthetic genes. When judiciously placed in Protospacer Adjacent Motifs (PAM), watermarks can disable CRISPR/Cas editing in the synthetic genes.[15] When designed in coding regions (CDS), watermarks can be used to identify native from watermarked mRNA molecules.

Whether inserted in coding or noncoding regions, the major downside of watermarks is the risk of unintended changes in the host physiology. Watermarking in coding regions is potentially less challenging as watermarks can be embedded in the CDS as silent mutation, taking advantage of the redundancy of the genetic code encompassing 61 codons for only 20 amino acids. However, while “silent” or synonymous mutations in CDS do not affect the amino acid sequence of the corresponding protein, they can alter cells at different levels. Codons can be classified as optimal and nonoptimal based on their frequency in the genome and the abundance of tRNAs with complementary anticodons.[16−18] It is now well established that cells use codon optimality to tune protein expression. Highly expressed genes, such as genes encoding the highly abundant glycolytic proteins, are enriched for optimal codons.[19,20] Furthermore, by tuning the translation rate, codon optimality regulates the cotranslational folding of polypeptides and plays a role in shaping proteins conformational states.[21−24] More recently, it has been shown that codon optimality also modifies mRNA structure, splicing, and stability.[25−28] Codon optimality preservation is therefore an important criterion to consider when introducing watermarks without causing undesirable changes in gene function. There is however little known about the impact of watermarking on cell physiology, and remarkably few studies are dedicated to Saccharomyces cerevisiae, a microbe intensively used in synthetic biology developments.[29,30] Heider and Barnekow demonstrated that watermarking of VAM7 in S. cerevisiae did not affect the vacuolar function of the corresponding protein.[6] Liss and co-workers expressed a watermarked GFP in S. cerevisiae and showed minimal impact on GFP protein by Western blotting.[7] In the Sc2.0 project in every ORF larger than 500 bp at least two 19–28 bp PCRtags were introduced, which were recoded approximately 33–60%. Every strain with a native chromosome replaced by a synthetic version showed no or minor fitness defects, and transcript profiling showed only few genes changed in expression.[9,12,31−36] Whether these transcript changes originated from the PCRtags was not always investigated, and it is unclear whether these PCRtags allow discrimination between native and synthetic mRNAs when both are present in the cell. Therefore, there remains a strong need for studies proposing a watermarking strategy with the ability to distinguish between native and synthetic DNA and mRNA, validated by a systematic, quantitative exploration of the impact of watermarking on transcription, translation, and general physiology.[37]

To fill this knowledge gap, using S. cerevisiae as eukaryotic paradigm, this study designed, implemented and experimentally validated a systematic approach to watermark DNA with minimal alteration of yeast physiology. The impact of simultaneously watermarking 13 genes encoding abundant proteins involved in the major pathway for sugar utilization (i.e., glycolysis and alcoholic fermentation) on metabolism, transcriptome, and enzyme activity was explored using batch cultures in tightly controlled bioreactors. Watermarked transcripts were segregated from native ones using the karyollelle specific expression detection method.[38] Finally, the ability of watermarks to protect synthetic genes from CRISPR/Cas9 DNA editing was evaluated.

Results and Discussion

Design and In Silico Validation of the Watermarking Strategy

The presence of watermarks in the CDS of glycolytic genes shall enable discrimination of the watermarked versus native DNA and mRNA sequences with a minimal effect on transcript and protein levels, activity of enzymes in the glycolytic pathway and ultimately, yeast physiology. Finding the optimal trade-off between robust watermark detection by sequencing and minimal physiological impact was therefore the main design principle of the watermarking strategy. On the basis of current RNA sequencing resolution (Illumina platform with an error rate of <1%), at least five nucleotide substitutions were required to distinguish watermarked from native sequences using random single nucleotide polymorphisms (SNPs). Codon replacement was performed on the amino acids encoded by four to six alternative triplets (A, G, P, T, V, L, R, S), favoring triplets for which only the third base pair of the triplet was different from the original codon. The codon with the most similar percentage of abundance when referring to the codon usage table of S. cerevisiae (Table S1) was chosen, avoiding triplets leading to more than 20% variation in abundance when possible. The structure of the 5′ region of the mRNA is important for translation efficiency. Not only does the folding energy at the 5′ end affect translation initiation, but the presence of nonoptimal codons close to the initiation site can stall ribosomes, thereby hampering translation initiation.[39,40] Furthermore, as translation initiation is considered as translation limiting step,[41] and following the example of Annaluru and colleagues,[9] the first 101 nt of the CDS were preserved. The optimal distribution of watermarks over the remaining CDS stretch was tested with two in silico approaches using ADH1. In the first approach, watermarks were colocalized in two 100 nt regions, one in the middle of the CDS, and the other located 10 nt upstream of the stop codon (Figure A). In the second approach, base pair substitutions were equally spread over the CDS sequence, every 85 nt (Figure B). In both approaches 11 watermarks were introduced, which resulted in an overall change in codon usage of 0.57 for the first and 0.47 for the second method (Figure , Figure S1. See Methods section for calculation of the change in codon usage). In order to evaluate the discriminatory potential of RNA sequencing with these two strategies, 100 bp paired-end sequence reads were simulated for both watermarked and native ADH1 copies (see Methods section). These data were processed using the k-mer method developed by Gehrmann et al.,[38] only considering reads containing watermarks (see Methods section), to selectively quantify watermarked and native reads. On average 52.5% of the reads were captured when using the first approach with clustered watermarks, while 99.4% were detected using the second approach, with watermarks spread over the CDS (Figure ). The Pearson correlation coefficient between generated and measured reads was above 0.99 for both methods, indicating that both methods are able to retrieve the variation in abundance across the samples, required for differential expression. The second approach resulted in a better sequence coverage and slightly lower codon usage change. However, the first approach is less labor intensive when manual design is performed, and is less likely to affect cotranslational folding,[21] as a shorter part of the CDS undergoes codon usage change. The first method was therefore selected as watermarking strategy (detailed in Box 1 and Figure S2) and used to edit in silico the CDS of 13 genes of glycolysis and alcoholic fermentation (HXK2, PGI1, PFK1, PFK2, FBA1, TPI1, TDH3, PGK1, GPM1, ENO2, PYK1, PDC1, and ADH1, see example for watermarking for FBA1 and ENO2 in Figure S3). This resulted in a reasonably low change in codon usage of the watermarked genes of 0.5 on average (Table S2). Using simulated data, we compared the performance of the k-mer method with traditional alignment and found that the k-mer method was able to achieve a higher read retrieval rate than alignment indicating a more accurate transcription estimate (Table S3). The watermarked CDS were synthesized with flanks compatible with Golden Gate assembly (plasmids pGGKp137 to pGGKp150, Table S4).

Figure 1

Comparison of two watermarking strategies. (A) First strategy with clustered watermarks. (B) Second strategy with watermarks distributed over the whole coding region. The tables in panel A and B represent the % of sequencing reads that can be captured by the two watermarking strategies, calculated from in silico simulated 100 bp paired-end sequencing reads.

Strain Construction Strategy and Confirmation

In the SwYG strain,[13] the set of genes involved in glycolysis and fermentation was reduced from 26 to 13 and relocalized to a single locus (Single Locus Glycolysis, SinLoG) on chromosome IX. The SwYG strain is a perfect platform to rapidly remodel glycolysis and alcoholic fermentation and test multiple (heterologous) variants. SwYG was therefore used as starting strain to express the watermarked genes. Using simultaneous Cas9-mediated genome editing and in vivo assembly, the entire glycolytic and fermentation pathways composed of 13 watermarked genes were integrated in one step in the CAN1 locus on chromosome V. The watermarked genes were framed by the native, standardized corresponding promoters and terminators (800 bp and 300 bp respectively). Three helper elements, two Autonomously Replicating Sequences (ARS) and a selection marker were included in the SinLoGs design (Figure ). Two active ARSs (ARS418 and ARS1211) were added on both ends of the ca. 35 Kb long SinLoGs to minimize the risk of perturbing DNA replication of this long DNA stretch. A selection marker was used to facilitate screening for correct integration and removal of the SinLoGs. The native SinLoG, present in the SGA1 locus on chromosome IX, was then removed using the Cas9 endonuclease, resulting in strain IMX1770 (Figure and Figure ). To obtain an isogenic control strain, the same procedure was followed to construct a strain with native SinLoG, framed by the same promoters and terminators as the watermarked genes, and integrated in the same CAN1 locus on chromosome V (Figure and 3). This control strain was named IMX1771. The genome of both strains was sequenced, confirming the presence of a single, correctly assembled glycolytic pathway at the targeted chromosomal location. Sequencing revealed the absence of mutations in the coding regions of the glycolytic and marker genes but identified a few mutations in the promoter and terminator regions of the glycolytic and selection marker expression cassettes (Table S5). In IMX1770, a single Single Nucleotide Variation (SNV) was found in the promoter of PFK1 and HIS3 and in the terminator of PGK1 and ENO2. In IMX1771, a single SNV was identified in the HIS3 terminator, and a short TA stretch was missing in the promoter of GPM1.

Figure 2

Construction of SinLoG (Single Locus Glycolysis) strains IMX1770 and IMX1771 using the glycolysis swapping strategy.[13] (A) A newly designed glycolysis is integrated in the CAN1 locus by simultaneous CRISPR/Cas9-aided editing of CAN1 and in vivo assembly of glycolytic expression cassettes and helper fragments (ARS418, ARS1211 and the selection marker HIS3). The > and < signs next to the gene names indicate the directionality of transcription and letters indicate the synthetic homologous recombination (SHR) sequence which was used for assembly. (B) Subsequently, the Single Locus Glycolysis present in the SGA1 locus was excised by double editing using CRISPR/Cas9 and replaced by the URA3 selection marker. The set of genes integrated in CAN1 is then the sole set of glycolytic genes present in the newly constructed strain and is essential for growth on glucose.

Figure 3

Strain construction workflow. The Switchable Yeast Glycolysis (SwYG) strain, IMX1338, served as parental strain to introduce in chromosome V a SinLoG (Single Locus Glycolysis) with native ORFs and standardized promoters/terminators (IMX1747) as well as with watermarked ORFs and standardized promoters/terminators (IMX1717). From both strains the native SinLoG in chromosome IX with variable promoters and terminators was removed (resulting in strain IMX1771 and IMX1770, respectively). After removal of URA3 from strain IMX1771 (native ORFs) the SinLoG with watermarked ORFs was introduced in chromosome IX, resulting in a strain with double glycolysis (IMX2028).

A third strain, IMX2028 was constructed. IMX2028 harbored a double SinLoG, one located on chromosome IX carrying the watermarked genes and another on one carrying the native yeast genes on chromosome V (IMX2028, Figure and Table ). Unfortunately, genome sequencing revealed the deletion of a large region of the mitochondrial DNA (Figure S4). The strain IMX2028 was constructed to evaluate the performance of the k-mer method to discriminate between watermarked and native genes when present in the same strain. Despite IMX2028 respiration deficiency, the watermarked and native SinLoG that this strain carried were essentially faithful to the in silico design (Table S5), which made this strain still valuable for differential quantification of watermarked and native genes.

Table 1

List of Strains Used in This Study

strain name	genotype	source and brief description
Strains Characterized in This Study
IMX1770	MATaura3–52his3–1leu2–3,112MAL2–8cSUC2hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4pyk1pgi1tpi1tdh3pfk2::(pTEF-cas9-tCYC1 natNT1)pgk1gpm1fba1hxk2pfk1adh1pdc1eno2glk1::Sphis5Δ::(pGAL1-I-SceI-tCYC1)can1::(ARS418CAN1 AHFBA1_*AH HTPI1_*HPPGK1_*PQADH1_*QNPYK1_*NOTDH3_*OAENO2_*ABHXK2_*BCPGI1_*CDPFK1_*DJPFK2_*J BPHIS3BP LGPM1_*LMPDC1_*M ARARS1211AR CAN1) sga1::KlURA3	This study; Prototrophic strain with watermarked single locus glycolysis (WMG strain); Derived from IMX1717
IMX1771	MATaura3–52his3–1leu2–3,112MAL2–8cSUC2hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4pyk1pgi1tpi1tdh3pfk2::(pTEF-cas9-tCYC1 natNT1)pgk1gpm1fba1hxk2pfk1adh1pdc1eno2glk1::Sphis5Δ::(pGAL1-I-SceI-tCYC1)can1::(ARS418CAN1 AHFBA1AH HTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2J BPHIS3BP LGPM1LMPDC1M ARARS1211AR CAN1)sga1::KlURA3	This study; Prototrophic strain with native single locus glycolysis (NG strain); Derived from IMX1747
IMX2028	MATa ura3–52his3–1leu2–3,112MAL2–8cSUC2hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4glk1::Sphis5Δ::(pGAL1-I SceI-tCYC1)can1::(ARS418CAN1 AHFBA1AH HTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2J BPHIS3BP LGPM1LMPDC1M ARARS1211AR CAN1)sga1::(ARS418sga1 AHFBA1_*AH HTPI1_*HPPGK1_*PQADH1_*QNPYK1_*NOTDH3_*OAENO2_*ABHXK2_*BCPGI1_*CDPFK1_*DJPFK2_*J BPHIS3BP LGPM1_*LMPDC1_*M ARARS1211AR sga1)	This study; Prototrophic strain with native and watermarked single locus glycolysis (DG strain); Derived from IMX1748
Strains Used As Starting Point or Intermediate in the Construction of the above Strains
CEN.PK113–7D	MATaMAL2–8cSUC2	Control strain;[46,47] Starting strain for all construction work;[13] Contains a SinLoG in Chr. IX, with variable promoter and terminator length; Uracil auxotroph
IMX589	MATaura3–52his3–1leu2–3,112MAL2–8cSUC2glk1::Sphis5hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4sga1::(FBA1GHTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2JKAmdSYMKLGPM1LMPDC1-SYNMF)pyk1pgi1tpi1tdh3pfk2::(pTEF-cas9-tCYC1 natNT1)pgk1gpm1fba1hxk2pfk1adh1pdc1eno2
IMX1338	MATaura3–52his3–1leu2–3,112MAL2–8cSUC2hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4sga1::(FBA1GHTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2JKAmdSYMKLGPM1LMPDC1-SYNMF)pyk1pgi1tpi1tdh3pfk2::(pTEF-cas9-tCYC1 natNT1)pgk1gpm1fba1hxk2pfk1adh1pdc1eno2glk1::Sphis5Δ::(pGAL1-ISceI-tCYC1)	This study; Derived from IMX589; Contains a SinLoG in Chr. IX; Uracil and histidine auxotroph
IMX1717	MATaura3–52his3–1leu2–3,112MAL2–8cSUC2hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4sga1::(FBA1GHTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2JKAmdSYMKLGPM1LMPDC1-SYNMF)pgk1gpm1fba1hxk2pfk1adh1pdc1eno2glk1::Sphis5Δ::(pGAL1-I SceI-tCYC1)can1::(ARS418AHFBA1_*AH.HTPI1_*HPPGK1_*PQADH1_*QNPYK1_*NOTDH3_*OAENO2_*ABHXK2_*BCPGI1_*CDPFK1_*DJPFK2_*J.BPHIS3BP.LGPM1_*LMPDC1_*M.ARARS1211AR)	This study; Derived from IMX1338; Strain with native (variable prom and term length) and watermarked SinLoG in Chr. IX and V respectively; Histidine auxotroph
IMX1747	MATaura3–52his3–1leu2–3,112MAL2–8cSUC2hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4 sga1::(FBA1GHTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2JKAmdSYMKLGPM1LMPDC1-SYNMF)pyk1pgi1tpi1tdh3pfk2::(pTEF-cas9-tCYC1 natNT1)pgk1gpm1fba1hxk2pfk1adh1pdc1eno2glk1::Sphis5Δ::(pGAL1-I SceI-tCYC1)can1::(pGAL1-I-SceI-tCYC1)can1::(ARS418CAN1 AHFBA1AH HTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2J BPHIS3BP LGPM1LMPDC1M ARARS1211AR CAN1)	This study; Derived from IMX1338; Strain with native SinLoG with variable prom and term length in Chr. IX and native SinLoG with standardized prom and term length in Chr. V; Histidine auxotroph
IMX1748	MATaura3–52his3–1leu2–3,112MAL2–8cSUC2hxk1::KlLEU2tdh1tdh2gpm2gpm3eno1pyk2pdc5pdc6adh2adh5adh4pyk1pgi1tpi1tdh3pfk2::(pTEF-cas9-tCYC1 natNT1)pgk1gpm1fba1hxk2pfk1adh1pdc1eno2glk1::Sphis5Δ::(pGAL1-I-SceI-tCYC1)can1::(ARS418CAN1 AHFBA1AH HTPI1HPPGK1PQADH1QNPYK1NOTDH3OAENO2ABHXK2BCPGI1CDPFK1DJPFK2J BPHIS3BP LGPM1LMPDC1M ARARS1211AR CAN1)	This study; Derived from IMX1771; Strain with watermarked SinLoG in Chr. IX and native SinLoG with standardized prom and term length in Chr. V; Uracil auxotroph

Watermarks Do Not Affect Yeast Physiology

To evaluate the impact of DNA watermarking on yeast physiology, the watermarked strain IMX1770 and its isogenic control IMX1771 were grown in aerobic batch bioreactors and their growth kinetics were compared. Both strains were prototrophic, meaning that they fully relied on glucose, the sole carbon and energy source catabolized via glycolysis, to produce the required cellular building blocks and therefore to grow. The two strains displayed identical growth rates (0.33 ± 0.004 h–1 and 0.32 h–1 ± 0.002 h–1 for IMX1770 and IMX1771, respectively) as well as glucose and O2 uptake rates, ethanol and CO2 production rates and yields (Figure and Table ). Both strains passed the diauxic shift and grew equally well using the ethanol, which was produced during fermentative growth on glucose, as carbon and energy source (Figure S5). Watermarking of glycolytic and alcoholic fermentation genes therefore did not alter metabolic fluxes and the overall physiological responses during fast respiro-fermentative on glucose and full respiratory growth on ethanol.

Figure 4

Physiological characterization of strains with native (IMX1771) and watermarked (IMX1770) glycolytic genes during aerobic batch cultures in bioreactors. (A) Left panel, biomass concentration (gram dry biomass per liter), central panel, glucose concentration (mM); right panel, ethanol concentration (mM). Three independent culture replicates are represented for each strain. Shades of blue with square symbols, IMX1770, shades of red with round symbols, IMX1771. (B) Specific enzyme activities of the 12 reactions encoded by the 13 glycolytic enzymes (Pfk1 and Pfk2 form an enzyme complex) of the strains with native (IMX1771, red bars) and watermarked glycolysis (IMX1770, blue bars). Samples were taken in midexponential phase. Bars represent the average and standard deviation of measurements from three independent batch cultures for each strain. Stars indicate enzyme activities that are significantly different between the two strains (Student t test, p-value threshold 0.05, two-tailed test, homoscedastic).

Table 2

Physiological Characterization in Bioreactor of Yeast Strains with Native and Watermarked Glycolysisa

	IMX1771	IMX1770
Yields
Ysx (gdry weight/gglucose)	0.12 ± 0.01	0.13 ± 0.01
Ys,glycerol (mol/mol)	0.07 ± 0.01	0.05 ± 0.01
Ys,ethanol(mol/mol)	1.47 ± 0.05	1.46 ± 0.02
Ys,acetate(mol/mol)	0.05 ± 0.00	0.06 ± 0.00
Specific rates
μmax (h–1)	0.32 ± 0.00	0.33 ± 0.01
qglucose (mmol/g–1·h–1)	–14.6 ± 0.7	–14.7 ± 0.9
qglycerol (mmol/g–1·h–1)	0.96 ± 0.03	0.76 ± 0.17
qethanol(mmol/g–1·h–1)	21.4 ± 0.3	21.4 ± 1.7
qacetate(mmol/g–1·h–1)	0.79 ± 0.05	0.88 ± 0.06
Carbon balances (%)	105 ± 2	103 ± 3

Physiological characterization of strains with native (IMX1771) and watermarked (IMX1770) glycolytic genes during aerobic batch cultures in bioreactors. Data represent the average and standard deviation of measurements from three independent batch cultures for each strain. Statistical analysis (Student t test, p-value threshold 0.05, two-tailed test, homoscedastic) revealed no significant differences between the two strains.

Watermarking might affect protein folding and consequently function. However, as yeast glycolysis is characterized by an overcapacity of its enzymes, mild variations of glycolytic enzymes activities might not be detectable by growth kinetics. The 12 specific activity assays of the 13 enzymes (Pfk1 and Pfk2 are subunits of a hetero-octameric phosphofructokinase[42]) encoded by the watermarked genes were therefore assayed in vitro. The specific activity of these 13 enzymes was, with the exception of Gpm1, remarkably similar between watermarked and native strains (p-values above 0.05; Student t test, two-tailed, homoscedastic). For all enzymes, specific activities were remarkably similar to protein abundance, including a 1.6-fold decrease in specific activity and protein abundance for Gpm1 (Figure , Figure S8 and S9). Watermarking therefore did not affect or marginally affect protein expression and functionality (Figure ).

To further explore the potential impact of watermarking on yeast physiology, the transcriptome of IMX1771 and IMX1770 grown in aerobic batch reactors was compared. The transcriptional response of these two strains was remarkably similar (Figure ). The native and watermarked glycolytic genes were the only differentially expressed genes between the two strains. This differential expression reflects the absence of the native genes and therefore their lack of expression in IMX1770, and the absence and lack of expression of the watermarked genes in IMX1771. However, expression levels of the native and watermarked genes in IMX1771 and IMX1770, respectively, were highly similar (Figure A).

Figure 5

Genome-wide transcriptome analysis of IMX1770 and IMX1771. The x-axis represents the log fold change in expression, and the y-axis represents the −log p-value. Each point represents a transcript. A negative log fold change reflects higher expression in the native strain than in the watermarked strain, and vice versa. The horizontal, dashed line represents the FDR corrected p-value threshold of 0.05, and the vertical dashed lines represent a log fold change threshold of 1. Red points indicate significantly differentially expressed transcripts (FDR-corrected p-value above 0.05 and Log fold change higher than 1).

Figure 6

Glycolysis and fermentation transcript levels of S. cerevisiae grown in aerobic batch cultures in bioreactors. (A) Watermarked transcript levels of IMX1770 (blue) and native transcript levels of IMX1771 (red). (B) Watermarked and native transcript levels of IMX2028. Bars represent the average and standard deviation of three independent cultures replicates. Samples were taken in midexponential phase (Table S9). No significant change in expression was found between watermarked and native genes (Student t test, p-value threshold 0.05, two-tailed test, homoscedastic) between IMX1770 and IMX1771 and within IMX2028.

Physiological characterization of IMX2028 confirmed the respiration deficiency suggested by the absence of mitochondrial DNA. The k-mer method was able to selectively quantify expression of the native and watermarked genes. While expression of glycolytic and respiration genes might differ in IMX2028 as compared to IMX1770 and IMX1771 due to the mutations in mitochondrial DNA and associated respiration deficiency, the relative expression of glycolytic and fermentation genes, expressed from the same promoters in the native and watermarked SinLoG, was not expected to differ between the native and watermarked genes in this strain. Accordingly, and in agreement with the similarity of the expression levels between IMX1770 and IMX1771, transcript levels of native and watermarked genes in IMX2028 were identical (Figure B).

Watermarking of 13 highly expressed genes of central carbon metabolism, essential for glucose utilization, had therefore no impact on yeast transcriptome and physiology.

Watermarking Enables Selective CRISPR/Cas9 Genome Editing

DNA binding and editing by CRISPR/Cas9 requires the presence at the targeted site of a specific PAM recognition sequence.[15] A single nucleotide variation in this sequence can abolish Cas9 ability to introduce a double strand DNA break.[15] This feature is particularly interesting when considering selective editing of identical or highly similar sequences. If strategically designed, watermarks can enable targeted editing of a watermarked gene, leaving the native copy intact or conversely, prevent editing of the watermarked gene while cutting the native copy. Guide RNAs (gRNAs) selectively targeting the native copies of PYK1 and TPI1 for CRISPR/Cas9 editing were designed (Table S6), inserted into expression vectors and transformed to IMX1717, a double SinLoG strain and direct ancestor of IMX1770 (Figure ). Double-stranded DNA fragments of 120 nt were supplied during transformation to repair via homologous recombination the break induced by CRISPR/Cas9. As the sequence of the native and watermarked genes is identical with the exception of watermarks, a single primer set designed just outside the open reading frame can be used to amplify both copies of PYK1 or TPI1 in a single PCR reaction using IMX1717 genomic DNA as template. Ran on a gel, the PCR products of this reaction would lead to a single band corresponding to both the native and watermarked copies of PYK1 or TPI1. Selective editing would lead to the appearance of a second, smaller band on gel, corresponding to the edited copy of PYK1 or TPI1. Out of 15 colonies of IMX1717 transformed with the gRNA targeting PYK1, three displayed two bands demonstrating editing of a single PYK1 copy (Figure ). Five out of 15 colonies of IMX1717 transformed with the gRNA targeting TPI1 showed selective editing (Figure ). For two transformants per gene (TPI1 and PYK1) showing two bands on the gel, sequencing the largest band confirmed the presence of the watermarked sequence only, confirming selective editing of the native TPI1 and PYK1 (Figure S6). It has recently been shown that cells can use chromosomal DNA with high homology to repair a CRISPR/Cas mediated DNA break, leading to loss of heterozygosity.[43] In the present case, it means that cells could repair the induced DNA break in the targeted, native gene copy with its watermarked homologue, resulting in two copies of the watermarked gene, but a single PCR product and therefore a single band on gel. Sequencing of the unique PCR product of four colonies in which editing of PYK1 and TPI1 was considered unsuccessful revealed that, for all tested colonies, the PYK1 and TPI1 genes were cut by CRISPR/Cas but repaired by (part of) the watermarked allele. Editing of the targeted, native genes by CRISPR/Cas9 was therefore highly efficient (100% of the tested colonies), however the DNA break was repaired either by the supplied repair DNA fragment or by the watermarked homologue.

Figure 7

Diagnostic PCR for selective editing of native glycolytic genes. Separation of PCR products resulting from outside–outside amplification to identify edited (nonwatermarked) and nonedited (watermarked) loci for PYK1 (A) and TPI1 (B) from transformants of IMX1717 (double SinLoG). (A) Lanes 1–15 show the PCR results of amplification of the PYK1 locus of randomly picked colonies. Successful editing of the locus results in a DNA fragment with a length of 670 bp. No editing of the locus results in a DNA fragment with a length of 2177 bp. Primers 11915 and 4667 were used. Lanes 1, 5, and 15 display bands of both sizes revealing selective editing. (B) Lane 1–15 show the PCR results of amplification of the TPI1 locus of randomly picked colonies. Successful editing of the locus results in a DNA fragment with a length of 378 bp. No editing of the locus results in a DNA fragment with a length of 1125 bp. Primers 3514 and 6406 were used. Lanes 9–11, 13, and 15 display bands of both sizes revealing selective editing. A negative control is indicated with “C-“ (IMX1338, SinLog). In the lanes indicated with “L”, GeneRuler DNA ladder mix was loaded. 1% (w/v) agarose in TAE.

Conclusion

The present study offers an innocuous watermarking strategy for coding regions that enables the discrimination of DNA and mRNA by sequencing through a k-mer approach and facilitates selective editing of watermarked and nonwatermarked sequences. While the design of watermarked genes was performed manually in the present study, it can easily be automated when a larger number of genes is concerned, with software similar to for example GeneDesign.[44]

The set of genes chosen to test the watermarking strategy encodes highly abundant proteins that are generally considered to operate at overcapacity, which means that the capacity of the enzymes is considerably larger than the actual in vivo flux. While this overcapacity might obscure physiological responses, a closer inspection of transcript levels and enzyme activities confirmed the watermarks' harmlessness for S. cerevisiae. Only one of the 13 tested genes showed an activity of the enzyme encoded by the watermarked allele significantly decreased (ca. 1.6-fold) as compared to the activity of the enzyme encoded by the native allele (phosphoglucomutase encoded by GPM1), which could be explained by a similar decrease in Gpm1 protein abundance. As native and watermarked transcript levels were identical for GPM1, the lower enzyme abundance in the watermarked strain might result from a slightly lower translation efficiency. Neither the watermarking specifics (type of codon substitution, change in codon usage, etc.) of this particular protein, nor information from literature hinted toward the mechanism underlying this decreased protein abundance. A recent study, combining measurements of protein synthesis rate with ribosome footprinting data confirmed that Gpm1, like most glycolytic proteins, has a fast synthesis rate.[45] Applying the same approach to watermarked strains could help characterizing the impact of nucleotide substitution on translational efficiency.

Another particularity of glycolytic genes is their high codon optimality (on average ca. 90% of optimal codons according to Hanson and Coller[24]). While one could argue that this set of genes might not be representative of the yeast genome, to the best of our knowledge there is no evidence that such genes are more or less robust toward changes in codon frequency than genes with lower levels of codon optimality. For future studies it would be interesting to explore if codon optimality affects genes sensitivity to watermarking.

Methods

Strains and Cultivation Conditions

The S. cerevisiae strains used in the study belong to CEN.PK family[46−48] and are listed in Table . Liquid cultures were grown in 500 mL shake flasks filled with 100 mL of medium at 30 °C with 200 rpm agitation. Complex media (further referred to as YPD) contained 10 g·L–1 yeast extract, 20 g·L–1 peptone and 20 g·L–1 glucose. Synthetic minimal medium (further referred as SMG) consisted of 3 g g·L–1 KH2PO4, 0.5 g·L–1 MgSO4·7H2O, 5 g·L–1 (NH4)2SO4, 1 mL·L–1 of a trace element solution, and 1 mL·L–1 of a vitamin solution as previously described[13] and supplemented with 20 g·L–1 glucose. For solid medium, 20 g·L–1 agar was added prior autoclaving. When selection in SMG was required, (NH4)2SO4 was replaced with 3 g·L–1 K2SO4 and 2.3 g·L–1 filter-sterilized urea, and the medium was supplemented with 200 mg·L–1 of G418, hygromycin B or 10 mM acetamide.[49,50] For the counterselection purpose, 1 mg·mL–1 5-FOA (Zymo Research, Irvine, US) was added to SMG supplemented with uracil (150 mg·L–1).[49] For plasmid propagation, E. coli XL1-Blue cells (Agilent Technologies, Santa Clara, CA, USA) were grown in Lysogeny broth (LB) medium supplied with 100 mg Lampicillin or 25 mg·L–1 chloramphenicol at 37 °C with 180 rpm agitation. Yeast and bacterial frozen stocks were prepared by addition of 30% (v/v) glycerol to exponentially growing cultures. Strain aliquots were stored at −80 °C.

Molecular Biology Techniques

PCR reactions for diagnostic purposes were performed using DreamTaq DNA polymerase Master Mix (Thermo Fisher Scientific, Walthman, MA, USA) according to the manufacturer’s instructions. For high fidelity PCR reactions, Phusion High-Fidelity DNA polymerase (Thermo Fisher Scientific) was used following the supplier’s manual. Oligonucleotides of desalted or PAGE quality, depending on the purpose, were purchased from Sigma-Aldrich (St Louis, MO, USA). DNA fragments were resolved in agarose gels and purified using PCR cleanup kit from the reaction mixture (Sigma-Aldrich, St Louis, MO, USA) or excised from the agarose gel and purified using Zymoclean gel purification kit (Zymo Research, Irvine, CA, USA) when required. Circular templates were removed by applying DpnI enzyme restriction according to the producer’s manual (Thermo Fisher Scientific).

Plasmids were isolated from E. coli cultures using Sigma GenElute Plasmid kit (Sigma-Aldrich, St Louis, MO, USA). E. coli transformations were performed using chemical competent XL-1 Blue cells (Agilent Technologies, Santa Clara, CA, USA) according to the manufacturer instructions. Golden Gate Assembly was performed as previously described[51] using equimolar concentrations of 20 fmol for each fragment. For a 10 μL reaction mixture 1 μL T4 DNA ligase buffer (Thermo Fisher Scientific), 0.5 μL T7 DNA ligase (NEB New England Biolabs, Ipswich, MA), and 0.5 μL of either FastDigest Eco31I (BsaI) or BsmBI (NEB) were added.

Gibson Assembly was performed using Gibson Assembly Master Mix (New England Biolabs, Ipswich, MA) according to the manufacturer’s protocol.

All plasmids are reported in Table S4 and primers in Table S7.

In Silico Design of the Watermarks

Watermarks were introduced in the genes of interest according to the guidelines described in Box 1 using the Clone Manager software.

Watermarked CDS were ordered as a synthetic gene from GenArt (Thermo Fisher, Regensburg). The list of synthesized plasmids encoding watermarked CDS (pGGKp137 to pGGKp150) can be found in Table S4.

The change in codon usage in a gene caused by watermarking was calculated aswhere i represents each codon substitution in a gene.

Construction of Libraries Encoding Transcriptional Units of Watermarked Glycolysis

The sequences of the watermarked genes, promoters (800 bp) and terminators (300 bp) were ordered from GeneArt (Thermo Fisher, Regensburg, Germany). For compatibility with Golden Gate Cloning, the sequences were ordered flanked with BsaI and BsmBI restriction sites. The promoters and terminators were delivered by GeneArt subcloned in the entry vector pUD565 and for the watermarked genes the subcloning into pUD565 was done in house using BsmBI Golden Gate cloning. An exception was made for pTDH3, pPGK1, tPGK1, tENO2, and tADH1 which were amplified from genomic DNA of CEN.PK113–7D using primers with flanks containing BsaI restriction sites listed in Table S7.

Subsequently, the assembly of the promoter, gene and terminator was done in the preassembled vector pGGKd012 using Golden Gate cloning as described in the previous section. pGGKd012 was assembled from the Yeast toolkit[51] plasmids pYTK-002, 047, 072, 078, 081, and 083 (Table S4).

Correct plasmid assembly was verified by enzyme digestion with either BsaI, BsmBI (New England Biolabs) or FastDigest enzymes (Thermo Fisher Scientific) following the manufacturer’s instructions.

Construction of gRNA Plasmids Used in the Study

The guide RNA (gRNA) plasmids pUDR413 and pUDR529 for the yeast strain construction were designed and constructed according to Mans et al. (2015).[52] gRNA targets were selected using the Yeastriction tool[52] in case of pUDR413, or designed manually for the K. lactis URA3 target in plasmid pUDR529. For pUDR413, the 2 μm fragment was constructed in two parts using the primer 6131 and 5975 and primer 6296 together with 5941 using pROS12 as a template. For pUDR529, the 2 μm fragment was obtained by PCR using primer 14549 and pROS12 as a template. The backbone for pUDR413 was amplified with primers 6005 and 6006 using pROS13 as a template, while for pUDR529 same primer pair was used to amplify the backbone from pROS12. For both plasmids, 100 ng of each purified fragment was used in the Gibson Assembly and correct plasmid assembly was verified with the primers 3841, 5941, and 6070 in case of pUDR413 and 4034 and 5941 for pUDR529.

The guide RNA plasmids for selective native copy gene removal of TPI1 and PYK1, named pUDR531 and pUDR532 respectively (Table S4), were constructed as described in Mans et al. (2015)[52] with the modifications regarding the design of the gRNA. gRNAs were designed manually to target the native CDS containing a PAM which was removed in the watermarked copy of the CDS. Each gRNA was ordered as a primer (Table S7, primers 14515, 14517, 14519, 14521). The 2 μm fragment for four gRNA plasmids was obtained by PCR using corresponding gRNA primer (Table S7, primers 14515, 14517, 14519, 14521) and pROS13 as a template. The backbone for the four plasmids was obtained by amplification with primers 6005 and 6006 using pROS12 as a template. For the assembly, 100 ng of purified backbone and gRNA fragments were used in the Gibson Assembly and correct plasmid assemblies were verified with the primers 3841 and 5941 in combination with gRNA specific primers listed in Table S7.

Construction of SwYG Strains with Native and Watermarked Glycolysis and Double Glycolysis Strain

A schematic overview of the strain construction approach is shown in Figure . All yeast transformations were performed according to Gietz and Woods (2002).[53] For highly efficient targeted integration CRISPR/Cas9 mediated editing was applied. To this end, 350 ng of a plasmid carrying a corresponding guide RNA (further gRNA) was transformed into the yeast strain together with a purified PCR fragment (150 fmol) containing 60 bp homology to the integration site and acting as donor DNA[4] (see primers listed in Table S7). gRNA plasmids and the donor DNA were specific for each strain construction step and will be specified below. When donor DNA was consisting of multiple fragments, 60 bp sequences for homologous recombination (SHR) were flanking each of the fragments to enable in vivo assembly by homologous recombination. PCR fragments for the native SinLoG genes and for ARS418, ARS1211 and HIS3 were obtained using CEN.PK113–7D genomic DNA as a template, while fragments encoding the watermarked SinLoG were amplified from plasmids encoding the corresponding transcriptional units (Table S4, Table S7).

To obtain a double auxotrophic host strain named IMX1338, the Schizosaccharomyces pombe HIS5 gene previously inserted in the glk1 locus of IMX589 containing the SinLoG in chromosome IX[13] was replaced by the I-SceI expression cassette (pGAL1 – I-SceI – tADH1), which was amplified from the plasmid pUDC073 (primers 10708 and 10709). The replacement was mediated by a Cas9 gRNA plasmid assembled in vivo from two PCR fragments amplified from the pMEL10 plasmid using primers 6005 and 6006 in combination with 10904 (gRNA primer). Transformants were selected on SMG, and the gRNA plasmid with URA3 marker was removed by two sequential restreaks on SMG with 5-FOA. The correct genotype was confirmed by diagnostic PCR using primers 6190 + 1525 and 1553 + 6189 and later by whole genome sequencing. To construct IMX1717 and IMX1747, IMX1338 was transformed with the p426-SNR52p-gRNA.CAN1.Y-SUP4t plasmid targeting the CAN1 locus,[54] and PCR fragments of the 13 native or watermarked SinLoG glycolytic genes together with ARS418, ARS1211, and the HIS3 marker gene (Figure ). Transformants were selected on SMG media and after strain confirmation by PCR (Table S7, Figure S7) the gRNA plasmid encoding the KlURA3 marker was removed. As the next step, the SinLoG with variable length of promoters and terminators was removed from the SGA1 locus in the strains IMX1717 and IMX1747. To this end, both strains were transformed with plasmid pUDR413 and 1 μg of KlURA3 repair fragment amplified with primers 13273 and 13274 introducing homology flanks to the SGA1 site (Figure ). Transformants were selected on SMG supplemented with G418 and after strain confirmation by PCR using primers 11898 + 7479, 11898 + 2363, and 170 + 7479, the plasmid was removed. For the construction of the strain IMX2028 containing the native SinLoG in Ch V and watermarked SinLoG in Ch IX, first, the intermediate strain IMX1748 was constructed by removing the KlURA3 gene from IMX1771. This was done by transformation with plasmid pUDR529 encoding a gRNA for the KlURA3 gene and a repair fragment amplified with the primers 4223 and 4224 and containing homology to the SGA1 locus. Colonies were selected on YPD media supplemented with Hygromycin B and correct strain construction was confirmed by PCR using primers 4223 and 4224. After KlURA3 marker removal, IMX1748 was transformed with plasmid pUDR314 and the mixture of fragments for the watermarked SinLoG, ARS418, ARS1211, and KlURA3 marker gene resulting in strain IMX2028. Correct assembly of the fragments was confirmed by PCR.

Selective CRISPR/Cas9 Genome Editing

For selective CRISR/Cas9 genome editing, IMX1717 (Table ) was transformed with 1 μg of a 120 bp repair fragment with homology to the beginning and end of the gene and with 1 μg of plasmids pUDR531 or pUDR532 containing a gRNA for TPI1 and PYK1 respectively as described in the section Construction of gRNA plasmids and Table S4. Cells were plated on YPD with Hygromycin B. Repair fragments (120 nt-long) and diagnostic primers are listed in Table S7-G.

Whole Genome Sequencing and Data Analysis

Yeast genomic DNA was isolated using the Qiagen Genomic DNA Buffer Set and Genomic-tip 100/G tips (Qiagen, Hilden, Germany) according to the manufacturer’s manual. The incubation step with zymolyase was performed for 11 h and the incubation step for digestion with proteinase K was performed overnight. The concentration of the genomic DNA mixture was measured with the BR ds DNA kit (Invitrogen, Carlsbad, CA, USA) using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and the purity was verified with a Nanodrop 2000 UV–vis Spectrophotometer (Thermo Fisher Scientific).

IMX1770, IMX1771, and IMX2028 genomes were sequenced on an Illumina MiSeq Sequencer (Illumina, San Diego, CA, USA) using the MiSeq Reagent Kit v3 with 2 × 300 bp read length. Extracted DNA was mechanically sheared to an aimed average size of 550 bp with the M220 ultrasonicator (Covaris, Wolburn, MA, USA). DNA libraries were prepared using the TruSeq DNA PCR-Free Library Preparation Kit (Illumina) according to the manufacturer’s manual. Quantification of the libraries was done by qPCR using the KAPA Library Quantification Kit for Illumina platforms (Kapa Biosystems, Wilmington, MA, USA) on a Rotor-Gene Q PCR cycler (Qiagen). The genome of CEN.PK113–7D, the in silico constructed watermarked and reference (native) SinLoG sequences and the KlURA3 repair fragment were used as a reference to map sequence reads of genomic DNA onto using the Burrows-Wheeler Alignment tool (BWA).[55] The sequence alignment was further processed using SAMtools.[56] Coverage of the sequence reads was also calculated using the Magnolya algorithm.[57]

All Illumina sequences are available at NCBI (http://www.ncbi.nlm.nih.gov/) under the bioproject accession number PRJNA554743.

Batch Cultivations of IMX1770, IMX1771, and IMX2028

Batch cultivations were performed in biologically independent triplicates in 2-Liter fermenters (Applikon, Delft, The Netherlands) with a working volume of 1.4 L. Cells from exponentially growing SMG shake flask cultures were inoculated into the fermenters containing SMG supplied with 0.2 g·L–1 antifoam Emulsion C (Sigma-Aldrich, St Louis, MO) at an OD660 of 0.4. The fermenters were sparged with dried compressed air at a rate of 700 mL/min (Linde, Gas Benelux, The Netherlands). The broth was stirred constantly at 800 rpm, kept at a constant temperature of 30 °C and at a pH of 5 by automatic addition of 2 M KOH by an Applikon ADI 1030 Bio Controller.

Optical density was measured every hour at 660 nm with a Jenway 7200 spectrophotometer (Staffordshire, United Kingdom). For extracellular metabolite analysis 1 mL of the broth was centrifuged for at least 10 min at 13 000 rpm and the supernatant was analyzed with high-performance liquid chromatography (HPLC) using an Agilent 1100 (Agilent Technologies, Santa Clara, CA, USA) with an Aminex HPX-87H ion-exchange column (BioRad, Veenendaal, The Netherlands) operated with 5 mM sulfuric acid as mobile phase at a flow rate of 0.6 mL/min. The carbon dioxide and oxygen concentration in the gas outflow were analyzed by a Rosemount NGA 2000 analyzer (Baar, Switzerland), after cooling of the gas by a condenser (2 °C) and drying using a PermaPure Dryer (model MD 110-8P-4; Inacom Instruments, Veenendaal, The Netherlands). Biomass dry weight was measured 5–6 times by filtering (pore size 0.45 μm, Gelman Laboratory, Ann Arbor, MI, USA), as described previously.[58] Sampling for RNA was done directly from the reactor in liquid nitrogen as described by Piper et al.(59) The cells were stored at −80 °C for maximally 2 weeks until further processing and RNA was extracted as previously described.[23] An equivalent of 48 mg dry weight per sample was used.

At the same time points as the samples that were taken for RNA isolation, approximately 62.5 mg dry weight was sampled for the enzyme assays, stored at −20 °C in 4 mL aliquots and further process as previously described.[58] Optical densities of the cultures at the moment of sampling for RNA analysis and enzyme assays can be found in Table S9.

Determination of In Vitro Enzyme Activities

On the day of the enzyme assays, frozen samples were thawed and prepared for assays as described by Postma et al.(58) Assays were performed using a U-3010 spectrophotometer (Hitachi, Tokio, Japan) at 30 °C and 340 nm as described by Jansen et al. (2005),[60] with the exception of Pfk, which was performed according to Cruz et al. (2012).[61] Reported activities are based on at least two technical replicates, measured with different cell extract concentrations. When necessary, cell extracts were diluted in 100 mM monopotassium phosphate buffer and with 2 mM magnesium chloride (pH 7.5), or in demineralized water when triose phosphate isomerase activity was measured. The protein concentration of the cell extracts was determined as described by Lowry et al.(62) using bovine serum albumin as a standard.

RNA Sequencing Simulation

To evaluate the watermarking methods, and to compare the markerQuant tool with traditional alignment, we generated artificial RNA-Seq reads from the native and watermarked sequences. Using the polyester R package,[63] we simulated two conditions, in which the second condition has a 4-fold expression of each transcript compared to the first condition. We used the standard error_rate parameter of 0.005. Generated reads were paired-end, each end 100 bp in length.

RNA Sequencing and Data Analysis

RNA libraries and sequencing were performed by Novogene Bioinformatics Technology Co., Ltd. (Yueng Long, Hong Kong). Sequencing was performed using HiSeq 150 bp paired-end reads system using 250–300 bp insert strand-specific library. As described by Novogene, library preparation involved mRNA enrichment using oligo (dT) beads, followed by random fragmentation of the mRNA. cDNA was synthesized from mRNA using random hexamer primers and a second strand synthesis was done applying a custom second strand synthesis buffer (Illumina), dNTPs, RNase H, and DNA polymerase I. After adaptor ligation, double-stranded cDNA library was finalized by size selection and PCR enrichment and samples were sequenced.

Obtained data had an average of 23.08 M reads per sample (Table S8). To quantify the abundance of glycolytic genes with and without watermarks, a similar scheme as the k-mer algorithm of Gehrmann and co-workers[38] was applied. Briefly, for each transcript, we identify sequence markers of 21bp that are unique in the transcript relative to the entire transcriptome and genome. With an exact matching algorithm, these markers can uniquely identify the transcript of origin of a read in RNA-Seq data. In contrast to the previous work,[38] we did not remove overlapping markers (we did not remove redundant markers) but merged them into larger sequences in which any 21 bp k-mer would uniquely identify the transcript of origin. This allowed us to recover a higher percentage of reads per transcript. Gaps in these merged sequences that are not unique relative to the genome and transcriptome were ignored in the marker quantification step. As in previous work, we used an Aho-Corasick exact string-matching algorithm to quantify transcripts. Differential expression was performed using DE-Seq2.[64] All RNA-seq sequences are available at NCBI GEO (https://www.ncbi.nlm.nih.gov/geo/) under GEO accession number: GSE135470.

RNA Sequencing Data Analysis Implementation and Code Availability

The marker discovery and quantification tools were developed in scala, and the entire pipeline is implemented in python using Snakemake.[65] In addition to the k-mer method, a traditional alignment pipeline is also implemented in the markerQuant utility. All code, including an example data set, is available at https://github.com/thiesgehrmann/markerQuant.

Label Free Quantification (LFQ) by Shot-Gun Proteomics

Cultivation and Sampling

For proteomics analysis, the yeast strains pregrown to exponential phase in SMG in shake flask were used to inoculate fresh SMG flasks. Five ml of these cultures in midexponential phase were centrifuged for 10 min at 5000g at 4 °C and the cell pellet was directly stored at −80 °C. Cultures were performed in biological triplicates for strains with watermarked and native glycolysis. To verify if the difference in Gpm1 activity observed in bioreactor between the strains with watermarked and native glycolysis was also present in shake flask culture, Gpm1 activity was assayed in cell samples from the shake flasks. This additional analysis confirmed the lower specific activity of Gpm1 in the watermarked strain (Figure S9).

Protein Extraction and Trypsin Proteolytic Digestion

Cell pellets were resuspended in lysis buffer composed of 100 mM TEAB containing 1% SDS and phosphatase/protease inhibitors. Yeast cells were lysed by glass bead milling and thus shaken 10 times for 1 min with a bead beater alternated with 1 min rest on ice. Proteins were reduced by addition of 5 mM DTT and incubated for 1 h at 37 °C. Subsequently, the proteins were alkylated for 60 min at room temperature in the dark by addition of 50 mM acrylamide. Protein precipitation was performed by addition of four volumes of ice-cold acetone (−20 °C) and proceeded for 1 h at −20 °C. The proteins were solubilized using 100 mM ammonium bicarbonate. Proteolytic digestion was performed by Trypsin (Promega, Madison, WI), 1:100 enzyme to protein ratio, and incubated at 37 °C overnight. Solid phase extraction was performed with an Oasis HLB 96-well μElution plate (Waters, Milford, USA) to desalt the mixture. Eluates were dried using a SpeedVac vacuum concentrator at 45 °C. Dried peptides were resuspended in 3% ACN/0.01% TFA prior to MS-analysis to give an approximate concentration of 250 ng per μL.

Large-Scale Shotgun Proteomics

An aliquot corresponding to approximately 250 ng protein digest was analyzed using an one-dimensional shot-gun proteomics approach.[66] Briefly, the samples were analyzed using a nanoliquid-chromatography system consisting of an EASY nano LC 1200, equipped with an Acclaim PepMap RSLC RP C18 separation column (50 μm × 150 mm, 2 μm), and an QE plus Orbitrap mass spectrometer (Thermo). The flow rate was maintained at 350 nL/min over a linear gradient from 6% to 26% solvent B over 45 min, followed by back equilibration to starting conditions. Data were acquired from 5 to 60 min. Solvent A was H2O containing 0.1% formic acid, and solvent B consisted of 80% acetonitrile in H2O and 0.1% formic acid. The Orbitrap was operated in data depended acquisition mode acquiring peptide signals from 385 to 1250 m/z at 70 K resolution. The top 10 signals were isolated at a window of 2.0 m/z and fragmented using a NCE of 28. Fragments were acquired at 17 K resolution.

Database Search, Label Free Quantification, and Visualization

Data were analyzed against the proteome database from Saccharomyces cerevisiae (Uniprot, strain ATCC 204508/S288C, Tax ID: 559292, July 2018) using PEAKS Studio X (Bioinformatics Solutions Inc.)[2] allowing for 20 ppm parent ion and 0.02 m/z fragment ion mass error, 2 missed cleavages, acrylamide as fixed and methionine oxidation and N/Q deamidation as variable modifications. Peptide spectrum matches were filtered against 1% false discovery rates (FDR) and identifications with ≥2 unique peptides. Changes in protein abundances between both strains IMX1770 and IMX1771 using the label free quantification (LFQ) option provided by the PEAKS Q software tool (Bioinformatics Solutions Inc.).[67] Protein areas were normalized to the total ion count (TIC) of the respective analysis run before performing pairwise comparison between the above-mentioned strains. LFQ was performed using protein identifications containing at least 2 unique peptides, which peptide identifications were filtered against 1% FDR. The significance method for evaluating the observed abundance changes was set to ANOVA. The abundances of the glycolytic enzymes were further visualized as bar graphs using Matlab2018b. The area of the biological triplicates were averaged and standard deviations were represented as error bars.

Data Availability

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (www.proteomexchange.org) via the PRIDE partner repository with the dataset identifier PXD016914.