Robuste Bisulfit-freie Einzelmolekül-Echtzeitsequenzierung von Methyldesoxycytidin auf der Grundlage eines neuartigen hpTet3-Enzyms

Das Vorhandensein des fünften Nukleosids 5-Methyldesoxycytidin (mdC), entweder in Promotorregionen oder im Genkörper beeinflusst den Transkriptionsstatus des entsprechenden Gens.¹ In der Regel führt das Vorkommen von mdC in Promotorregionen zur Inaktivierung des betreffenden Gens, während unmethylierte Promotoren eine aktivere Transkription anzeigen. Die Identifizierung von mdC in Genen ermöglicht die Charakterisierung des Transkriptionszustands des betreffenden Gens, was für die Identifizierung und Charakterisierung von Tumorzellen von entscheidender Bedeutung ist.^{2, 3} In Tumorzellen werden Onkogene oft fälschlicherweise aktiviert, während Tumorsuppressorgene irrtümlich ausgeschaltet werden. Daher ist die Sequenzierung von mdC mit minimalem Ausgangsmaterial von großem Interesse, um einen neuen Bereich der Tumordiagnostik, die sogenannte Flüssigbiopsie, zu etablieren.⁴

Bislang wird die mdC-Sequenzierung überwiegend mit einer Bisulfit-Behandlung durchgeführt (Abbildung 1a). Genomische DNA, die bei >60 °C mit Bisulfit behandelt wird, wandelt alle unmethylierten Cytidine in Uracil um, während mdC intakt bleibt. Die Bestimmung der mdC-Positionen im Genom erfolgt nach PCR und Sequenzierung durch einen Vergleich der erhaltenen Sequenzen mit einem Referenzgenom. Ein wesentliches Problem dieser Methode besteht jedoch darin, dass ein großer Teil der genomischen Ausgangs-DNA die harten Bedingungen der Bisulfit-Behandlung aufgrund starker Fragmentierung nicht übersteht. Diese Einschränkung wird durch eine umfangreiche PCR-basierte Amplifikation der nicht abgebauten DNA gemildert. Ein weiterer Nachteil ist, dass das Bisulfit-Sequenzierungsprotokoll umständlich und fehleranfällig ist. Mildere Methoden wie EM-seq, die derzeit entwickelt werden, verwenden das desaminierende Enzym APOBEC3 A (A3 A), das ebenfalls dC zu dU desaminiert.⁵ Durch die Desaminierung aller dC-Basen zu dU wird die Komplexität des Genoms jedoch von einem vierstelligen Code auf einen Code aus drei Nukleobasen (dA, dG und dU sowie das verbleibende mdC) reduziert, was die Sequenzzuordnung insbesondere bei repetitiven Elementen erschwert.

Ein alternativer Ansatz zur Sequenzierung von mdC ist die Sequenzierung der dritten Generation, bei der die Sequenzen direkt ausgelesen werden, ohne dass ein PCR-Schritt erforderlich ist. Derzeit ermöglichen alle Einzelmolekül-Sequenzierungswerkzeuge der dritten Generation, wie Nanopore oder SMRT-Sequenzierung, das direkte Auslesen von mdC.⁶

Diese Methoden befinden sich jedoch noch im Anfangsstadium, und die Unterschiede zwischen den für dC und mdC erhaltenen Sequenziersignalen sind oft minimal. Dies erfordert eine umständliche Dekonvolution der Daten, für die ein erheblicher Bedarf an Bioinformatik besteht.^{7, 8} Daher ist die Entwicklung von milden mdC-Sequenzierungsmethoden, welche die Cytidin-Desaminierung umgehen, für die Entwicklung neuer Instrumente für die frühe Tumordiagnose äußerst wünschenswert.

Wir stellten die Hypothese auf, dass die derzeitigen Einschränkungen bei der Sequenzierung von mdC durch die quantitative Oxidation von mdC zu 5-Carboxycytidin (cadC, Abbildung 1b) überwunden werden könnten. Im Gegensatz zu dC verfügt cadC über eine zusätzliche Carboxylgruppe, die unter neutralen pH-Bedingungen negativ geladen ist. Dies bewirkt vermutlich einen signifikanten Signalunterschied zwischen dem neutralen mdC und dem negativ geladenen cadC. Das Konzept, cadC anstelle von mdC durch SMRT-Sequenzierung nachzuweisen (Abbildung 1c), wurde bereits vorgeschlagen,⁹ bislang jedoch nicht umgesetzt.

