Key aurophilic motif for robust quantum-tunneling-based characterization of a nucleoside analogue marker†

A quantum sequencer offers a scalable electrical platform for single-molecule analysis of genomic events. A thymidine (dT) analog exhibiting uniquely high single-molecule conductance is a key element in capturing DNA synthesis dynamics by serving as a decodable marker for enzymatic labeling of nascent strands. However, the current design strategies of dT analogs that focus on their molecular orbital energy levels require bulky chemical modifications to extend the π-conjugation, which hinders polymerase recognition. We report herein a polymerase-compatible dT analog that is highly identifiable in quantum sequencing. An ethynyl group is introduced as a small gold-binding motif to differentiate the nucleobase–gold electronic coupling, which has been an overlooked factor in modifying nucleobase conductance. The resulting C5-ethynyl-2′-deoxyuridine exhibits characteristic signal profiles that allowed its correct identification at a 93% rate while maintaining polymerase compatibility. This study would expand the applicability of quantum sequencing by demonstrating a robust nucleoside marker with high identifiability.

For robust calculation of CV values, the signals with more than 10 data points (1 ms) were adopted for further statistical analysis. In addition, signals that are out of the other thresholds (gray box) were excluded from the analysis because of their abnormalities. (b) Schematic flow of base-assignment processes. The signal data were divided into five segments, four of which were adopted as a training data set and the remainder of which was used as a test data set. The probability density functions (PDFs) of the signal features for each of the nucleosides were estimated by kernel density estimation. Based on the PDFs, the probabilities to be assigned to each of the nucleobase was calculated for each signal in the test data set. When the signal coverage = n (n = 1-50), the final probabilities were estimated by naïve Bayes estimation from each of the probabilities calculated for the n signals randomly picked from the test data set. Then, the nucleoside with the largest probability was assigned as the most probable assignment. Upper and lower limits of the box plots are corresponding to the 90th and 10th percentiles, respectively. Fig. S5b and S5c showed that the oligonucleotide containing EtdU provides signals with high Ave and/or CV values more frequently than the oligonucleotide lacking EtdU. This indicates that EtdU in the oligonucleotide also exhibits statistically characteristic signals regarding Ave and CV compared with canonical dT and dG. Plots of correct assignment rates for each nucleoside against the signal coverage in 1D-statistical processing of Ave. (c) Plots of correct assignment rates for each nucleoside against the signal coverage in 2D-statistical processing of Ave and CV. Correct assignment rates were converged rapidly within fifty-signal coverage by the 2D signal processing compared with 1D signal processing. These observations indicated that the current fluctuation represented by CV is also an important signal feature in addition to the signal intensity for the accurate identification of the dT analog.

Supplementary tables
2D-histograms & box plots ( Fig. 4c and S5) The number of current plateaus

Preparation of gap electrodes
The sensor plate with a gold-nano wire was purchased from SCREEN Holdings Co. Ltd.. The gold gap electrodes were prepared by mechanically controlled break junction (MCBJ) method as described in the literature. 3 Briefly, the nano-gap devices were fabricated as follows. First, spincoating polyimide film was applied as an insulating layer on thin-silicon substrate. A nano-gold wire with a width of several tens of nanometers was formed by electron-beam lithography. SiO2 film was deposited on the gold-wire by chemical vapor deposition. After patterning, dry-etching was performed to form free-standing gold-wire. According to the SEM observation, the etched depth is about 1 μm. The free-standing gold bridge in SiO2 plate was subjected to repetitive threepoint bending. After several ten times of bending, the gold junction was mechanically broken to generate a pair of gap electrode, which was monitored by the sudden drop of current flow in the electric circuit. The size of gap electrodes was precisely controlled in real time during the current measurement by fine-tuning of lifting bar using a piezoelectric element.

Current measurements of nucleosides, oligonucleotides and p-disubstituted benzenes
Tunnel current measurements were conducted using 1 µM deionized aqueous solutions of the molecules, where concentration indicates the strand concentration for the oligonucleotide. In preparation of nucleoside and p-disubstituted benzene samples, the compounds were firstly dissolved in DMSO (1 mM) and then the solution was diluted with deionized water by 1000-fold S12 to adjust final concentration to 1 µM. For electrical measurements, a 20-µL sample solution was dropped onto the center of the sensor plate where the gap electrodes were generated. No electrophoretic force was applied to the sample solution and molecules depended on Brownian motion to access gap electrodes. Gap size was set to 0.6 nm and adjusted by the piezoelectric element during all measurements. Bias voltage between a pair of electrodes was 0.1 V. The time of each set of the measurement was typically 1 h. Measurements were carried out more than three sets for each molecule using different gold gap sensors (Table S1).

