Live-seq enables temporal transcriptomic recording of single cells – Nature

Biological materials

RAW264.7, 293T and HeLa cells were obtained from ATCC. RAW264.7 cells with Tnf-mCherry reporter and relA-GFP fusion protein (RAW-G9 clone) were kindly provided by I.D.C. Fraser (National Institutes of Health). The IBA cell line derived from the stromal vascular fraction of interscapular brown adipose tissue of young male mice (C57BL6/J) was kindly provided by C. Wolfrum’s laboratory (ETH Zurich). Primary ASPCs were isolated from subcutaneous fat tissue of female C57BL/6J mice at age between 7 and 10 weeks as previously described⁵⁰. The cell lines have not been authenticated or tested for mycoplasma contamination. All these cells were cultured in high-glucose DMEM medium (Life Technologies) supplemented with 10% foetal bovine serum (FBS) (Gibco) and 1× penicillin/streptomycin solution (Life Technologies) in a 5% CO₂ humidified atmosphere at 37 °C and maintained at less than 80% confluence before passaging.

Oligos

Biotinylated-olig0-dT30VN (IDT):

/5Biosg/AAGCAGTGGTATCAACGCAGAGTACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTVN.

TSO (Exiqon): AAGCAGTGGTATCAACGCAGAGTACATrGrG+G.

Iso-TSO (Exiqon): iCiGiCAAGCAGTGGTATCAACGCAGAGTrGrG+G.

Biotinylated-TSO (Exiqon): /5Biosg/AAGCAGTGGTATCAACGCAGAGTACATrGrG+G.

Hairpin-TSO (Exiqon):

GGCGGCGCGCGCGCGCCGCCAAGCAGTGGTATCAACGCAGAGTACATrGrG+G.

Biotinylated-ISPCR oligo (IDT): /5Biosg/AAGCAGTGGTATCAACGCAGAGT.

Hairpin-TSO was annealed in IDT Duplex Buffer before use.

FluidFM setup

A FluidFM system composed of a FlexAFM-near-infrared scan head and a C3000 controller driven by the EasyScan2 software (Nanosurf), and a digital pressure controller unit (ranging from −800 to +1,000 mbar) operated by a digital controller software (Cytosurge) was used. A syringe pressure kit with a three-way valve (Cytosurge) was used in addition to the digital pressure controller to apply underpressure lower than −800 mbar and overpressure higher than +1,000 mbar. FluidFM Rapid Prototyping probes with a microchannel height of 800 nm were obtained from Cytosurge. A 400-nm wide triangular aperture was custom-milled by a focused ion beam near the apex of the pyramidal tips and imaged by scanning electron microscopy as previously described³. The FluidFM probes were plasma treated for 1 min (Plasma Cleaner PDG-32G, Harrick Plasma) and coated overnight with SL2 Sigmacote (Sigma-Aldrich) in a vacuum desiccator³.

Optical microscopy setup

The FluidFM scan head was mounted on an inverted AxioObserver microscope equipped with a temperature-controlled incubation chamber (Zeiss), and coupled to a spinning disc confocal microscope (Visitron) with a Yokogawa CSU-W1 confocal unit and an EMCCD camera system (Andor). Phase-contrast and fluorescence images were acquired using ×10 and ×40 (0.6 numerical aperture, NA) objectives and a ×2 lens switcher using VisiView software (Visitron). Microscopy images were analysed using the AxioVision and ImageJ softwares.

Enhanced Smart-seq2