Hier zeigen wir, dass cadC die beiden kinetischen Parameter, die bei der SMRT-Sequenzierung ermittelt werden, nämlich die Interpulsdauer (IPD) und die Pulsbreite (PW), signifikant erhöht und es uns ermöglicht, einen Deep-Learning-Algorithmus zu trainieren. Unsere Ergebnisse zeigen, dass die Genauigkeit der neuen Methode alle aktuellen mdC-Sequenzierungsmethoden übertrifft und den Weg für eine sanfte epigenetische mdC-Sequenzierung ebnet.

Grundlage dieser Methode ist ein neu entwickeltes verkürztes, aber robustes Tet3-Enzym, das in E. coli überexprimiert werden kann und mdC im Genom mit über 99 % Ausbeute zu cadC oxidiert.

Die Implementierung der neuen Tet-basierten Technologie erfolgte unter Zuhilfenahme der SMRT-Sequenzierung, bei der eine Polymerase die Nukleotide in der Vorlage mit fluoreszenzmarkierten eingehenden Triphosphaten paart. Ein Detektor misst das Fluoreszenzsignal des im aktiven Zentrum der Polymerase gebundenen Triphosphats in Echtzeit, bevor die Fluoreszenzmarkierung während der Bildung der Phosphodiesterbindung abgespalten wird. Da das zu sequenzierende DNA-Fragment in eine zirkuläre Struktur eingebettet ist (Abbildung 1c), bewegt sich die Polymerase bei der Sequenzierung mehrfach entlang der zirkulären Vorlage, so dass jede Base (einschließlich der cadC-Base) wiederholt gelesen wird, was eine große Anzahl von Datenpunkten für jede Base und ein hohes Maß an Genauigkeit liefert. Neben dem Fluoreszenzsignal, das die ankommende Base identifiziert, zeichnet die SMRT-Sequenzierung auch die Zeit auf, die die Polymerase benötigt, um die Phosphodiesterbindung zu bilden (PW-Wert), sowie die Zeit zwischen den einzelnen Inkorporationsereignissen (IPD-Wert), wodurch mehrere Parameter für jede zu sequenzierende Base, einschließlich cadC, welches durch Tet-induzierte Oxidation aus mdC gebildet wird, zur Verfügung stehen. Das erste Problem, mit dem wir uns befassten, war die Entwicklung eines robusten Enzyms, das in der Lage ist, mdC quantitativ zu cadC zu oxidieren. Dies wurde mit Hilfe eines Ten-Eleven-Translocation-Enzyms erreicht, einer Fe²⁺ und α-Ketoglutarat-abhängigen Monooxygenase.¹⁰ Die Oxidation von mdC zu cadC wurde bislang unter Zuhilfenahme der Enzyme Tet1 und Tet2 durchgeführt, deren Überexpression jedoch mit Schwierigkeiten verbunden ist.^{5, 11-14}

Im Gegensatz dazu ist das im Gehirn von Erwachsenen am weitesten verbreitete Tet-Enzym Tet3,¹⁵ welches erhebliche Mengen von mdC zu 5-Hydroxymethyl-dC (hmdC) oxidiert.¹⁶ Wir haben ein strukturgeleitetes Design auf der Grundlage der Kristallstruktur von menschlichem Tet2¹⁷ durchgeführt (siehe SI) und identifizierten eine stark (um ~72 %) verkürzte Maus-Tet3-Variante (hpTet3), die aus nur 465 Aminosäuren (52 kDa) besteht. Bei dieser hpTet3-Variante ersetzten wir die Region mit geringer Komplexität innerhalb der katalytischen Domäne (cd) durch eine Glycin-Serin-Brücke (Abbildung 2). Das resultierende hpTet3-Protein konnte in E. coli überexprimiert werden. Ähnliche Ansätze von zwei anderen Gruppen^{18, 19} weichen in Bezug auf die vorgenommenen Verkürzungen deutlich von unserem Ansatz ab (siehe Hintergrundinformationen für weitere Einzelheiten). Das mit einem N-terminalen Strep-Tag fusionierte hpTet3 wurde zunächst durch Affinitätschromatographie über ein StrepTrap XT-Material gereinigt und anschließend eine Heparin-Säule verwendet, um Chaperon-Verunreinigungen (Hsp40 und Hsp70) zu entfernen.