Picking and analysis of tunnel current signals
A single signal pulse (S) was defined as the current region from the point which exceeded the first threshold (η ! ( )) to the point which fell below the second threshold (η " ( )). To determine those thresholds, the adaptive threshold analysis was adopted from the literature. 4 In this process, the baseline current (µ( )) and the standard deviation of baseline current (σ( )) that were used to define the η1 and η2 values were recursively determined according to the following equations; where ( ) is the current value at the -th position of the data point and is a regulatory parameter that controls the influence of ( ). In this experiment, 1/30 was adopted as the value. Using the µ( ) and σ( ) values, the thresholds at the -th position of the data point (η ! ( ) and η " ( )) were defined as follows.
η ! ( ) = µ( ) + 6 × σ( ) Those values were not updated in the signal region (S) and restarted to be updated just after the data points went below η " ( ). The averaged intensity from the baseline current (Ave) and the coefficient of variation of the signal (CV) were calculated from the data points included in each S.
Those parameters were adopted as the signal features to characterize the nucleosides in this study.
The signals with more than 10 data points (corresponding to 1 ms) were selected for statistical analysis of Ave and CV because short signals are not appropriate for calculation of CV. The signals whose features were out of the thresholds shown in Fig. S3a were also excluded from the further S13 statistical analysis and the base-assignment procedures due to their abnormalities or difficulties to be discriminated from noise signals.

