A single-cell survey of the small intestinal epithelium

Mäuse

Alle Arbeiten an Mäusen wurden in Übereinstimmung mit den Institutional Animal Care and Use Committees (IACUC) und den einschlägigen Richtlinien des Broad Institute und des Massachusetts Institute of Technology durchgeführt, mit den Protokollen 0055-05-15 bzw. 0612-058-18. Für alle Experimente wurden die Mäuse nach dem Zufallsprinzip den Behandlungsgruppen zugewiesen, nachdem Geschlecht und Alter der 7-10 Wochen alten weiblichen oder männlichen Wildtyp-C57BL/6J- oder Lgr5-EGFP-IRES-CreERT2 (Lgr5-GFP)-Mäuse, die vom Jackson Laboratory (Bar Harbour) bezogen wurden, oder Gfi1beGFP/+ (Gfi1b-GFP)-Mäuse43 ermittelt worden waren. Die Mäuse wurden unter spezifisch pathogenfreien Bedingungen in den Tiereinrichtungen des Broad Institute, Massachusetts Institute of Technology oder der Harvard T. H. Chan School of Public Health untergebracht.

Salmonella enterica und H. polygyrus-Infektion. C57BL/6J-Mäuse (Jackson Laboratory) wurden mit 200 Larven des dritten Stadiums von H. polygyrus oder 108 Salmonella enterica infiziert, die im Massachusetts General Hospital (Charlestown) unter spezifisch-pathogenfreien Bedingungen gehalten wurden, gemäß Protokoll 2003N000158. H. polygyrus wurde wie zuvor beschrieben44 vermehrt. Die Mäuse wurden 3 und 10 Tage nach der Infektion mit H. polygyrus eingeschläfert. Für Salmonella enterica wurden Mäuse mit einem natürlich streptomycinresistenten SL1344-Stamm von S. Typhimurium (108 Zellen) infiziert, wie zuvor beschrieben44, und 48 Stunden nach der Infektion eingeschläfert.

Zelldissoziation und Kryptenisolierung

Kryptenisolierung. Der Dünndarm von C57BL/6J-Wildtyp-, Lgr5-GFP- oder Gfi1b-GFP-Mäusen wurde isoliert und in kaltem PBS abgespült. Das Gewebe wurde der Länge nach geöffnet und in kleine Fragmente von etwa 2 mm Länge geschnitten. Das Gewebe wurde 90 Minuten lang in 20 mM EDTA-PBS auf Eis inkubiert, wobei es alle 30 Minuten geschüttelt wurde. Anschließend wurde das Gewebe kräftig geschüttelt, und der Überstand wurde als Fraktion 1 in einem neuen konischen Röhrchen aufgefangen. Das Gewebe wurde in frischem EDTA-PBS inkubiert, und alle 30 Minuten wurde eine neue Fraktion gesammelt. Die Fraktionen wurden so lange gesammelt, bis der Überstand fast ausschließlich aus Krypten bestand. Die letzte Fraktion (angereichert mit Krypten) wurde zweimal in PBS gewaschen, 3 Minuten lang bei 300 g zentrifugiert und 1 Minute lang bei 37 °C mit TrypLE express (Invitrogen) dissoziiert. Die Einzelzellsuspension wurde dann durch einen 40-μm-Filter gegeben und für FACS für scRNA-seq (unten) gefärbt oder für die Organoidkultur verwendet. Wir bestätigten die Robustheit dieser Methode, indem wir weitere Einzelzellisolierungsmethoden testeten – entweder „ganz“ (Abschaben der Epithelauskleidung) oder „zottenangereichert“ (Fraktion 1; siehe oben) – und stellten fest, dass aufgrund der hohen Sterblichkeitsrate (durch Anoikis) von post-mitotischen differenzierten Zellen (deren Hauptbestandteil reife Enterozyten sind) die zottenangereicherte Einzelzellsuspension die Zusammensetzung der Dünndarmzelltypen getreu wiedergibt (Daten nicht gezeigt).

Isolierung von Follikel-assoziierten Epithelien. Epithelzellen aus den Follikel-assoziierten Epithelien wurden isoliert, indem kleine Schnitte (0,2-0,5 cm) mit Peyer’schen Flecken aus dem Dünndarm von C57Bl/6J oder Gfi1beGFP/+ Mäusen entnommen wurden.

Zellsortierung