The overall workflow follows the same steps as Smart-seq2 (ref. ²⁵), but several modifications were introduced. RNA, Live-seq samples or single cells were transferred into a PCR tube with 4.2 μl of lysis buffer, which contains 1 μl of dNTP (10 mM each), 1 μl of biotinylated-oligo-dT30VN oligo (10 μM), 2 μl of 0.2% Triton-X-100, 0.1 μl of Recombinant ribonuclease inhibitor (40 U μl⁻¹, no. 2313, Takara) and 0.1 μl ERCC RNA Spike-In Mix (10⁷ dilution, no. 4456740, ThermoFisher). After briefly spinning, the PCR tubes/plates were heated at 72 °C for 3 min and cooled down on ice for 1 min. Then, 5.2 μl of reverse transcription mix (1.29 μl of H₂O, 0.06 μl of MgCl₂ (1 M), 2 μl of Betaine (5 M), 0.08 μl of biotinylated-TSO oligo (100 μM) (or other TSOs as listed in Supplementary Fig. 1c), 2 μl of Maxima H Mimus RT buffer (5×), 0.25 μl of recombinant ribonuclease inhibitor and 0.1 μl of Maxima H Minus Reverse Transcriptase (200 U μl⁻¹, EP0751, ThermoFisher)) were added. The reverse transcription step was performed in a thermal cycler using the following program: (1) 42 °C for 90 min; (2) 50 °C for 2 min; (3) 42 °C for 2 min, go to step (2) for four cycles; (4) 85 °C for 5 min and (5) end. The reverse transcription products were directly used for PCR by adding 15 μl of PCR mixture, which contained 12.5 μl of KAPA HiFi HotStart ReadyMix (2×, 07958935001, KAPA), 0.25 μl of Biotinylated-ISPCR oligo (10 μM) and 2.25 μl of H₂O. The PCR program involves the following steps: (1) 98 °C for 3 min; (2) 98 °C for 20 s; (3) 67 °C for 15 s; (4) 72 °C for 6 min, go to step (2) for x cycles (for Live-seq, x = 24; for scRNA-seq, x = 18); (5) 72 °C for 5 min and (6) end. The PCR products were purified twice with AMPure XP beads (A63882, Beckman Coulter) according to the manufacturer’s instructions with 0.6 volume of beads per 1 volume of PCR product. The concentration and profile of the purified DNA were measured by Qubit double-stranded DNA HS Assay Kit (no. Q32854, ThermoFisher) and fragment analyser (Advanced Analytical), respectively. Then 1 ng of DNA was used for tagmentation in the following reaction: 1 μl TAPS buffer (50 mM TAPS-NaOH, pH 8.3, 25 mM MgCl₂), 0.5 μl of Tn5 (made in-house, 12.5 μM), H₂O to 5 μl in total. The tagmentation reaction was assembled on ice, initiated by placing the mix on a thermal cycler at 55 °C for 8 min after which it was kept on ice. Then 1.5 μl of SDS (0.2%) was added to inactivate the Tn5 enzyme. PCR reagent mixture was added directly to the tagmented DNA, which contained 1.5 μl of dNTP (10 mM each), Nextera XT index primers (or other compatible indexing primers, with a final concentration of 300 nM for each primer), 10 μl of KAPA HiFi Fidelity buffer (5×), 1 μl of KAPA HiFi DNA Polymerase (1 U µl⁻¹, KK2101, KAPA) and H₂O to 50 μl in total. Then, the DNA was amplified using the following program: (1) 72 °C for 3 min; (2) 98 °C for 30 s; (3) 98 °C for 10 s; (4) 63 °C for 30 s; (5) 72 °C for 1 min, go to step (3) for ten cycles; (6) 72 °C for 5 min and (7) end. The PCR products were purified twice with AMPure XP beads (A63882, Beckman Coulter) according to the manufacturer’s instructions with 0.6 volume of beads per 1 volume of PCR product. The concentration and profile of the purified DNA were measured using a Qubit dsDNA HS Assay Kit and fragment analyser, respectively, after which the libraries could be sent for sequencing using Illumina sequencers. We thereby noted that the ‘cDNA’ yield from the negative control (0 pg input RNA) (Extended Data Fig. 1b) was mainly derived from oligo sequences, because sequencing of the respective library revealed mainly reads with a very high A/T content, which poorly aligned to the mouse/human genome (Extended Data Fig. 1e, f).

Cytoplasmic biopsies

Cell samples and reagents

For extraction experiments, the cells were seeded and cultured for at least 18 h onto 50-mm tissue-culture-treated low µ-dishes (Ibidi). Shortly before the experiment, the culture medium was replaced with 5 ml of CO₂-independent growth medium supplemented with 10% FBS, 1× penicillin/streptomycin solution and 2 mM glutamine. For the extraction of LPS-stimulated RAW cells, the culture medium was replaced with 5 ml of CO₂-independent growth medium supplemented with 10% FBS, 1× penicillin/streptomycin solution, 2 mM glutamine and 100 ng ml⁻¹ LPS, and the cells were extracted between 4 and 5 h after the addition of LPS. For the extraction of differentiating ASPCs, induction medium (ASPC adipogenic differentiation section) was replaced shortly before the experiment with 5 ml of CO₂-independent growth medium supplemented with 10% FBS, 1× penicillin/streptomycin solution, 2 mM glutamine, 1 μM dexamethasone, 0.5 mM 3-isobutyl-1-methylxanthine, 1 μM rosiglitazone and 167 nM insulin. Sampling buffer for preloading was prepared by supplementing a 0.2% solution of Triton-X 100 in nuclease-free water with 2 U μl⁻¹ recombinant RNase inhibitors (Clonetech). Lysis buffer for extract transfer was prepared as that used in the enhanced Smart-seq2 protocol. Whereas cells were maintained at 37 °C for time-lapse microscopy before and after extractions, the extraction procedures were all performed at room temperature.

Extraction process

The probe reservoir was loaded with 10 µl of mineral oil (Sigma-Aldrich) and a pressure of Δ + 1,000 mbar was applied to flow the oil into the microchannel. Once the probe microchannel was filled, the probe was shortly immersed in nuclease-free water, and then kept in air with the residual water carefully blotted off the probe holder with a kimwipe tissue. A 1.0-μl drop of sampling buffer was deposited onto an AG480F AmpliGrid (LTF Labortechnik). The cantilever was introduced into the drop using the micrometre screws to displace the atomic force microscope. Once the cantilever was located inside the drop, underpressure (−800 mbar) was applied for the suction of roughly 0.5 pl of sampling buffer into the probe. The cantilever was then withdrawn from the drop using the micrometre screws.