Change point analysis and the detection of current plateaus for oligonucleotides
The multistep signals of the oligonucleotide may correspond to more than one nucleotide in the sequence. To extract the current plateaus that might correspond to each of the single nucleotides trapped between the gap, package "changepoint" 1 of R language was adopted (Fig. S5a). Ave and CV values were independently calculated for each of the plateaus and assembled to construct 2D histograms (Fig. 4c). Specifically, the "cpt.meanvar" function was used with Pruned Exact Liner Time method (PELT) 5 to detect the change points of the intensity and the variance in the current traces. In this process, the change points were assigned to the points with less than 2% of the type Ⅰ error rate by specifying "Asymptotic" as a penalty argument with pen.value = 0.02. Then, the current plateaus with more than 10 data points were adopted for the further statistical analysis of Ave and CV.
Base-assignment to the nucleoside signals ( Fig. 3 and S7) Base-assignment to the nucleoside signals were carried out based on the two-dimensional characterization of the signal features (Ave and CV) and one-dimensional characterization of the signal feature (Ave) in Fig. 3c and S7. Throughout this process, the Ave (pA) and CV values were converted to the base 10 logarithms. The kernel density estimation of probability density functions of the signal features was conducted based on "gaussian_kde" of the scipy.stats module in the "SciPy 1.0" package. 6 The general descriptions of the kernel density estimation are provided at "Univariate and bivariate kernel density estimation" of this section. The detailed procedures for discriminating the m nucleosides (C1,···, Ck,···, Cm) are described below.
First, the procedures to omit noise signals was conducted. For this purpose, the signals were categorized into two groups, the signals from background measurements and nucleoside samples.
2D or 1D probability density functions of the signal features (PDFs) were independently constructed by kernel density estimation for background and nucleoside groups. Then, the probabilities to be classified as a background noise and a nucleoside signal were calculated for each signal from the nucleoside group based on the PDFs. The signals that were classified as the background noise with no less than 70% probabilities were excluded from the subsequent baseassignment processes.
The signal data sets of each of the nucleosides (C1,···, Ck,···, Cm) were independently divided into five segments, four of which were grouped as a training data set, and the rest of which was a test data set. From the training data, the 2D or 1D probability density functions of signal features S14 were independently estimated for each of the nucleosides using kernel density estimation.
Subsequently, the test data signals of each of the nucleosides were randomly picked and grouped by every n signals ( ! ,···, # ) as one signal set where n represents the signal coverage. Then, the probability for each set of the signals to be classified as the k-th nucleoside, ( $ | ! ,···, # ), was calculated according to the naive Bayes estimation as described below; where is the normalization constant, ( $ ) is the prior probability (0.5 for binary, 0.33 for ternary and 0.2 for quinary), and ( $ | % ) is the probability to be classified as $ for the i-th signal ( % ) of the signal set. Finally, the nucleoside with the largest ( $ | ! ,···, # ) was assigned to the signal set. The random grouping of the test signals into signal sets and base-assignment to all the signal sets were repeated n times when the signal coverage = n.
A series of the base-assignment processes described above were repeated five times until all possible segmentations of the five data segments into training and test data sets were tested. Finally, for each of the k nucleosides, the percentage of correct and incorrect assignments was calculated for each nucleoside from the total number of correct and incorrect base-assignments and used collectively to provide a confusion matrix.

Univariate and bivariate kernel density estimation [Univariate kernel density estimation]
Estimation of 1D density function ( ( )) of Ave is performed by the kernel density estimation described by the following equation.
Here, is a variate (Ave in this study), is the number of signals, ( ) is a kernel function, ℎ is a bandwidth parameter of the kernel function and ' ( ) = ( /ℎ)/ℎ . In this study, ( ) corresponds to the gaussian function represented as follows. The kernel density estimation can be extended to the multivariate cases. The general description of the probability density function obtained from multivariate kernel density estimation ( S ( )) are the following equations.
where ℎ $ is a bandwidth parameter of the k-th variate, ( ) is a kernel function, and d is the number of the variables (d = 2, in the case of 2D processing of Ave and CV).

Primer extension
The primer (0.4 µM) and template (0.4 µM) were mixed in a 10 µL of reaction mixture containing 2×reaction buffer for Ex Taq polymerase (Takara Bio) and 1 unit of Ex Taq polymerase (Takara Bio). The primer and template were denatured by heating at 95 ºC for 2 minutes and then annealed by cooling to 54 ºC at the rate of -0.1 ºC/sec. A primer extension reaction was started by adding a 10 µL of 400 µM dNTP mix containing dTTP (Sigma-Aldrich) or EtdUTP (abcam) as a genetic alphabet "T" and the reaction mixtures were kept at 54 ºC. After a defined time, the reaction was stopped by adding a 20 µL of loading buffer (7 M urea and 20 mM of EDTA in TBE buffer). A 10 µL of the resulting mixture was analyzed by 15% denaturing polyacrylamide gel electrophoresis. The FAM-modified primers and extended products were detected by a gel imager.

PCR amplification of human TNFα gene
Human TNFα gene (702 bp) was amplified with PCR using pBluescript SKⅡ hTNFα as a template. PCR was performed in 20 μL total reaction volume with primers (1 µM), dNTP mix containing dTTP or EtdUTP (250 µM), 1×reaction buffer for Ex Taq polymerase (Takara Bio) and 1 unit of Ex Taq polymerase (Takara Bio). The reaction mixture was first heated at 95 ºC for 2 minutes and followed by 30 cycles of denaturing (95 ºC for 30 seconds), annealing (50 ºC for 30 seconds) and elongation (72 ºC for 60 seconds) steps. protecting group was cleaved using a standard DMTr off procedure. Then, the CPG beads were S17 treated with 28% NH3 aqueous solution at room temperature for 1 hour to cleave off the oligonucleotide. The resulting solutions were heated at 55 ºC for 16 hours. After cooling down, the samples were lyophilized. To deprotect the TIPS protecting group from the EtdU precursor, the oligonucleotide was dissolved in DMF (200 µL) and heated at 50 ºC for 30 minutes. Then, 50

Synthesis and purification of oligonucleotides
µL of 1 M tetrabutylammonium fluoride (TBAF) in THF was added to the solution. After the mixture was incubated at 45 ºC for 15 minutes, the reaction was quenched by adding 250 µL of 2 M triethyl ammonium acetate (TEAA) buffer. Then, the solution was supplied to a purification step by the reverse phase HPLC. The product was eluted by 100 mM TEAA (pH 7.0)/ACN (ACN: 2% for 5 min, linear gradient of ACN: 2%-25% for 23 minutes). The elution of the oligonucleotide was monitored by UV absorbance at 260 nm and the fraction containing the target product was lyophilized. The resulting dry oligonucleotide was dissolved in deionized water and characterized by MALDI-TOF MS (Table S2).