Für plattenbasierte scRNA-seq-Experimente wurde ein FACS-Gerät (Astrios) verwendet, um eine einzelne Zelle in jede Vertiefung einer 96-Well-PCR-Platte zu sortieren, die 5 μl TCL-Puffer mit 1 % 2-Mercaptoethanol enthält. Für die Isolierung von EpCAM+ wurden die Zellen auf 7AAD- (Life Technologies), CD45- (eBioscience), CD31- (eBioscience), TER-119- (eBioscience) und EpCAM+ (eBioscience) angefärbt; für spezifische Epithelzellen wurde auch auf CD24+/- (eBioscience) und c-Kit+/- (eBioscience) gefärbt. Zur Anreicherung spezifischer Epithelzellpopulationen des Darms wurden Zellen aus Lgr5-GFP-Mäusen isoliert, mit den oben genannten Antikörpern angefärbt und auf GFP-high (Stammzellen), GFP-low (TAs), GFP-/CD24+/c-Kit+/- (sekretorische Linien) oder GFP-/CD24-/EpCAM+ (Epithelzellen) gated. Um eine bessere Paneth-Zell-Wiederfindung zu erreichen, haben wir höhere Parameter für die seitliche Streuung und die Vorwärtsstreuung in Kombination mit CD24+/c-Kit+ zugelassen, um die Paneth-Zell-Wiederfindung in EpCAM+-Zellen zu überprüfen. Für die Tuft-2-Isolierung wurden Epithelzellen von drei verschiedenen Mäusen wie oben beschrieben angefärbt, jedoch unter Verwendung von EpCAM+/CD45+ zur Sortierung von 2.000 Einzelzellen. Wir verwendeten ein nachsichtiges Sortiergatter, um sicherzustellen, dass wir eine ausreichende Anzahl dieser seltenen Tuft-2-Zellen erhielten, was zu einer höheren Kontaminationsrate von T-Zellen führte, die wir in unserer Einzelzellanalyse durch unüberwachtes Clustering entfernten.

Für die Sortierung von scRNA-seq in voller Länge wurde die 96-Well-Platte mit einem Microseal F dicht verschlossen und 1 Minute lang bei 800 g zentrifugiert. Die Platte wurde sofort auf Trockeneis eingefroren und bei -80 °C aufbewahrt, bis sie für das Lysat-Clean-up bereit war. Die Zellen der Bulk-Population wurden in ein Eppendorf-Röhrchen mit 100 μl TCL-Lösung mit 1 % 2-Mercaptoethanol sortiert und bei -80 °C gelagert.

Für die tröpfchenbasierte scRNA-seq wurden die Zellen mit denselben Parametern wie für die plattenbasierte scRNA-seq sortiert, jedoch in ein Eppendorf-Röhrchen mit 50 μl 0,4 % BSA-PBS und auf Eis gelagert, bis sie zur GemCode-Plattform für Einzelzellen weitergeleitet wurden.

Plattenbasierte scRNA-seq

Einzelzellen. Die Bibliotheken wurden nach einem modifizierten SMART-Seq2-Protokoll16 erstellt. Kurz gesagt, wurde das RNA-Lysat mit RNAClean XP Beads (Agencourt) gereinigt, gefolgt von einer reversen Transkription mit Maxima Reverse Transcriptase (Life Technologies) und einer Ganztranskriptionsamplifikation (WTA) mit KAPA HotStart HIFI 2 × ReadyMix (Kapa Biosystems) für 21 Zyklen. Die WTA-Produkte wurden mit Ampure XP-Beads (Beckman Coulter) gereinigt, mit dem Qubit dsDNA HS Assay Kit (ThermoFisher) quantifiziert und mit einem hochempfindlichen DNA-Chip (Agilent) ausgewertet. RNA-seq-Bibliotheken wurden aus gereinigten WTA-Produkten mit dem Nextera XT DNA Library Preperation Kit (Illumina) erstellt. Auf jeder Platte wurden die Populations- und Nicht-Zell-Kontrollen nach der gleichen Methode wie für die Einzelzellen bearbeitet. Die Bibliotheken wurden auf einem Illumina NextSeq 500 sequenziert.

Massenproben. Bulk-Populationsproben wurden durch Extraktion von RNA mit dem RNeasy Plus Micro Kit (Qiagen) gemäß den Empfehlungen des Herstellers verarbeitet und dann mit dem modifizierten SMART-Seq2-Protokoll nach der Lysatreinigung, wie oben beschrieben, fortgesetzt.

Tropfen-basierte scRNA-Seq

Einzelzellen wurden mit der GemCode Single Cell Platform unter Verwendung der GemCode Gel Bead, Chip und Library Kits (10X Genomics, Pleasanton) gemäß dem Herstellerprotokoll verarbeitet. Kurz gesagt, wurden einzelne Zellen in 0,4% BSA-PBS sortiert. Jedem Kanal wurden 6.000 Zellen mit einer durchschnittlichen Wiederfindungsrate von 1.500 Zellen hinzugefügt. Die Zellen wurden dann im GemCode-Gerät in Gel-Beads in Emulsion aufgeteilt, wo die Zelllyse und die barcodierte reverse Transkription der RNA stattfand, gefolgt von Amplifikation, Scherung und 5′-Adapter- und Probenindexanbindung. Die Bibliotheken wurden auf einem Illumina NextSeq 500 sequenziert.

Immunofluoreszenz und smFISH