Als nächstes untersuchten wir die katalytische Aktivität des hpTet3-Enzyms. Zu diesem Zweck methylierten wir zunächst mit der CpG-spezifischen Methyltransferase M.SssI die genomische DNA des Lambda-Phagen und verdauten sie bis auf die Ebene der einzelnen Nukleoside. Die Quantifizierung aller im Genom vorhandenen Nukleoside, insbesondere potenziell oxidierter Nukleoside wie 5-Formyl-dC (fdC), 5-Hydroxymethyl-dU (hmdU) und 8-Oxo-dG, erfolgte mittels quantitativer Triple-Quadrupol-Massenspektrometrie (UHPLC-QQQ-MS) unter Verwendung eines vollständigen Satzes mit stabilen Isotopen markierter interner Standards für dA, dT, dG, dC, mdC, hmdC, fdC, cadC, hmdU und 8-oxo-dG (Abbildung 3) nach unserer zuvor beschriebenen Methode.²⁰ In einem zweiten Experiment behandelten wir die genomische DNA vor dem Verdau mit hpTet3 (zu den Behandlungsbedingungen siehe SI) und wiederholten das Quantifizierungsexperiment mit den genannten Isotopenstandards. Dadurch konnten wir sehr genaue quantitative Daten erhalten. Wie in Abbildung 3a dargestellt, haben wir festgestellt, dass bei der Oxidation das Signal für mdC verschwindet und stattdessen ein starkes Signal für cadC zu erkennen ist. Die Oxidationszwischenstufen hmdC und fdC wurden nicht gefunden, was eine vollständige Umwandlung beweist. Wir wiederholten die Studie mit verschiedenen eukaryotischen Genomen (Abbildung 3b) und stellten fest, dass in allen Fällen das Signal für mdC (und 5-Hydroxymethyl-dC) bei der Oxidation mit hpTet3 vollständig verschwand und ein neues und starkes Signal für cadC entstand.

In allen Fällen konnten wir kein restliches mdC nachweisen. Stattdessen wurde cadC in Konzentrationen nachgewiesen, die denjenigen von mdC im Ausgangsmaterial entsprechen, z. B. 4,23 % in HEK293T gDNA (siehe Hintergrundinformationen Table S1 für andere Proben). Die quantitative Analyse der hpTet3-Reaktion ergab eine Oxidationsausbeute von mdC zu cadC von 99,96 % in genomischer DNA.

Um unerwünschte oxidative Schäden abzuschätzen, haben wir den Gehalt an 8-oxo-dG quantifiziert und keinen signifikanten Anstieg festgestellt. Bei hmdU wurde ein erwartungsgemäßer Anstieg auf 48 hmdU-DNA-Läsionen pro Genom beobachtet.²¹ Diese Daten zeigen, dass die Oxidation von mdC zu cadC mit hpTet3 sehr effizient ist und nur eine geringe Anzahl von hmdU-Läsionen erzeugt (SI Abb. S3 & Table S2).

In der Folge wurde die Möglichkeit einer Sequenzierung von mdC mittels SMRT-Sequenzierung der dritten Generation evaluiert. Für die SMRT-Sequenzierung haben wir drei Modellgenome aus Lambda-Phagen-DNA (dam^-, dcm^-) hergestellt. Das erste Genom (LMD-dC) enthielt kein mdC. Im zweiten Genom haben wir alle CpG-Stellen mit M.SssI enzymatisch methyliert (LMD-mdC). Für das dritte Genom oxidierten wir das LMD-mdC-Genom mit hpTet3, um alle mdCs in cadCs umzuwandeln (LMD-cadC). Um die Reinheit dieser Modell-DNAs zu überwachen, verdauten wir die drei Genome bis auf Nukleosid-Ebene und analysierten die Nukleosid-Zusammensetzung mit UHPLC-QQQ-MS. Die erhaltenen Daten (SI-Tabelle S1) bestätigten die hohe Methylierungseffizienz von M.SssI (99,99 %) und die hohe Oxidationseffizienz zu cadC mit hpTet3 (99,96 %).