Next, the preloaded probe was immersed in the cell sample experimental medium, the cell to be extracted was visualized by light microscopy and the tip of the FluidFM probe was placed above the cytoplasm of the selected cell. The tip of the probe was then inserted into the cytoplasm through a forward force spectroscopy routine driven by the Z-piezo. The probe was then maintained inside the cell at constant force. Although lower forces may be sufficient, a force setpoint of 500 nN was used to ensure the full insertion of the probe aperture into all the cell types. Underpressure larger than −800 mbar was applied to aspirate the cellular content into the probe, whereby the harvested cytoplasmic fluid immediately mixed with the preloaded sampling buffer. The pressure-assisted flow of the intracellular content into the FluidFM probe was interrupted by switching the pressure back to zero. We collected cytoplasmic extracts of 1.1 pl on average, ranging from 0.1 up to 4.4 pl. The probe was then retracted out of the cell, shortly immersed in nuclease-free water and then kept in air with the residual water carefully blotted off the probe holder with a kimwipe tissue.

A 1.0-μl drop of lysis buffer was deposited onto an AG480F AmpliGrid (LTF Labortechnik). The cantilever was introduced into the drop, and overpressure (more than 1,000 mbar) was applied to release the extract. The microchannel was then rinsed three times by suction and release of lysis buffer into the probe. The 1.0-μl drop was then pipetted into a PCR tube containing an additional 3.2 μl of lysis buffer, and the solution was briefly centrifuged and stored at −80 °C until further processing.

All steps were monitored in real time by optical microscopy in brightfield. The entire procedure took roughly 15 min per sample, with around 5 min for loading the sampling buffer and approaching a selected cell, 5 min to extract the cytoplasmic biopsy and 5 min to transfer the extract into the lysis buffer and then to the PCR tube. We performed 43 experiments, collecting 10 to 20 biopsies per experiment.

Alternated human/mouse sampling

For the assessment of cross-contamination between samples that were extracted sequentially with the same probe, alternated human/mouse cell sampling was conducted, with HeLa and IBA cells seeded on separate dishes. The cells were extracted alternatively from one or the other cell type as described above, using the same FluidFM probe.

Sequential sampling of individual RAW cells

For the sequential sampling of the same RAW cell, up to 24 cells were extracted as described above, with all the cells monitored within one vision field. The first extractions were performed within a 1 h time window, 0.5 to 1.5 h before the addition of LPS. The cells were then incubated at 37 °C under time-lapse monitoring with intervals of 5 min. After 30 min, LPS was added to the dish for stimulation and the cells were further monitored at 5 min intervals for another 30 min. The time-lapse intervals were then increased to 30 min, and the cells were further monitored for 3.5 h. The temperature was then switched back to room temperature, and the same cells were extracted a second time, in a time interval of 4 to 5 h after the addition of LPS.

Sequential sampling of individual ASPC cells

For the sequential sampling of ASPCs, barcoded ASPCs (ASPC genetic barcoding section) were mixed with ASPCs at a 1:1,000 ratio and the cells were then seeded to a final density of 200,000 cells per cm² into a well with a growth area of 0.22 cm² (two-well insert, Ibidi). The cells were grown in the insert for 2 days. The insert was then removed and the growth medium exchanged for CO₂-independent medium. The entire culture area was first imaged at ×10 magnification, in brightfield and at 488 nm, and the images were assembled to create a map of the GFP-expressing, barcoded ASPCs in the entire culture area. Selected GFP-expressing ASPCs were then extracted for a first time, during a time window of 3 to 4 h. Following extraction, the medium was exchanged for induction medium (ASPC adipogenic differentiation), and the cells were incubated in a 5% CO₂ humidified atmosphere at 37 °C. After 2 days, the medium was exchanged for CO₂-independent medium supplemented with 1 μM dexamethasone, 0.5 mM 3-isobutyl-1-methylxanthine, 1 μM rosiglitazone and 167 nM insulin. The entire culture area was imaged again at ×10 magnification, in brightfield and at 488 nm, and the images assembled to create the map. The cells that had been extracted were relocalized from their respective position in the created map, and were extracted a second time. After the extractions, the medium was exchanged for the maintenance medium (ASPC adipogenic differentiation section), and the cells were incubated for 2 days in a 5% CO₂ humidified atmosphere at 37 °C. The medium was then exchanged for complete culture medium, and the cells were incubated for another 3 to 5 days with the medium exchanged every 2 days. Lipid staining was then performed using 5 μg ml⁻¹ BODIPY 558/568 (Invitrogen) for 20 min. Nucleus staining with 5 μg ml⁻¹ 4,6-diamidino-2-phenylindole (DAPI) was performed at the same time. The entire culture area was then imaged to create the map, before imaging individual barcoded cells in brightfield, at 405, 488 and 640 nm.

Molecular recording of RAW cells

For the molecular recording to predict a cell’s downstream phenotype, RAW cells within one vision field were extracted as described above, within a 1 h time window. The cells were then monitored at 37 °C under time-lapse imaging with intervals of 5 min. After 30 min, LPS was added to the dish for stimulation (0.5 to 1.5 h after extraction), and the cells were further monitored with 5 min intervals for another 30 min. The time-lapse intervals were then increased to 30 min, and the cells were further monitored for at least 8 h.