Immunofluoreszenz. Die Färbung des Dünndarmgewebes wurde wie zuvor beschrieben34 durchgeführt. Kurz gesagt, wurde das Gewebe 14 Stunden lang in Formalin fixiert, in Paraffin eingebettet und in 5 μm dicke Schnitte geschnitten. Die Schnitte wurden mit Standardverfahren entparaffiniert, über Nacht bei 4 °C mit primären Antikörpern und anschließend 30 Minuten lang bei Raumtemperatur mit sekundären Antikörpern inkubiert. Die Objektträger wurden mit Slowfade Mountant + DAPI (Life Technologies, S36964) montiert und versiegelt.

smFISH. Ein RNAScope Multiplex Flourescent Kit (Advanced Cell Diagnostics) wurde gemäß den Empfehlungen des Herstellers mit den folgenden Änderungen verwendet. Die Kochzeit für das Target Retrieval wurde auf 12 Minuten und die Inkubation mit Protease IV bei 40 °C auf 8 Minuten eingestellt. Die Objektträger wurden mit Slowfade Mountant+DAPI (Life Technologies, S36964) aufgezogen und versiegelt.

Kombinierte Immunofluoreszenz und smFISH. Dazu wurde zunächst smFISH wie oben beschrieben durchgeführt, mit folgenden Änderungen. Nach Amp 4 wurden die Gewebeschnitte in Waschpuffer gewaschen, über Nacht bei 4 °C mit primären Antikörpern inkubiert, dreimal in 1× TBST gewaschen und dann 30 Minuten bei Raumtemperatur mit sekundären Antikörpern inkubiert. Die Objektträger wurden mit Slowfade Mountant + DAPI (Life Technologies, S36964) eingeklebt und versiegelt.

Bildanalyse

Die Bilder der Gewebeschnitte wurden mit einem konfokalen Mikroskop Fluorview FV1200 unter Verwendung von Kalman und sequenzieller Laseremission aufgenommen, um Rauschen und Signalüberlappung zu reduzieren. Mit der konfokalen Software FV10-ASW 3.1 Viewer wurden zu jedem Bild Maßstabsbalken hinzugefügt. Die Bilder wurden mit der Software Image J45 überlagert und visualisiert.

Antikörper und Sonden

Darmorganoidkulturen

Nach der Isolierung der Krypten wurde die Einzelzellsuspension in Matrigel (BD Bioscience) mit 1 μM Jagged-1-Peptid (Ana-Spec) resuspendiert. Etwa 300 Krypten, eingebettet in 25 μl Matrigel, wurden in jede Vertiefung einer 24-Well-Platte ausgesät. Nach dem Aushärten wurde das Matrigel in 600 μl Kulturmedium (Advanced DMEM/F12, Invitrogen) mit Streptomycin/Penicillin und Glutamatax inkubiert und mit EGF (100 ng ml-1, Peprotech), R-Spondin-1 (600 ng ml-1, R&D) ergänzt, Noggin (100 ng ml-1, Prepotech), Y-276432 Dihydrochlorid-Monohydrat (10 μM, Tochris), N-Acetyl-1-Cystein (1 μM, Sigma-Aldrich), N2 (1X, Life Technologies), B27 (1X, Life Technologies) und Wnt3A (25 ng ml-1, R&D Systems). Frisches Medium wurde am dritten Tag ersetzt, und die Organoide wurden durch Dissoziation mit TrypLE passagiert und am sechsten Tag in neuem Matrigel im Verhältnis 1:3 resuspendiert. Für ausgewählte Experimente wurden die Organoide zusätzlich mit RANKL (100 ng ml-1, Biolegends) behandelt. Die behandelten Organoide wurden dissoziiert und mit beiden Methoden der scRNA-seq unterzogen.

Quantitative PCR

Die cDNA von 16 ganztranskriptom-amplifizierten Einzelzellen von tuft-1, tuft-2 und random EpCam+ aus den scRNA-seq-Platten mit voller Länge wurde für die relative qPCR verwendet. Die Genexpression wurde mittels quantitativer Echtzeit-PCR auf einem LightCycler 480 Instrument II (Roche) unter Verwendung des LightCycler 480 SYBR green mix (Roche) mit den folgenden Primersätzen analysiert: HPRT1-F, GTTAAGCAGTACAGCCCCAAA; HPRT1-R, AGGGCATATCCAACAACAAACTT; UBC-F, CAGCCGTATATCTTCCCAGACT; UBC-R, CTCAGAGGGATGCCAGTAATCTA; tslp-F, TACTCTCAATCCTATCCCTGGCTG; Tlsp-R, CCATTTCCTGAGTACCGTCATTTC; Alpi-F, TCCTACACCTCCATTCTATGG, Alpi-R, CCGCCTGCTGCTTGTAG; Dclk1-F, GGGTGAGAACCATCTACACCATC; Dclk1-R, CCAGCTTCTTAAAGGGCTCGAT. Die qPCR-Primer wurden für eine Exon-Exon-Grenze in allen Transkripten entworfen.