Nach der Bibliotheksvorbereitung und Sequenzierung auf einer Sequel IIe-Plattform erhielten wir 385.942 Sequenzen für das unmodifizierte Lambda-Genom, 433.815 Sequenzen für das mdC-haltige Genom und 456.646 Sequenzen für das cadC-haltige Genom. Als Nächstes führten wir die Extraktion von Alignment-Merkmalen (IPD- und PW-Werte) mit ccsmeth (https://github.com/PengNi/ccsmeth) durch, wie in Abb. 4a dargestellt. Wir extrahierten mittlere IPD- und PW-Werte für jede CpG-Position im Lambda-Genom innerhalb eines 21-k-mers (±10 bp um das interessierende Cytosin) (Abb. 4b). Hinsichtlich der IPD-Werte wurde ein großer Unterschied im 21-k-mer zwischen den Situationen dC, mdC und cadC beobachtet (Abb. 4b links). Im Vergleich zu dC zeigten mdC und cadC unterschiedliche Signalmuster in der Nähe der xdC-Position (k-mer-Positionen 8–19). Während die Signalmuster für dC und cadC ähnlich sind, sind die normalisierten Zeitwerte für cadC stark erhöht. Außerdem wurden Unterschiede in den PW-Mustern beobachtet (Abb. 4b rechts). Insbesondere zeigt cadC einen starken Zeitanstieg an der k-mer-Position 18, die 7 Positionen stromabwärts von der cadC-Position liegt. Diese Daten zeigen, wie komplex die Unterschiede im Fußabdruck zwischen dC, mdC und cadC sind, insbesondere außerhalb von CpG-Dyaden. Denn die PW- und IPD-Unterschiede zeigen sich nicht nur am Nukleotid selbst, sondern auch mehrere Nukleotide entfernt, stromaufwärts oder stromabwärts (SI Abb. S4). Um sicherzustellen, dass diese Muster nicht das Ergebnis überlappender Signale benachbarter CpG-Stellen sind, haben wir IPD und PW für 21-k-mere mit einem, zwei oder mehr als zwei CpGs aufgetragen (SI Abb. S5). In allen Fällen sind die extrahierten Muster für IPD und PW ähnlich, was auf vielfältige und komplexe Wechselwirkungen zwischen der Polymerase und dem Templat hinweist.

Als Nächstes untersuchten wir, ob die komplexen, aber starken kinetischen PW- und IPD-Daten für cadC zum Training eines KI-basierten rekurrenten neuronalen Netzwerks (RNN) verwendet werden können. In Anlehnung an die Trainingsprozedur von ccsmeth trainierten wir ein mdC-Modell auf der Grundlage der LMD-dC- und LMD-mdC-Datensätze, sowie ein cadC-Modell auf der Grundlage der LMD-dC- und LMD-cadC-Datensätze (Abb. 4a). Die aus dem KI-Modell gewonnenen Schlüsselparameter sind in Hintergrundinformationen Abb. S6 dargestellt. Zu unserer Freude stellten wir fest, dass die kinetischen Sequenzierungsdaten von cadC in Kombination mit dem trainierten Algorithmus ein cadC-Modell ergaben, das die Leistung des kanonischen ccsmeth und unseres mdC-LMD-Modells in allen Aspekten übertraf. Zur Erkennung von mdC benötigte das RNN 181 Trainingsrunden, um eine Modellgenauigkeit von 0,945 zu erreichen. Die Präzision erreichte einen Wert von 0,962 und der Recall (Anzahl der beschreibbaren CpG-Dyaden) lag bei 0,956. Für cadC benötigte das Modell nur 16 Trainingsschritte, um eine Genauigkeit von 0,987, eine Präzision von 0,987 und einen Recall von 0,988 zu erreichen. Dies sind fantastische Werte, insbesondere angesichts der wenigen erforderlichen Trainingsschritte (SI Abb. S6).

Anschließend testeten wir die neuen Modelle, indem wir sie auf die einzelnen Datensätze für LMD-dC, LMD-mdC und LMD-cadC anwandten. Für die Erkennung von mdC verglichen wir unsere LMD-mdC- und LMD-cadC-Modelle mit dem ccsmeth-Standardmodell, das auf menschlicher DNA trainiert wurde.²² Anhand unserer UHPLC-QQQ-MS-Messungen wissen wir, dass das modifizierte genetische Material, das den LMD-mdC- und LMD-cadC-Modellen zugrunde liegt, 99,99 % mdC bzw. 99,96 % cadC in den CpG-Dyaden enthält (SI Tabelle S1). Im LMD-mdC-Modell sollte die Methylierungshäufigkeit einer CpG-Dyade folglich fast 1 betragen, während sie in der nicht-methylierten LMD-dC-DNA nahe null liegen sollte. Die Sequenzierungsergebnisse sind in den Abbildungen 4c und S7 zusammengefasst und in den SI-Tabellen S3 und S4 im Detail dargestellt.

Als Referenz haben wir LMD-dC und LMD-mdC mit Bisulfit sequenziert, um unsere mit der individuellen ccsmeth-Prozedur erzielten Ergebnisse damit zu vergleichen. Die Bisulfit-Daten zeigen, dass 98,2 % aller CpGs eine Methylierungsfrequenz von >=90 % aufweisen. Allerdings weisen nur 23,82 % aller CpGs eine Methylierungshäufigkeit von >=95 % auf.

Bei der Analyse der SMRT-seq-Daten mit unserem auf mdC trainierten Modell stellten wir fest, dass 76,75 % der CpG-Stellen in der LMD-mdC-Probe eine Methylierungshäufigkeit von >95 % aufwiesen. Darüber hinaus wiesen 21,1 % der CpGs in der LMD-dC-Probe eine Methylierung von mehr als 10 % auf. Die Zahlen waren sogar noch niedriger, wenn wir das kanonische ccsmeth-Modell anwandten, bei dem nur 30 % der CpG-Stellen mit einer Methylierungshäufigkeit von >95 % berichtet wurden (SI Tabelle S3).

Mit dem neuen cadC-Modell, das auf hpTet3-oxidierter Lambda-DNA trainiert wurde, erkennen wir ein scharfes Signal (blau) von nahezu 100 % für die cadC-Häufigkeit in CpG-Dyaden (Abbildung 4c). Genauer gesagt haben 94,68 % aller CpGs aus der LMD-dC-Probe eine Modifikationsfrequenz von <10 % und 99,95 % aller CpGs aus der LMD-cadC-Probe eine Modifikationfrequenz von >=95. Tatsächlich sagt unser Modell voraus, dass 94,63 % aller CpGs eine Methylierungsfrequenz von >=98 % haben.

Diese Ergebnisse zeigen, dass das neue Modell in der Lage ist, alle cadCs im Lambda-Genom mit bisher unerreichter Genauigkeit vorherzusagen. Dies beruht auf der Tatsache, dass die Oxidation von mdC zu cadC mit dem neuen hpTet3-Enzym sehr charakteristische PW- und IPD-Werte liefert. In Kombination mit dem KI-generierten cadC-Modell ermöglicht dies eine sehr genaue und empfindliche Sequenzierung von mdC. Im Gegensatz zu anderen bisulfitfreien Methoden^{5, 11, 23} ist unser Ansatz einstufig und PCR-frei, was mögliche Verzerrungen reduziert.

Wir haben festgestellt, dass eine sequenzabhängig abweichende Oxidationsaktivität des hpTet3-Enzyms (SI Abb. S8 & S9 und Hintergrundinformationen Table S5 & S6) gering ist, insbesondere wenn hohe Enzymkonzentrationen verwendet werden. Dies ist ein weiterer Grund für die hohe Genauigkeit (99,9 %) unserer 5mdC-Sequenzierungsmethode.

Zusammenfassend beschreiben wir die Entwicklung einer neuen hpTet3-Variante, die in E. coli überexprimiert werden kann. hpTet3 ist stabil und oxidiert mdC zu cadC mit hoher Effizienz. Der Sequenzeinfluss auf die Oxidation ist gering. Vor allem bei höheren Enzymkonzentrationen sind wir in der Lage, mdCs mit hoher Effizienz zu cadC zu oxidieren, auch in Nicht-CpG-Kontexten. Somit kann das Enzym in Kombination mit einer Desaminierungsreaktion und PCR prinzipiell auch die mdC-Sequenzierung in Nicht-CpG-Kontexten unterstützen. Wir haben die SMRT-Sequenzierung von hpTet3-behandelter DNA durchgeführt und sehr charakteristische IPD- und PW-Werte festgestellt, die eine präzise Lokalisierung von cadC in allen CpG-Kontexten ermöglichten. Die SMRT-Sequenzierung von mdC über cadC wurde nach weiteren Manipulationen mit noch nie dagewesener Genauigkeit erreicht. Eine weitere Stärke der Methode ist die Möglichkeit der Sequenzierung mit besonders hoher Leselänge, die auch für repetitive Genomelemente hochpräzise Daten liefert. Das neue hpTet3 kann außerdem dazu beitragen, andere 5mdC-Sequenzierungsstrategien, wie TAPS und EM-Seq, zu verbessern. In Analogie dazu nennen wir unsere Methode Enzymatische Methyl-Oxidations-Sequenzierung (EMox-seq).

Die in dieser Veröffentlichung besprochenen Daten sind in der Supporting Information verfügbar. Sequenzierdaten wurden im Gene Expression Omnibus des NCBI hinterlegt²⁴ hinterlegt und sind über die GEO-Serien-Zugangsnummer GSE256446 zugänglich. (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE256446).