LPS stimulation of RAW cells

LPS (no. L4391-1MG) was prepared in PBS at 100 μg ml⁻¹ as a stock solution. For stimulation of RAW cells, LPS solution was added to the 5 ml of CO₂-independent medium to reach a final concentration of 100 ng ml⁻¹.

ASPC genetic barcoding

The genetic barcoding of ASPCs was performed using the pLARRY system⁵¹. LARRY Barcode Library v.1 was a gift from F. Camargo (Addgene no. 140024). The DNA was prepared directly by collecting cells directly from a Luria–Bertani agar plate rather than from a liquid culture to maximally preserve the overall library complexity, which was further confirmed by next-generation sequencing.

Barcode-bearing lentivirus was produced in 293T cells using the third-generation lentivirus packing system. The virus-containing medium was collected 48 h posttransfection. ASPCs were transduced by the virus-containing medium and fresh medium with Polybrene (final concentration at 10 μg ml⁻¹), followed by centrifugation at 1,500g for 30 min after which cells were returned to the cell incubator. The medium was exchanged 12 h later and cells were maintained at a density lower than 80% confluency.

ASPC adipogenic differentiation

To differentiate ASPCs into adipocytes, confluent cells were exposed to adipogenic induction medium (complete culture medium supplemented with 1 μM dexamethasone, 0.5 mM 3-isobutyl-1-methylxanthine, 167 nM insulin and 1 μM rosiglitazone (DMIR), all from Sigma). After 2 days, the induction medium was removed and the maintenance medium (the complete culture medium supplemented with 167 nM insulin) added. At day 4, the maintenance medium was removed and complete culture medium was added. Lipid droplets were stained using BODIPY 558/568 (Invitrogen) at 5 μg ml⁻¹ for 20 min after which cells were subjected to imaging.

Recovery of genetic barcodes from Live-seq-sampled ASPCs

The barcodes were retrieved from the cDNA of Live-seq samples. The barcode region was enriched from 1 ng cDNA by PCR using the primers BC-FOR1 (5′-ctgagcaaagaccccaacgagaa-3′) and BC-REV1 primers (5′-gctggcaactagaaggcacag-3′) and using the following program: (1) 98 °C for 30 s; (2) 98 °C for 10 s; (3) 63 °C for 30 s; (4) 72 °C for 1 min, go to step (2) for 15 cycles; (5) 72 °C for 5 min and (6) end. Then 1 μl of PCR product was subjected to a second round of PCR with For-MEDS-A (5′-tcgtcggcagcgtcagatgtgtataagagacagcgttgctaggagagaccatatg-3′) and Rev-MEDS-B primers (5′-gtctcgtgggctcggagatgtgtataagagacaggtcgacaccagtctcattcagc-3′), and using the following program: (1) 98 °C for 30 s; (2) 98 °C for 10 s; (3) 63 °C for 30 s; (4) 72 °C for 1 min, go to step (2) for 15 cycles; (5) 72 °C for 5 min and (6) end. The result PCR product was purified with AMPure XP beads (Beckman) using a 2.5 volume of beads per 1 volume of PCR product ratio, and eluted with 20 μl of nuclease-free water. Then 10 μl of the purified product was indexed using the Nextera index kit (Illumina), using the following program: (1) 98 °C for 30 s; (2) 98 °C for 10 s; (3) 63 °C for 30 s; (4) 72 °C for 1 min, go to step (2) for four cycles; (5) 72 °C for 5 min and (6) end. The resulting product was purified with AMPure XP beads (Beckman) using a 2.5 volume of beads per 1 volume of PCR product ratio, and eluted with 20 μl of nuclease-free water. The concentration was then measured using a Qubit dsDNA HS Assay Kit and subjected to sequencing in the Gene Expression Core facility at the EPFL.

For the data analysis, we created a .fasta file representing the vector sequence, replacing the barcodes by N strings. All samples were then aligned on this vector using STAR v.2.7.9a (ref. ⁵²) with default parameters and these two extra options: ‘–outFilterMultimapNmax 1’ (for removing several mapped reads) and ‘–outFilterMismatchNmax 2 –outFilterScoreMinOverLread 0 –outFilterMatchNminOverLread 0 –outFilterMatchNmin 0’ (for authorizing alignment on Ns, and a maximum of two mismatches). We then extracted the sequences from the aligned .bam file at the position of the two barcodes and counted their occurrence to find the most prevalent ones. A Hamming distance of the barcodes between cells was then calculated using the function StrDist of the DescTools package⁵³ in R using the parameter (method = ‘hamming’, mismatch = 1, gap = 1).

Determination of the extracted volumes

For each extraction, the volume of preloaded sampling buffer and the volume of cytoplasmic extract mixed with the sampling buffer were measured on brightfield images using the AxioVision software. The area occupied by the aqueous solutions confined in the cantilever was multiplied by the channel height of 0.8 μm, and the volume of the hollow pyramidal tip (90 fl) was added. To determine the volume of extracted cytoplasmic fluid, the volume of preloaded sampling buffer (FluidFM probe before extraction) was subtracted to the volume of mixed sampling buffer and extract (FluidFM probe after extraction).