Computergestützte Analyse

Vorverarbeitung von tröpfchenbasierten scRNA-seq-Daten. De-Multiplexing, Alignment auf das mm10-Transkriptom und UMI (Unique Molecular Identifier)-Collapsing wurden mit dem Cellranger-Toolkit (Version 1.0.1) von 10X Genomics durchgeführt. Für jede Zelle quantifizierten wir die Anzahl der Gene, für die mindestens ein Read kartiert wurde, und schlossen dann alle Zellen mit weniger als 800 entdeckten Genen aus. Die Expressionswerte Ei,j für das Gen i in der Zelle j wurden berechnet, indem die UMI-Zählungen für das Gen i durch die Summe der UMI-Zählungen in der Zelle j geteilt wurden, um die Unterschiede in der Abdeckung zu normalisieren, und dann mit 10.000 multipliziert wurden, um TPM-ähnliche Werte zu erhalten, und schließlich log2(TPM + 1) berechnet wurde. Die Batch-Korrektur wurde mit ComBat46 durchgeführt, wie im R-Paket sva47 implementiert, unter Verwendung des standardmäßigen parametrischen Anpassungsmodus. Die Ausgabe war eine korrigierte Expressionsmatrix, die als Eingabe für die weitere Analyse verwendet wurde.

Die Auswahl der variablen Gene erfolgte durch Anpassung eines verallgemeinerten linearen Modells an die Beziehung zwischen dem quadrierten Variationskoeffizienten und dem mittleren Expressionsniveau im logarithmischen Raum und die Auswahl von Genen, die signifikant (P < 0,05) von der angepassten Kurve abwichen48.

Vorverarbeitung der SMART-Seq2 scRNA-seq-Daten. BAM-Dateien wurden mit dem von Illumina bereitgestellten Softwarepaket Bcl2Fastq v2.17.1.14 in zusammengeführte, de-multiplexierte FASTQs konvertiert. Paired-End-Reads wurden dem UCSC mm10-Maus-Transkriptom mit Bowtie49 mit den Parametern „-q –phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6“ zugeordnet, was die Ausrichtung von Sequenzen mit einer Fehlpaarung ermöglicht. Die Expression von Genen wurde anhand von TPM-Werten quantifiziert, die von RSEM50 v1.2.3 im Paired-End-Modus berechnet wurden. Für jede Zelle quantifizierten wir die Anzahl der Gene, für die mindestens ein Read kartiert wurde, und schlossen dann alle Zellen aus, bei denen entweder weniger als 3.000 Gene erkannt wurden oder die eine Transkriptomkartierung von weniger als 40 % aufwiesen. Anschließend identifizierten wir hochvariable Gene wie oben beschrieben.

Dimensionalitätsreduktion mittels PCA und t-SNE. Wir schränkten die Expressionsmatrix auf die oben erwähnten Teilmengen variabler Gene und hochwertiger Zellen ein und zentrierten und skalierten die Werte, bevor wir sie in die Hauptkomponentenanalyse (PCA) eingaben, die mit der R-Funktion prcomp aus dem stats-Paket für den SMART-seq2-Datensatz implementiert wurde. Für den Droplet-basierten Datensatz verwendeten wir eine randomisierte Annäherung an die PCA, die mit der Funktion rpca aus dem R-Paket rsvd implementiert wurde, wobei der Parameter k auf 100 gesetzt wurde. Diese Näherung mit niedrigem Rang wurde verwendet, weil sie für sehr breite Matrizen um mehrere Größenordnungen schneller zu berechnen ist. Da viele Hauptkomponenten nur einen geringen Teil der Varianz erklären, kann das Signal-Rausch-Verhältnis durch die Auswahl einer Teilmenge von n „signifikanten“ Hauptkomponenten erheblich verbessert werden. Nach der PCA wurden signifikante Hauptkomponenten mithilfe des Permutationstests51 ermittelt, der mit der Funktion permutationPA aus dem R-Paket jackstraw implementiert wurde. Dieser Test identifizierte 13 und 15 signifikante Hauptkomponenten in den 10X- und SMART-Seq2-Datensätzen von Abb. 1b bzw. Abb. 2a (erweiterte Daten). Nur die Scores dieser signifikanten Hauptkomponenten wurden als Input für die weitere Analyse verwendet.

Zur Visualisierung wurde die Dimensionalität der Datensätze mit Hilfe der Barnes-hut-Version von t-SNE52,53 weiter reduziert. Dies wurde mit der Funktion Rtsne aus dem R-Paket Rtsne mit 20.000 Iterationen und einer Perplexitätseinstellung implementiert, die je nach Größe des Datensatzes zwischen 10 und 30 variierte.

Identifizierung von Zelldifferenzierungstrajektorien mithilfe von Diffusionskarten

Vor der Dimensionalitätsreduktion mit Diffusionskarten wählten wir hochvariable Gene in den Daten wie folgt aus. Zunächst haben wir ein Nullmodell für die Grundvariabilität der Genexpression zwischen den Zellen in den Daten erstellt, wobei wir eine Potenzgesetz-Beziehung zwischen dem Variationskoeffizienten und dem Mittelwert der UMI-Zahlen aller exprimierten Gene verwendet haben, ähnlich wie in früheren Arbeiten54. Anschließend berechneten wir für jedes Gen die Differenz zwischen dem Wert seines beobachteten Variationskoeffizienten und dem durch das Nullmodell erwarteten Wert (CVdiff). Das Histogramm von CVdiff wies einen „fetten“ Schwanz auf. Wir berechneten den Mittelwert μ und die Standardabweichung σ dieser Verteilung und wählten alle Gene aus, bei denen der CVdiff > μ + 1,67σ betrug, was 761 Gene für die weitere Analyse ergab.

Wir haben die Dimensionalität mit Hilfe des Diffusionskartenansatzes22 reduziert. Kurz gesagt, wurde eine Zell-Zell-Übergangsmatrix mit einem Gauß-Kernel berechnet, wobei die Kernelbreite an die lokale Nachbarschaft jeder Zelle angepasst wurde55. Diese Matrix wurde nach der Normalisierung in eine Markov-Matrix umgewandelt. Die rechten Eigenvektoren vi (i = 0, 1, 2, …) dieser Matrix wurden berechnet und in der Reihenfolge der abnehmenden Eigenwerte λi (i = 0, 1, 2, …) sortiert, nachdem der „oberste“ Eigenvektor v0, der λ0 = 1 entspricht (was die Normalisierungsbedingung der Markov-Matrix widerspiegelt), ausgeschlossen wurde. Die übrigen Eigenvektoren vi (i = 1, 2, …) definieren die Einbettung der Diffusionskarte und werden als Diffusionskomponenten (DCk, k = 1, 2, …) bezeichnet. Wir stellten eine spektrale Lücke zwischen λ4 und λ5 fest und behielten daher DC1-DC4 sowohl für den ursprünglichen Datensatz (Erweiterte Daten, Abb. 4) als auch für die aus verschiedenen Darmregionen extrahierten Daten (Abb. 2c) bei.

Entfernen von kontaminierenden Immunzellen und Doubletten

Obwohl die Zellen vor der Sequenzierung mit EpCAM sortiert wurden, wurde im 10X-Datensatz eine kleine Anzahl kontaminierender Immunzellen beobachtet. Diese 264 Zellen wurden in einer ersten Runde des unüberwachten Clustering (dichtebasiertes Clustering der t-SNE-Karte unter Verwendung von dbscan56 aus dem R-Paket fpc) entfernt, da sie einen extrem deutlichen Cluster bildeten. Für den SMART-Seq2-Datensatz waren mehrere Zellen Ausreißer in Bezug auf die Bibliothekskomplexität, die möglicherweise mehr als einer einzelnen Zelle pro Sequenzierbibliothek entsprechen könnten („Doubletten“). Diese Zellen wurden dann entfernt, indem das oberste 1 %-Quantil der Verteilung der pro Zelle entdeckten Gene berechnet und alle Zellen in diesem Quantil entfernt wurden.

Cluster-Analyse

Um einzelne Zellen nach ihrer Expression zu clustern, verwendeten wir einen unüberwachten Clustering-Ansatz, der auf dem Infomap-Graph-Clustering-Algorithmus9 basiert und sich an Ansätzen für Einzelzell-CyTOF-Daten57 und scRNA-seq10 orientiert. Kurz gesagt, wir konstruierten einen k-Nächste-Nachbarn-Graphen auf den Daten, indem wir für jedes Zellpaar den euklidischen Abstand zwischen den Werten der signifikanten Hauptkomponenten verwendeten, um k nächste Nachbarn zu identifizieren. Der Parameter k wurde so gewählt, dass er mit der Größe des Datensatzes vereinbar ist. Konkret wurde k für den tropfenbasierten Datensatz mit 7 216 Zellen (Abb. 1b) und für den SMART-Seq2-Datensatz mit 1 522 Zellen (Erweiterte Daten, Abb. 2a) auf 200 bzw. 80 festgelegt. RANKL-behandelte Organoide enthielten 5.434 Zellen und k wurde auf 200 gesetzt; der Salmonellen- und H. polygyrus-Datensatz enthielt 9.842 Zellen und k wurde auf 500 gesetzt. Für Clusteranalysen innerhalb der Zelltypen, insbesondere der enteroendokrinen und der Büschelzellen-Untergruppen, verwendeten wir die Pearson-Korrelationsdistanz anstelle der euklidischen Distanz und setzten k = 15, k = 30 und k = 40 für die enteroendokrinen Subtypen (533 Zellen) und für die 166 bzw. 102 Büschelzellen in den 10X- und SMART-Seq2-Datensätzen. Der Nearest-Neighbour-Graph wurde mit der Funktion nng aus dem R-Paket cccd berechnet. Der k-Nächste-Nachbar-Graph wurde dann als Eingabe für Infomap9 verwendet, das mit der Funktion infomap.community aus dem R-Paket igraph implementiert wurde.