Determination of the cell volumes

To quantify the whole-cell volumes, IBA and ASPC cells were dissociated by trypsinization and RAW264.7 cells with a cell scraper. The dissociated cells were then imaged in brightfield with the ×40 objective and the ×2 lens switcher, and the diameter of the rounded cells was measured from the micrographs using the AxioVision software. Three diameters were measured and averaged for each cell (n = 277 for ASPC and n = 500 for IBA and RAW). The volumes were calculated using the formula for a sphere. For the longitudinal measurements of RAW cell volumes, the areas of selected semi-adherent cells were measured on brightfield images acquired with the ×40 objective and the ×2 lens switcher at several time points before and after extraction, and cell volumes were calculated assuming a spherical cell shape.

Time-lapse monitoring of mCherry expression

For live imaging of RAW-G9 cells expressing mCherry, the Ibidi dish was covered and the temperature was maintained at 37 °C. Time-lapse images were acquired at 5 min intervals for 1 h, then at 30 min intervals overnight. For each time point, two sequential frames were acquired, in brightfield and for mCherry (561 nm laser and 609/54 nm emission filter), using the ×40 objective and the ×2 lens switcher. At the end of the time-lapse recording, 2 µl of tricolour calibration beads (Invitrogen MultiSpeck Multispectral Fluorescence Microscopy Standards Kit, Life Technologies Europe B.V.) was added to the sample and at least 30 individual beads were imaged with the 561 laser. For image quantification, all the cells that moved out of view or focus, died or overlapped with other cells were excluded from analysis in each time-lapse frame. Cell boundaries were all manually defined in Fiji, and background intensities were measured for each time point and subtracted from all intensity values. We analysed 77 extracted cells, 122 LPS-stimulated cells that were not extracted and 23 cells that were not extracted and not stimulated with LPS. The fluorescence intensity of the calibration beads with subtracted background intensity was measured using Fiji, and the average of the 30 bead intensities was used to normalize the fluorescence measurements of the cells acquired in different experiments. The area under a curve was calculated from 3 to 7.5 h post-LPS treatment. A two-sided Wilcoxon rank-sum test was used to examine the differences between conditions.

Cell viability after extraction

To evaluate the postextraction viability of ASPC (n = 33) and IBA (n = 37) cells, the cells were stained between 2 and 4 h after extraction using a LIVE/DEAD Cell Imaging Kit 488/570 (Invitrogen), and following the manufacturer’s protocol. To evaluate the postextraction viability of RAW264.7 cells, the extracted cells (n = 72) were monitored by time-lapse microscopy during around 10 h at 30 min intervals. Cells were evaluated as dead or alive on the basis of their morphology, movements and expression of mCherry in response to LPS stimulation. The postextraction viability was assessed for 42 cells extracted once before stimulation, 30 cells extracted once after LPS stimulation and ten cells extracted twice, before and after stimulation.

The volumes that were extracted ranged between 0.7 and 3.3 pl (mean 1.4 pl) for ASPCs, between 0.4 and 2.9 pl (mean 1.0 pl) for IBA cells and between 0.2 and 3.5 pl (mean 1.1 pl) for RAW cells. The viability of all the cell types was calculated as an absolute value without normalization.

scRNA-seq

Cells were trypsinized and dissociated into a single-cell suspension and then kept on ice for all the downstream processing. After passing through a 40-μm cell constrainer and DAPI staining, live singlet cells were sorted into 96-well PCR plates with 4.2 μl of lysis buffer (mentioned above). At least three wells without cells were preserved as negative controls. The plate was quickly spun and further processed using the Smart-seq2 method²⁵.

To perform scRNA-seq on cells post-Live-seq extraction, cells were first cultured in a dish containing a silicone micro-insert (Ibidi, no. 80409) at a density of around 20 cells per well. The insert was then removed just before Live-seq sampling and all the cells were subjected to Live-seq extraction as described above, during a 1 h time window. To extract all the cells in this time, the extracted cytoplasm was not preserved. Then 1 and 4 h after the middle of the extraction time window (that is, 1 ± 0.5 and 4 ± 0.5 h postextraction), the cells, along with the control cells not extracted, were collected on ice and single cells were picked using a serial dilution approach. The downstream processing followed a similar workflow as the Smart-seq2 method.

Cell cycle analyses