Die entdeckten Cluster wurden Zelltypen oder Zwischenzuständen zugeordnet, wobei bekannte Marker für Subtypen von Darmepithelzellen verwendet wurden. (Erweiterte Daten: Abb. 1g, Erweiterte Daten: Abb. 2a). Für die Subanalyse der enteroendokrinen Zellen (EEC) (Abb. 3) wurde jede Gruppe von EEC-Vorläufer-Clustern mit durchschnittlichen paarweisen Korrelationen zwischen signifikanten Hauptkomponenten-Scores von r > 0,85 zusammengeführt, was zu vier Clustern führte. Wir bezeichneten diese vier Cluster als Vorläufer ‚A‘ auf der Grundlage hoher Ghrl-Konzentrationen oder als Vorläufer (früh), (mittel) oder (spät) (in dieser Reihenfolge) auf der Grundlage abnehmender Konzentrationen von Stamm- (Slc12a2, Ascl2, Axin2) und Zellzyklusgenen und zunehmender Konzentrationen bekannter EEC-Regulationsfaktoren (Neurod1, Neurod2 und Neurog3) (Erweiterte Daten, Abb. 5c, Ergänzende Tabelle 6). Für den SMART-Seq2-Datensatz wurden zwei Cluster mit hohen Konzentrationen von Stammzellmarkergenen (Erweiterte Daten, Abb. 2a) zu einem „Stamm“-Cluster und zwei weitere Cluster zu einem „TA“-Cluster zusammengeführt.

Bei der Clusteranalyse des Datensatzes des Follikel-assoziierten Epithels mit 4.700 Zellen waren die Mikrofaltenzellen sehr selten (0,38 %), weshalb die ClusterDP-Methode58 zu ihrer Identifizierung verwendet wurde, da sie empirisch besser abschnitt als der k-nearest-neighbour-Graph-Algorithmus für diesen Datensatz. Wie bei den k-nearest-neighbour-Methoden wurde ClusterDP mit signifikanten (P < 0,05) Hauptkomponenten-Scores (in diesem Fall 19) als Eingabe ausgeführt und mit den Funktionen findClusters und densityClust aus dem R-Paket densityClust mit den Parametern rho = 1.1 und delta = 0.25.

Extrahieren seltener Zelltypen für die weitere Analyse

Die anfängliche Clusterbildung des gesamten Darmdatensatzes (7.216 Zellen; Abb. 1b) zeigte einen Cluster aus 310 EEC-Zellen und 166 Büschelzellen. Die Büschelzellen wurden für die Subanalyse (Abb. 4a, b) unverändert übernommen, während die EEC-Zellen mit einem zweiten Cluster von 239 EEC-Zellen kombiniert wurden, die im regionalen Datensatz identifiziert wurden (Abb. 2a, rechts), was insgesamt 549 EEC-Zellen ergab. Eine Gruppe von 16 Zellen exprimierte die EEC-Marker Chga und Chgb gemeinsam mit Markern von Paneth-Zellen, einschließlich Lyz1, Defa5 und Defa22, und wurde daher als Dubletten interpretiert und aus der Analyse entfernt, so dass 533 EEC-Zellen übrig blieben, die die Grundlage für die Analyse in Abb. 3 bildeten. Für den Vergleich der Expressionsprofile von Enterozyten aus dem proximalen und distalen Dünndarm (Abb. 2b) wurden die 1.041 Enterozyten verwendet, die aus 11.665 Zellen des regionalen Datensatzes (Abb. 2a) identifiziert wurden.

Definieren von Zelltyp-Signaturen

Um maximal spezifische Gene für Zelltypen zu identifizieren, führten wir Tests auf differentielle Expression zwischen jedem Paar von Clustern für alle möglichen paarweisen Vergleiche durch. Anschließend wurden für ein bestimmtes Cluster die mutmaßlichen Signaturgene anhand des maximalen FDR-Q-Werts gefiltert und nach dem minimalen log2(fold change) geordnet. Die minimale fache Veränderung und der maximale Q-Wert stellen die schwächste Effektgröße über alle paarweisen Vergleiche hinweg dar; es handelt sich also um ein strenges Kriterium. Die in Abb. 1c, Abb. 2b (Erweiterte Daten), Abb. 8e (Erweiterte Daten) und den ergänzenden Tabellen 2-4 und 8 gezeigten Zelltyp-Signaturgene wurden unter Verwendung einer maximalen FDR von 0,05 und einer minimalen log2(fold change) von 0,5 ermittelt. Im Fall der Signaturen für post-mitotische Zelltypen überschritten alle Gene diesen Schwellenwert sowohl in den 3′- (Abb. 1c) als auch in den Volldaten (Erweiterte Daten, Abb. 2b) Datensätzen.

Im Fall der Signaturgene für Subtypen innerhalb der Zelltypen (Abb. 3b, Abb. 4b, Erweiterte Daten, Abb. 7b) wurde ein kombinierter P-Wert für die Subtypen ermittelt. 7b) wurde ein kombinierter P-Wert (über die paarweisen Tests) für die Anreicherung mit Hilfe der Fisher-Methode berechnet – ein milderes Kriterium als der maximale P-Wert – und ein maximaler FDR-Q-Wert von 0,01 wurde verwendet, zusammen mit einem Grenzwert von mindestens log2(fold change) von 0,25 für Büschelzell-Subtypen (Abb. 4b, Erweiterte Daten Abb. 7b, Ergänzende Tabelle 7) und von 0,1 für EEC-Subtypen (Abb. 3b, Ergänzende Tabelle 6). Alle Gene in der Büschelzellsignatur übertrafen diesen Grenzwert sowohl in den 3′- (Abb. 4b) als auch in den Volllängen-Datensätzen (Erweiterte Daten, Abb. 7b), während die EEC-Subtyp-Signaturen nur mit 3′ definiert wurden. Aufgrund der geringen Zellzahlen (n = 18) wurde der kombinierte Fisher-P-Wert auch für die in vivo-Mikrofalzzellsignatur verwendet, mit einem FDR-Cut-off von 0,001 (Abb. 5d, ergänzende Tabelle 8). Die Markergene wurden nach der minimalen log2(fold change) gereiht. Differentialexpressionstests wurden mit dem Mann-Whitney-U-Test (auch bekannt als Wilcoxon-Rangsummentest) durchgeführt, der mit der R-Funktion wilcox.test implementiert wurde. Für die Infektionsexperimente (Abb. 6) verwendeten wir ein zweiteiliges „Hürden“-Modell, um sowohl die technische Qualität als auch die Maus-zu-Maus-Variation zu kontrollieren. Dieses Modell wurde mit dem R-Paket MAST59 implementiert, und die P-Werte für die differentielle Expression wurden mit dem Likelihood-Ratio-Test berechnet. Die Korrektur der multiplen Hypothesentests erfolgte durch Kontrolle der FDR60 mit Hilfe der R-Funktion p.adjust.

Bewertung von Zellen mit Hilfe von Signatur-Gensätzen

Um einen Score für einen spezifischen Satz von n Genen in einer bestimmten Zelle zu erhalten, wurde ein „Hintergrund“-Gensatz definiert, um Unterschiede in der Sequenzierungsabdeckung und der Bibliothekskomplexität zwischen den Zellen zu kontrollieren, ähnlich wie in Ref. 12. Der Hintergrund-Gensatz wurde so ausgewählt, dass er den interessierenden Genen in Bezug auf das Expressionsniveau ähnlich ist. Konkret wurden die 10n nächsten Nachbarn im zweidimensionalen Raum, der durch die mittlere Expression und die Entdeckungshäufigkeit in allen Zellen definiert ist, ausgewählt. Der Signatur-Score für diese Zelle wurde dann definiert als die mittlere Expression der n Signaturgene in dieser Zelle abzüglich der mittleren Expression der 10n Hintergrundgene in dieser Zelle.

Schätzungen der Zelltyp-Stichprobenhäufigkeiten

Für jeden Zelltyp wird die Wahrscheinlichkeit der Beobachtung von mindestens n Zellen in einer Stichprobe der Größe k anhand der kumulativen Verteilungsfunktion eines negativen Binoms NBcdf(k, n, p) modelliert, wobei p die relative Häufigkeit dieses Zelltyps ist. Für m Zelltypen mit demselben Parameter p ist die Gesamtwahrscheinlichkeit, jeden Typ mindestens n Mal zu sehen, NBcdf(k; n, p)m. Eine solche Analyse kann mit benutzerdefinierten Parametern unter http://satijalab.org/howmanycells durchgeführt werden.

EEC-Dendrogramm

Durchschnittliche Expressionsvektoren wurden für alle 12 EEC-Subset-Cluster unter Verwendung von log2(TPM + 1)-Werten berechnet und auf die Untergruppe von 1.361 Genen beschränkt, die als signifikant variabel zwischen den EEC-Subsets identifiziert wurden (P < 0,05), wie oben beschrieben. Die durchschnittlichen Expressionsvektoren einschließlich dieser Gene wurden mit dem R-Paket pvclust (Spearman-Distanz, Ward.D2-Clustermethode) hierarchisch geclustert, das Bootstrap-Vertrauensschätzungen für jeden Dendrogrammknoten als empirischen P-Wert über 100.000 Versuche liefert (Extended Data Abb. 6a).