The cell cycle reporter vector pCSII-EF-miRFP709-hCdt(1/100) was a kind gift from V. Verkhusha (Addgene plasmid no. 80007; http://n2t.net/addgene:80007; RRID Addgene_80007). Lentivirus carrying this vector was transduced into RAW-G9 cells after which miRFP709-positive cells were sorted to enrich for transduced cells. The cells were seeded on an Ibidi dish and monitored for miRFP709 (640 nm laser and 700/75 nm emission filter) and in brightfield for 24 h at 60 min intervals with the ×40 objective and the 2× lens switcher. The experimental growth medium was then exchanged, LPS was added at a final concentration of 100 ng ml⁻¹, and the cells were monitored for another 24 h at 60 min intervals in brightfield, for miRFP709 and for mCherry (561 nm laser and 609/54 nm emission filter). Fluorescence intensities of the cells were measured as described above (Time-lapse monitoring of mCherry expression) for both fluorescent reporters. The mCherry intensities were further normalized using the beads calibration.

Nfkbia reporter analyses

The mouse Nfkbia promotor (+1,606 to −121) fragment was obtained by PCR (forward primer ttcaaaattttatcgatcagtgaaatccagaccagccgggcctac, reverse primer ggctgtgcggggctgagcgg) from mouse genomic DNA. The TagBFP CDS fragment was obtained by PCR (forward primer tgcagcctgcacccgctcagccccgcacagccACCatgagcgagctgattaaggagaac, reverse primer tgtaatccagaggttgattgtcgacgcggccgcttaattaagcttgtgccccagtttgc) from a TagBFP-bearing plasmid. A linearized lentivirus vector devoid of the EF1a promoter was obtained by PCR (forward primer gcggccgcgtcgacaatcaac, reverse primer cccggctggtctggatttcactgatcgataaaattttgaattttgtaatttgtttttgtaattc) using the pLV-vector as template (kindly provided by J. Han). All three fragments were then assembled using a Gibson Assembly Master Mix (NEB) according to the manufacturer’s instructions.

The lentivirus production from 293T and transduction into RAW-G9 cells were performed in the same way as implemented for ASPC barcoding, as described above. The resulting cells were then seeded into 96-well plates with each well containing no more than one cell to generate single clones. The single clones were then screened for their capacity to respond to LPS as well as the brightness and the broad distribution of BFP fluorescence at basal level. Cells from a selected clone were seeded in a 50-mm low-wall dish (Ibidi) and grown for 1 day at normal cell culture conditions. One culture area was then imaged with a ×40 objective and the ×2 lens switcher, at 37 °C, in brightfield and at 405 nm (BFP) and 561 nm (mCherry), at 30 min intervals for 10 h. The cells were first monitored for 1.3 h to record their basal level of Nfkbia-BFP and Tnf-mCherry, after which LPS was added to the sample and the cell were imaged for another 8.7 h to monitor their response to the stimulation. Fluorescence intensities of the cells were measured as described above (Time-lapse monitoring of mCherry expression section) for both fluorescent reporters. The mCherry and BFP intensities were further normalized using the beads calibration. The basal BFP intensity was averaged from the three first time-frames, acquired before LPS stimulation.

Live-seq and scRNA-seq data analysis

Dataset description

There were ten IBA cells and HeLa cells processed through Live-seq sampling in an alternating fashion. In addition, 588 Live-seq samples were prepared across five experimental replicates, with each of them containing both single and sequential sampling events. All the data from these 588 cells were used for the analyses presented in Figs. 2–5. Then 554 cells were processed using conventional scRNA-seq as part of three experimental replicates and the data linked to these cells are presented in Figs. 2–5. Finally, to evaluate the potential molecular perturbation of cytoplasmic sampling, scRNA-seq data of IBA cells 1 h (49 cells) and 4 h (43 cells) postextraction were generated using cells (70 cells) that were not subjected to such sampling as control.

Alignment and feature counting

Libraries were sequenced in either 75 single-read or 2 × 75 paired-end format using Nextera indexes on Illumina Nextseq500 or Hiseq4000 sequencers at the Gene Expression Core facility (at EPFL). Basecalls are performed using bcl2fastq (v.2.19.1). To keep consistency, only read 1 with 75 bp was used for further analysis. The reads were aligned to the human (hg19/GRCh38) or mouse (mm10/GRCm38) genomes using STAR (v.2.6.1c)⁵² with default settings and filtered for uniquely mapped reads. Then, the number of reads per feature (gene) was counted using HTseq (v.0.10.0)⁵⁴ with parameter ‘htseq-count -s no -m union -f bam’ and the gene annotation of Ensembl release 87 supplemented with ERCC, EGFP and mCherry features was used. The counts of all samples were merged into a single gene expression matrix, with genes in rows and samples in columns (Gene Expression Omnibus (GEO) accession number GSE141064, processed data).

To evaluate the cross-sample contamination, we sampled human and mouse cells alternatively. Reads were aligned to the mixed human: mouse reference genome (hg38 and mm10) using STAR and the number of reads per feature was counted by HTseq using the same settings as mentioned above. Digital gene expression matrices were generated for each species. We analysed the downstream data using R (v.3.5.0), plots generated using the R package ggplot2 (v.3.2.1).

Cell and gene filtering

ERCCs were used for technical evaluation by testing the correlation with the expected number of spike-in RNA molecules, and not included for further analysis. The top 20 genes found in negative controls (samples were merged) are probably due to misalignment/sequencing errors of the oligos (Supplementary Table 5), and were thus also removed from downstream analyses. Ribosomal protein-coding genes were removed as they confound the downstream differential gene expression analysis, consistent with previous findings⁵⁵. The downstream analysis followed the procedures of the Seurat R package (v.3.0). Samples showing low quality were filtered out, with quality cut-offs being: (1) the number of genes fewer than 1,000, (2) the mitochondrial read ratio more than >30% or (3) the uniquely mapped rate fewer than <30%. Then, the data were normalized to the total expression, multiplied by a scale factor of 10,000, after which a pseudo-count was added and the data were log transformed.

Feature selection, dimensionality, reduction clustering and others

The top 500 highly variable genes were chosen on the basis of the variance stabilizing transformation (vst) result (function, FindVariableFeatures(object, selection.method = ‘vst’)). Scaling was applied to all the genes using the function ‘ScaleData’. The scaled data of the 500 highly variable genes were used for principal component analysis (PCA) analysis. The first ten principal components were chosen for further clustering and t-SNE analysis on the basis of the ranking of the percentage of variance explained by each principal component (ElbowPlot function). Seurat-embedded graph-based clustering was applied. In both Live-seq and scRNA-seq data, postdifferentiated ASPCs were in early differentiation stage rendering their classification by treatment difficult to capture when in a dataset driven by more distinct cell types and states. Furthermore, the clustering of these cells was prone to be biased by batch as only one of the two batches contained both pre- and postdifferentiation ASPCs. Therefore, to capture ASPCs’ heterogeneity, ASPCs were clustered separately for both scRNA-seq and Live-seq data (default Seurat pipeline and data scaled for nCount and nGene), and the clustering found on the whole datasets, at resolution 0.2, was adapted. For data visualization, individual cells were projected on the basis of the ten principal component scores onto a two-dimensional map using t-SNE⁵⁶.

To evaluate the effect of varying sequencing depth on the clustering of Live-seq data, we down-sampled the raw counts to the desired number. The down-sampled matrices were analysed in the same way as described above. The clusters were consistent with the original analysis, indicating that the clusters were not driven by sequencing depth.

The DE genes were analysed as described below (Comparison of DE genes between Live-seq and scRNA-seq section). The pseudo-genes were filtered out⁵⁷. For the Gene Ontology analysis, the top first 100 DE genes, ordered by logFC, were loaded into the EnrichR package (v.3.0)³² to determine gene enrichment among the biological processes of Gene Ontology and predict the cell types using the Mouse Gene Atlas database. We carried out Gene Set Enrichment Analysis in R using the ClusterProfiler package⁵⁸ v.3.14.3 and msigdbr v.7.1.1 using default settings with all the gene expression changes. Both wikipathways-20200810-gmt-Mus_musculus.gmt and the Gene Ontology terms obtained from the R database org.Mm.eg.db (v.3.10.0) were used as reference.

The scores of cell cycle phases were calculated using the Seurat function ‘CellCycleScoring’ on the basis of canonical markers⁵⁹.

To integrate Live-seq and scRNA-seq data, the CCA and mutual nearest neighbours (MNNs)-based approaches embedded in Seurat v3 were used. Specifically, the top 500 highly variable genes of both datasets were chosen for PCA analysis independently. The first ten principal components of each were used to identify the anchors and for data integration. As another level of control, data stemming from cells belonging to each cluster in both Live-seq and scRNA-seq analyses were collapsed into a ‘bulk’ RNA-seq dataset after which the Pearson’s correlation between the bulk Live-seq and the bulk scRNA-seq datasets was determined. To evaluate the effect of varying sequencing depth on data integration, we down-sampled the raw counts to the desired number. The down-sampled matrices were then analysed in the same way as described above.

The cells sampled with Live-seq or scRNA-seq from the two cell types, RAW cells and ASPCs, were analysed on a cell type-by-cell-type basis following Seurat’s pipeline. For the ASPCs of Live-seq and scRNA-seq, the only batch containing both DMIR-treated and non-treated cells Live-seq the scRNA-seq were selected. The cells were filtered using the same filtering as previously described, then the data normalized and scaled for nCounts, nFeatures and batch (if more than one batch) were used to compute the PCA. Finally, the first ten principal components were used to compute the t-SNE with perplexity set to 10.

Comparison of DE genes between Live-seq and scRNA-seq

Differential gene expression analysis was conducted using edgeR⁶⁰ v.3.34.0. Only genes expressed in at least 5% of the data with a minimal count of 2 (filterByExpr() min.count = 2, min.prop = 0.05) and in at least 15% of the cells of one of the categories with a minimal count of 2 in the scRNA-seq data were considered. The dispersions and negative binomial glm were fitted (glmQLFit and glmQLFTest functions) for the model roughly 0 + categories + batch, the categories being the different groups to be tested. The quasi-likelihood F-tests were calculated for each group, with the null hypothesis being that there is no difference between the mean in the group of interest and the average over all the remaining categories (clusters). The P values were corrected using the Benjamini–Hochberg procedure. The genes with an FDR < 0.05, absolute logFC > 1, and expressed in at least 15% of one of the two groups were considered DE.

To compare scRNA-seq and Live-seq data, the differential expression analyses were performed separately for the two techniques. For each result, the log fold-expression changes of each gene were plotted against each other. The DE genes were highlighted (Bonferroni adjusted P value < 0.05 and an absolute logFC > 1). A linear model was fitted using lm() over (1) all the genes, (2) only the genes detected as DE in the scRNA-seq datasets or (3) the genes defined as DE in both the Live-seq and scRNA-seq datasets. Similar analyses were also applied per cell type basis, that is in RAW cells before and after LPS treatment and ASPC before and after differentiation.

Differential expression analysis per cell type was performed using edgeR package by contrasting between non-treated versus treated cells in the same fashion as described above. The genes with an FDR < 0.05, absolute logFC > 1 and expressed in at least 15% of one of the two groups of cells with at least two counts were considered DE.

For the two cell types mentioned above, cells of the scRNA-seq datasets were randomly selected to match the number of cells of the Live-seq data. The count matrices were then down-sampled per cell to have the same density distribution of the number of features as the corresponding Live-seq data. More precisely, the cells were ordered by the number of features in the scRNA-seq and Live-seq datasets and paired on the basis of this metric. The number of sampled reads (with replacement) of the scRNA-seq cells were defined so that the down-sampled data reached a similar number of features (absolute difference below 5) compared to its paired Live-seq cell. The down-sampled data were then (up-)sampled with replacement to match the library size of their paired Live-seq cell. Differential expression analysis was then performed as described above.

Live-seq and live-cell imaging integration

Among the 40 cells that were both subjected to Live-seq and tracked for LPS-induced Tnf-mCherry fluorescence, 17 of them passed the quality control as mentioned in the Live-seq section. For each cell in the time-course, we calculated the intercept (that is, basal expression) and slope (that is, extent of response) using a linear model between the time after LPS treatment and the natural log of the Tnf-mCherry fluorescence intensity. As we were mainly interested in the initial response and as the curve was linear from the first three to 7.5 h (Extended Data Fig. 9b), we only used values from this time window. To rank genes, we then constructed a linear model that predicts the intercept or slope of the mCherry response using the expression of a gene measured by Live-seq before treatment with LPS. Given the limited number of cells and thus to increase the statistical power of the models, only the 500 most variable genes from both the RAW-G9 Live-seq and scRNA-seq data were used, but removing genes with a dispersion lower than 0.1. The genes were then ranked on the basis of their respective R² values. We used two tests to assess the significance of each gene: an F-test on the overall model as implemented using R’s lm function, and a bootstrapping approach in which we randomly sampled cells with replacement to calculate an empirical P value. P values were corrected for multiple testing using the Benjamini–Hochberg procedure as implemented using the R’s p.adjust function.

Trajectory inference

To infer the trajectory from conventional scRNA-seq data, the dynverse R package (v.0.1.2)¹⁰ with multiple wrapped trajectory inference methods was used^{61,62,63,64,65}. The parameters shown below were used for the ‘answer_questions’ function- ‘multiple_disconnected = NULL, expect_topology = NULL, expected_topology = NULL, n_cells = 554, n_features = 1,000, memory = ‘100GB’, docker = TRUE’. The most suggested methods were chosen on the basis of the guidelines provided by dynverse. In addition, we also applied the Monocle DDRTree method⁶⁶ given its widespread use. All these methods were run in the docker with default parameters. For the RNA velocity analysis, annotated spliced, unspliced and spanning reads in the measured cells were generated in a single loom file using the command line ‘velocyto run_smartseq2 -d 1’ function. This also generates an HDF5 file containing detailed molecular mapping information that was used for the analysis model on the basis of gene structure. Genes were filtered to have a min.max.cluster.average of at least 5, 1 and 0.5 for the exonic, intronic and spanning read expression matrices, respectively. Three different algorithms were used following the velocyto (v.0.6) pipeline (http://pklab.med.harvard.edu/velocyto/notebooks/R/chromaffin2.nb.html): (1) cell kNN pooling with the gamma fit based on extreme quantiles, with function ‘gene.relative.velocity.estimates(deltaT = 1, kCells = 5, fit.quantile = 0.05)’, (2) relative gamma fit without cell kNN smoothing, with function ‘gene.relative.velocity.estimates(deltaT = 1, deltaT2 = 1, kCells = 1, fit.quantile = fit.quantile)’ and (3) velocity estimate based on gene structure. Here, the unfiltered intronic and spanning expression matrix was used to include more genes for the genome-wide model fit.

Analysis of the most variable genes from primary macrophage scRNA-seq data

We retrieved single-cell expression data of five macrophage subsets from the GEO database (GSE117081)⁶⁷. We calculated for each gene a standardized variance by modelling the relationship between the observed mean expression and variance using local polynomial regression, as implemented in the Seurat (v.3.1.4) FindVariableFeatures function. Genes of the KEGG NF-κB signalling pathway downstream TLR4 receptor were highlighted.

Bioethics

All mouse experiments were conducted in strict accordance with the Swiss law, and all experiments were approved by the ethics commission of the state veterinary office (licence number VD 3406, valid from 14 October 2018 to 14 January 2022). Mice were housed under specific pathogen-free conditions at 20–24 °C with 45–65% humidity and a 12 light/12 dark cycle.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.