Zelltyp-spezifische Transkriptionsfaktoren, GPCRs und Leucin-reiche Repeat-Proteine

Eine Liste aller Gene, die als Transkriptionsfaktoren in Mäusen identifiziert wurden, wurde von AnimalTFDB61 bezogen. Die Gruppe der GPCRs wurde aus der UniProt-Datenbank (http://www.uniprot.org/uniprot/?query=family%3A%22g+protein+coupled+receptor%22+AND+organism%3A%22Mouse+%5B10090%5D%22+AND+reviewed%3Ayes&sort=score) entnommen. Funktionelle Annotationen für jedes Protein (Erweiterte Daten, Abb. 2d) wurden von der British Pharmacological Society (BPS) und der International Union of Basic and Clinical Pharmacology (IUPHAR) (http://www.guidetopharmacology.org/GRAC/GPCRListForward?class=A) bezogen. Die Liste der Leucin-reichen Repeat-Proteine wurde aus ref. 62. Für die Zuordnung von Gennamen von Menschen zu Mäusen wurden die Orthologe von Menschen und Mäusen von Ensembl (neueste Version 86; http://www.ensembl.org/biomart/martview) und die Synonyme von Menschen- und Mausgenen vom NCBI (ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/) heruntergeladen. Für jedes menschliche Leucin-reiche Repeat-Gen wurden alle menschlichen Synonyme mit Hilfe der Ortholog-Liste auf das orthologe Gen in der Maus abgebildet, und die Gennamen der Maus wurden mit Hilfe der Synonym-Liste auf die in den Einzelzelldaten enthaltenen Gene abgebildet.

Dann wurden zelltyp-angereicherte Transkriptionsfaktoren, GPCRs und Leucin-reiche Repeat-Proteine identifiziert, indem die Liste der in jedem Zelltyp angereicherten Gene mit den oben definierten Listen der Transkriptionsfaktoren, GPCRs und Leucin-reichen Repeat-Proteine gekreuzt wurde. Zelltyp-angereicherte Gene wurden unter Verwendung des SMART-Seq2-Datensatzes als solche mit einer minimalen log2(fold change) von 0 und einer maximalen FDR von 0,5 definiert, wobei maximal 10 Gene pro Zelltyp in den erweiterten Daten (Abb. 2e, f) erhalten blieben (vollständige Listen sind in der ergänzenden Tabelle 5 enthalten). Darüber hinaus wurde ein umfangreicheres Panel von zelltypspezifischen GPCRs identifiziert (Erweiterte Daten, Abb. 2d), indem ein milderer Schwellenwert gewählt wurde. Dies wurde erreicht, indem jeder Zelltyp mit allen anderen Zellen verglichen wurde, anstatt der im vorherigen Abschnitt beschriebenen paarweisen Vergleiche, und alle GPCR-Gene ausgewählt wurden, die unterschiedlich exprimiert wurden (FDR < 0,001).

Prüfung auf Veränderungen der Zelltyp-Proportionen

Wir modellieren die entdeckte Anzahl jedes Zelltyps in jeder analysierten Maus als eine zufällige Zählvariable unter Verwendung eines Poisson-Prozesses. Die Entdeckungsrate wird dann modelliert, indem die Gesamtzahl der Zellen, die bei einer bestimmten Maus profiliert wurden, als Offset-Variable angegeben wird, wobei der Zustand jeder Maus (Behandlung oder Kontrolle) als Kovariate angegeben wird. Das Modell wurde mit dem R-Befehl glm aus dem Paket stats angepasst. Der P-Wert für die Signifikanz des durch die Behandlung hervorgerufenen Effekts wurde mit einem Wald-Test auf den Regressionskoeffizienten bewertet.

Für die Bewertung der Signifikanz der räumlichen Verteilungen der EEC-Untergruppen (Abb. 3e) umfasste der Vergleich mehr als zwei Gruppen. Unsere Nullhypothese war insbesondere, dass der Anteil jeder EEC-Untergruppe, die in den drei Darmregionen (Duodenum, Jejunum und Ileum) nachgewiesen wurde, gleich war. Um diese Hypothese zu testen, verwendeten wir die Varianzanalyse (ANOVA) mit einem χ2-Test auf die oben beschriebene Anpassung des Poisson-Modells, implementiert mit der anova-Funktion des stats-Pakets.

Anreicherung des Gensatzes und Gen-Ontologie-Analyse

Die Gen-Ontologie-Analyse wurde mit dem goseq R-Paket63 durchgeführt, wobei signifikant differentiell exprimierte Gene (FDR < 0.05) als Zielgene und alle Gene, die mit log2(TPM + 1) > 3 in mindestens zehn Zellen exprimiert wurden, als Hintergrund verwendet.

Datenverfügbarkeit

Codeverfügbarkeit