High-throughput total RNA sequencing in single cells using VASA-seq

Ethics statement

Experiments were performed in accordance with European Union guidelines for the care and use of laboratory animals and under the authority of appropriate United Kingdom governmental legislation. Use of animals in this project was approved by the Animal Welfare and Ethical Review Body for the University of Cambridge, and relevant Home Office license PPL (7677788) is in place.

Cell lines

HEK293T cells were passaged every second day and cultured in T75 flasks. The culture media was DMEM-F12 (Thermo Fisher Scientific) supplemented with 10% heat-inactivated FBS (Thermo Fisher Scientific) and 100 U ml⁻¹ of penicillin–streptomycin (Thermo Fisher Scientific). For passaging, the cells were washed with 10 ml of ice-cold 1× PBS (Lonza) twice. Then, 9 ml of PBS was added to the flask, and cells were detached by adding 1 ml of 10× trypsin-EDTA (Sigma-Aldrich) and incubated at 37 °C for 5 minutes. Trypsin-EDTA was then inactivated with 15 ml of DMEM-F12 + 10% FBS and incubated at 37 °C for 5 minutes. The cells were then pelleted at 300g for 3 minutes, and the supernatant was aspirated. After aspiration of the supernatant, the cells were washed twice in PBS and viability-assessed and counted before encapsulation.

mESCs were passaged every other day and cultured in 2i+LIF medium. In brief, DMEM/F-12 nutrient mixture without L‑glutamine (Thermo Fisher Scientific) and neurobasal medium without L‑glutamine (Thermo Fisher Scientific) in a 1:1 ratio, 0.1% sodium bicarbonate (Thermo Fisher Scientific), 0.11% bovine albumin fraction V solution (Thermo Fisher Scientific), 0.5× B-27 supplement (Thermo Fisher Scientific), 1× N-2 supplement (Cambridge Stem Cell Institute, made in-house), 50 µM 2-mercaptoethanol (Thermo Fisher Scientific), 2 mM L-glutamine (Thermo Fisher Scientific), 100 U ml⁻¹ of penicillin–streptomycin (Thermo Fisher Scientific), 12.5 µg ml⁻¹ of insulin zinc (Thermo Fisher Scientific), 0.2 µg ml⁻¹ of mLIF (Cambridge Stem Cell Institute), 3 µM CHIR99021 (Cambridge Stem Cell Institute) and 1 µM PD0325901 (Cambridge Stem Cell Institute). Culture dishes were coated with 0.1% gelatine in PBS for at least 30 minutes. Cells were detached with 500 μl per six-well of Accutase (Merck) for 3 minutes at 37 °C. The detached cells were transferred into 9.5 ml of washing medium (DMEM/F-12 with 1% bovine albumin fraction V solution) and centrifuged at 300g for 3 minutes. The supernatant was aspirated, and the cell pellet was resuspended in 2i+LIF medium and re-plated at 80,000 cells per six-well. For the encapsulation process, the cells were washed twice in PBS, viability-assessed and counted before dilution to the correct concentration.

Murine embryo collection and dissociation

Pregnant C57BL/6 female mice were purchased from Charles River Laboratories or obtained from natural mating of C57BL/6 mice in-house. Mice were maintained on a lighting regimen of 12-hour light/dark cycle with food and water supplied ad libitum. Detection of a copulation plug after natural mating indicated E0.5. After euthanasia of the females using cervical dislocation, the uteri were collected into PBS (Lonza) with 2% heat-inactivated FBS (Gibco, Thermo Fisher Scientific), and the embryos were immediately dissected and processed for scRNA-seq. Mouse embryos were dissected at timepoints E6.5, E7.5, E8.5 and E9.5, as previously reported⁵⁴. Embryos from the same stage were pooled in a LoBind tube (Eppendorf). E8.5 and E9.5 embryos were cut into pieces under stereomicroscopy before collecting into a tube. The pooled samples were centrifuged at 300g for 5 minutes. The supernatant was aspirated, and 100–200 µl of TrypLE Express (Gibco) dissociation reagent was added to the samples. The tube was incubated at 37 °C for a minimum of 7 minutes (or until completely dissociated) in an orbital shaker. Subsequently, 1 ml of FBS was added to the tube to inactivate TrypLE. The sample was repeatedly centrifuged and washed with PBS before finally being resuspended in PBS supplemented with 0.04% BSA and filtered through a 40-µm Flowmi Tip Strainer (Thermo Fisher Scientific).

VASA-plate: cell sorting in 384-well plates

Single cells were sorted into 384-well hardshell plates (BioRad) using a BD FACSJazz. Each well was pre-filled with 5 µl of mineral oil (Sigma-Aldrich) and 50 nl of CEL-seq2/SORT-seq^1,27 primer with a concentration of 0.25 µM. Plates were sealed (Greiner, SILVERseal sealer, 676090) and spun down at 2,000 revolution centrifugal force (r.c.f). for 1 minute (Eppendorf 5810R) before being stored at −80 °C.

VASA-plate: cell lysis and RNA fragmentation

All dispensions were carried out with a NanoDrop II (Innovadyne Technologies), all incubations with a GeneAmp PCR System 9700 Thermal Cycler (Applied Biosystems) and all spinning steps with an Eppendorf 5810R, unless otherwise specified. Next, 50 nl of lysis and fragmentation mix (3.4× First-Strand Buffer (Invitrogen), 1.2 mU of Thermolabile Proteinase K (New England Biolabs (NEB))) and 0.2% IGEPAL CA-630 (Sigma-Aldrich) were added to each well. Plates were sealed and spun down at 2,000 r.c.f. for 2 minutes. Lysis was carried out at 25 °C for 1 hour, followed by 55 °C at 10 minutes. Plates were snap-chilled on ice before fragmentation was carried out at 85 °C for 3 minutes. Plates were snap-chilled, spun down at 2,000 r.c.f. for 1 minute and stored on ice before next dispensation.

VASA-plate: RNA repair and poly(A) tailing

Next, 50 nl of RNA repair and poly(A)-tailing mix (0.6× First-Strand Buffer (Invitrogen), 20 mM DTT (Invitrogen), 7.5 nM ATP (NEB), 37.5 mU of E. coli Poly(A) Polymerase (NEB), 50 mU of T4 PNK (NEB) and 10 mM RNaseOUT (Invitrogen)) were added to each well. Plates were sealed and spun down at 2,000 r.c.f. for 2 minutes. Repair and tailing were carried out at 37 °C for 1 hour. Plates were snap-chilled, spun down at 2,000 r.c.f. for 1 minute and stored on ice before next dispensation.

VASA-plate: reverse transcription

Next, 50 nl of reverse transcription mix (2 mM (each) dNTP mix (Promega) and 0.8 U of SuperScript III (Invitrogen)) was added to each well. Plates were sealed and spun down at 2,000 r.c.f. for 2 minutes. Reverse transcription was carried out at 50 °C for 1 hour. Plates were snap-chilled, spun down at 2,000 r.c.f. for 1 minute and stored on ice before next dispensation.

VASA-plate: second-strand synthesis

Next, 1,100 nl of second-strand synthesis mix (1.14× Second-Strand Buffer (Invitrogen), 0.23 mM (each) dNTP mix (Promega), 0.35 U of E. coli DNA Polymerase I (Invitrogen) and 20 mU of RNaseH (Invitrogen)) was added to each well. Plates were sealed and spun down at 2,000 r.c.f. for 2 minutes. Second-strand synthesis was carried out at 16 °C for 2 hours, followed by 85 °C for 20 minutes. Plates were snap-chilled, spun down at 2,000 r.c.f. for 1 minute and stored on ice before pooling. The protocol for pooling and in vitro transcription (IVT) was the same as SORT-seq²⁷.

VASA-drop: design of the droplet generation device

The droplet generation device for compressible barcoded bead and single-cell co-encapsulation (Extended Data Fig. 1c) was modified from previous designs^5,59. The flow-focusing junction (80 µm) was narrowed to generate smaller droplets (0.55 nl) at high throughput (115 Hz).

VASA-drop: design of droplet picoinjection devices

The design of both droplet picoinjector devices is based on the findings of a previous study²⁸. Several key features were added to the architecture of previous designs to ameliorate the robustness of the injections in large droplets containing compressible barcoded beads:

1. 1.

   Emulsion-diluting oil inlet, number 2 (Extended Data Fig. 1d), which reduces the packing of the emulsion to eliminate fragmentation of densely packed droplets before being re-injected in the picoinjection channel. This design feature allows for packed droplets to arrange into an evenly spaced monolayer that reduces fluctuations in volume of droplets after picoinjection.

2. 2.

   Smooth narrowing of the reinjection chamber facilitating the ordering of droplets before spacing, which reduced droplet break-up.

3. 3.

   Deepening of the outlet junction before the outlet, number 5 (Extended Data Fig. 1d,e, deep blue color), which stabilizes droplets and reduces droplet merging, which was observed during the rapid transition from the shallow microfluidic channel to a wide tubing or collection tip.

VASA-drop: photolithography of microfluidic molds

The channel layout for the microfluidic chips was designed using AutoCAD (Autodesk) and printed out on a high-resolution film photomask (Micro Lithography Services). The designs in Extended Data Fig. 1 are deposited on https://openwetware.org/wiki/DropBase:Devices and can be found in the supplementary file ‘SI_VASAdrop CAD designs_5masks.dxf’. The microfluidic devices were fabricated following standard hard and soft lithography protocols that can be performed in local cleanrooms or outsourced to contract manufacturing companies. First, microfluidic molds were patterned on a 3-inch silicon wafer (MicroChemicals) using high-resolution film masks (Micro Lithography Services) and SU-8 2075 photoresists (Kayaku Advanced Materials). An MJB4 mask aligner (SÜSS MicroTec) was used to UV expose all the SU-8 spin-coated wafers. The thickness of the structures (corresponding to the depth of channels in the final microfluidic devices) was measured using a DektakXT Stylus profilometer (Bruker).

We used the following settings for photolithography:

VASA-drop: soft lithography

To manufacture PDMS microfluidic devices, 20–30 g of silicone elastomer base and curing agent (Sylgard 184, Dow Corning) were mixed at a 10:1 (w/w) ratio in a plastic cup and de-gassed in a vacuum chamber for 30 minutes. PDMS was then poured on a master wafer with SU-8 structures and cured in the oven at 65 °C for at least 4 hours. Next, the inlet holes were punched using two types of biopsy punchers with plungers (Kai Medical Laboratory): a 1.5-mm-diameter punch was used to make the inlet for the cell delivery tip, number 2 (Extended Data Fig. 1c); outlet for droplet collection tip, inlet number 5 (Extended Data Fig. 1c,d); and the inlets for droplet reinjection, number 1 (Extended Data Fig. 1d,e); whereas other inlets were made using a 1-mm-wide biopsy puncher. The patterned PDMS chip was then plasma bonded to a 52 mm × 76 mm × 1 mm (length × width × thickness) glass slide (VWR) in a low-pressure oxygen plasma generator (Femto, Diener Electronics). Next, the hydrophobic modification of microfluidic channels was performed by flushing the device with 1% (v/v) trichloro(1H,1H,2H,2H-perfluorooctyl)silane (Sigma-Aldrich) in HFE-7500 (3M) and baked on a hot plate at 75 °C for at least 30 minutes to evaporate the fluorocarbon oil and silane mix.

Although we have not used these commercial suppliers, we propose the following list of contract manufacturers for users who may not have access to photolithography/soft lithography: Fivephoton Biochemicals, Darwin Microfluidics, uFluidix, Flowjem and Microfactory.

VASA-drop: cell loading and droplet collection/re-injection chamber manufacturing

Cell injection chamber

The cells were loaded in a cell loading tip pre-filled with mineral oil (Sigma-Aldrich). To manufacture the cell loading tip, a low-retention pipette tip with 200-µl volume capacity (Axygen) was cut at the top, in parallel to the rim and under the filter. A solidified 3-mm-thick piece of PDMS (Dow Corning) was punched from a slab of PDMS with a 5.0-mm sampling tool (EMS-Core). The circular piece of PDMS was then biopsy-punched with a 1-mm-wide biopsy puncher (Kai Medical Laboratory) in the middle. The circular piece of PDMS was pushed inside the tip while remaining parallel to the upper rim of the tip. A 1-ml glass syringe (SGE) was then pre-filled with 1 ml of mineral oil and connected to a 30-cm-long tubing (Portex, Smiths Medical) that can be inserted to a hole in the middle of the circular PDMS piece in the tip. Next, the tip was pre-filled with mineral oil by manually pushing the syringe, and the cell-containing solution was further aspirated with care as to not introduce any air bubble in the system. The tip can then be connected to the cell-encapsulation PDMS device, inlet number 2 (Extended Data Fig. 1c,g), and injection rates are modulated by a Nemesys syringe pump (Cetoni).

Droplet collection and reinjection chamber

A second type of tip chamber was designed to collect, incubate and re-inject droplets for each microfluidic step. To this end, a 5-mm-thick PDMS piece was punched from a slab of PDMS with an 8.0-mm sampling tool (EMS-Core) and re-punched in its center using a 1-mm-wide biopsy puncher (Kai Medical Laboratory), and a 30-cm-long tubing (Portex, Smiths Medical) was connected to the latter punched hole. The resulting piece of PDMS was then inserted into a 1-ml filterless pipette tip (Sigma-Aldrich) with a parallel orientation to the rim. Unsolidified PDMS (Dow Corning, 1:10 (w/w) ratio, de-gassed) was then deposited into the space between the rim and the circular PDMS piece at the top. The tip was then incubated at 65 °C for at least 4 hours and connected to a 1-ml glass syringe (SGE) pre-filled with mineral oil. The tip was then pre-filled with mineral oil by manually pushing the connected syringe. To collect the droplets after the initial encapsulation or at the end of the first picoinjection, the tip can be connected to the outlet of the devices, inlet number 5 (Extended Data Fig. 1c,d), and the syringe is disconnected to allow the evacuation of mineral oil as the tip gets loaded. For each of the two droplet picoinjection steps, the mineral oil can be pushed using a Nemesys syringe pump (Cetoni) to re-inject droplets into the picoinjectors, inlet number 1 (Extended Data Fig. 1d,e). For each of the re-injection and collection steps, the PDMS-punched holes on the microfluidic device need to be primed with 5% (w/w) 008-FluoroSurfactant (RAN Biotechnologies) in HFE-7500 (3M) to avoid a trapped air bubble to perturbate the stability of re-injection or the integrity of the emulsions. After droplet collection during the encapsulation and first picoinjection, the tip can be closed by inserting the narrower end of the tip into a glass-bonded PDMS plug, which closes the system and allows for incubation of the tip in the water bath (Extended Data Fig. 1f). The glass-bonded PDMS plug was fabricated before the experiment by punching a 8-mm-thick piece of PDMS with a 1.5-mm biopsy puncher that was then bonded to microscopy glass using an oxygen plasma.

VASA-drop: microfluidic device operation

Polyacrylamide beads manufacturing

Barcoded polyacrylamide beads were manufactured following a previously described protocol⁵⁹. In brief, a polyacrylamide mix was used to generate 60 µm of water-in-oil emulsions using a single-inlet flow-focusing device and collected in a 1.5-ml LoBind tube (Eppendorf) containing 200 µl of mineral oil (Sigma-Aldrich). The droplets were solidified overnight at 65 °C, de-emulsified using a 20% 1H,1H,2H,2H-perfluoro-1-octanol (Alfa Aesar) in HFE-7500 (3M) solution and stored at 4 °C for up to 6 months.

Co-encapsulation of cells and barcoded beads

A detailed protocol⁵⁹ was used as a reference for droplet generation. First, the microfluidic droplet generation chip was installed on the stage of an inverted microscope (Olympus XI73). Next, two pieces of polyethylene tubing (Portex, Smiths Medical) were connected to two 1-ml gas-tight syringes (SGE) and filled with PBS (Lonza). The tubing was manually filled with PBS, and a small, 1-cm-long air bubble was left at the end tip of each tubing. The bead suspension and lysis mix were manually aspirated to the tubing, and the small air bubble provided a separation between the reagents and the PBS buffer. Then, 150 µl of cell suspension was manually aspirated into the cell loading tip pre-filled with mineral oil (Sigma-Aldrich). A fourth 2.5-ml glass syringe (SGE) was filled with 5% (w/w) 008-FluoroSurfactant (RAN Biotechnologies) in HFE-7500 (3M). Next, all three tubings and the cell chamber with cell suspension were inserted to the corresponding inlets of the droplet generation chip (Extended Data Fig. 1c). Four Nemesys syringe pumps (Cetoni) were used to flow each component, and the droplet formation was monitored using ×4 or ×10 objectives (Olympus) and a fast camera (Phantom Miro eX4) connected to the inverted microscope. After the device was primed and droplet generation was stabilized, the collection chamber was connected to the outlet.

Microfluidic device operation—picoinjection

Before starting the picoinjection of droplets containing single-cell lysates, the electrode sections, numbers 6 and 7 (Extended Data Fig. 1d,e) of the devices, were pre-filled with filtered 5 M NaCl as previously described⁶⁰. The picoinjection chip was filled with 5% (w/w) 008-FluoroSurfactant (RAN Biotechnologies) in HFE-7500 (3M) using a pre-filled 2.5-ml glass syringe (SGE) connected to a piece of tubing (Portex, Smiths Medical). The reaction mix was primed, and the tip containing the emulsions (with fluorinated oil evacuated by pushing the glass syringe until the emulsions reached the exit of the tip) was primed and connected to the device. Next, flows of droplet emulsions, the reaction mix, the emulsion-diluting oil, number 2 (Extended Data Fig. 1d,e), and the droplet-spacing oil, number 3 (Extended Data Fig. 1d,e) were applied using the Nemesys syringe pumps (Cetoni). The droplets were diluted in a first instance in the re-injection chamber and then spaced with the second stream of oil in a flow-focusing re-injection junction. 5% (w/w) 008-FluoroSurfactant (RAN Biotechnologies) in HFE-7500 (3 M) was used for both diluting and spacing of droplets. The function generator (AIM & Thurlby Thandar Instruments) was set to generate square waves of 2.5V amplitude and 10kHz frequency, which were further amplified 100 times to 250 V by a Trek 601C-1 amplifier, which enabled coalescence-activated injection of the reagent into the droplets. The droplets were collected in a 1-ml collection tip connected at the outlet, number 5 (Extended Data Fig. 1d,e).

VASA-drop: polyacrylamide bead barcoding

The bead barcoding procedure was performed as previously described⁵⁹ with the inDrop v3 barcoding scheme⁶¹. In brief, the solidified barcoded beads were filtered and dispensed in four 96-well plates containing the first barcode from the inDrop v3 design, and the bead-bound adapter was extended using a Bst 2.0 DNA polymerase (NEB) after annealing the barcoded oligonucleotides. The reaction was then stopped, and the second strand was removed using a sodium hydroxide treatment. The second barcode was added in a similar fashion, and the beads were stored for up to 6 months at 4 °C.

VASA-drop: cell encapsulation in water-in-oil emulsions

For the cultured cells and the embryos, we used a loading concentration of 450 cells per µl in 1× PBS (Lonza) with 15% OptiPrep (Sigma-Aldrich). The lysis mix was made fresh before each encapsulation, as follows: 0.5 mM dNTPs each (Thermo Fisher Scientific), 0.52% IGEPAL-CA630 (Sigma-Aldrich), 40 mM UltraPure Tris-HCl pH 8 (Life Technologies), 3.76× First-Strand Buffer (Invitrogen), 3 mM magnesium chloride (Ambion) and 6 U ml⁻¹ of Thermolabile Proteinase K (NEB). The barcoded PAAm beads were prepared for encapsulation as previously described⁵. The lysis mix and bead suspensions were loaded in the tubing of two individual 1-ml SGE glass syringes filled with PBS (Lonza) and separated by an air bubble from the reagents in the tubing. The cells were loaded into a cell injection container pre-filled with mineral oil (Sigma-Aldrich). The injection flow rates for the droplet encapsulation device (Extended Data Fig. 1c) were as follows: the cell suspension was flown at 85 µl per hour, number 2 (Extended Data Fig. 1c); the bead suspension was flown at 65 µl per hour, number 3 (Extended Data Fig. 1c); the lysis solution was flown at 75 µl per hour, number 1 (Extended Data Fig. 1c); and 5% (w/w) 008-FluoroSurfactant (RAN Biotechnologies) in HFE-7500 (3M) was flown at 450 µl per hour using a 2.5-ml glass syringe (SGE), number 4 (Extended Data Fig. 1c). All flow rates for each microfluidic manipulation were controlled using Nemesys pumps (Cetoni). The average droplet size was ~0.55 nl for these flow rates and a microfluidic device depth of 80 µm. The droplets were collected for approximately 1 hour in a 1-ml pipette tip (Greiner) pre-filled with mineral oil at the outlet, number 5 (Extended Data Fig. 1c), and connected to a tubing via a PDMS connector (Extended Data Fig. 1g). The collection tip was then closed by connecting a 1-ml SGE glass syringe pre-filled with mineral oil to the tubing and the tip was then connected to a glass-bonded PDMS plug (Extended Data Fig. 1f).

VASA-drop: cell lysis and RNA fragmentation

The tip container was further left at room temperature (23 °C) for 20 minutes to allow for cell lysis to occur, and the tip was then placed under a High-Intensity UV Inspection Lamp (UVP) that was switched on for 7 minutes for barcode photocleavage (Extended Data Fig. 1h). The container was then submerged in a water bath (Grant JB) placed at 85 °C for 6 minutes and 30 seconds. After incubation, the container was immediately submerged in an ice bucket filled up with half proportions of ice and water.

VASA-drop: first picoinjection for RNA repair and poly(A) tailing

The droplets were re-injected in the first picoinjector device with the shorter re-injection channel (Extended Data Fig. 1d) to perform coalescence-induced merging with a poly(A) solution consisting of 26.6 mM Tris-HCl, pH 8 (Invitrogen), 15.8 mM DTT (Invitrogen), 0.83× First-Strand Buffer (Invitrogen), 0.19 mM ATP (NEB), 3.15 kU ml⁻¹ of T4 Polynucleotide Kinase (NEB), 250 U ml⁻¹ of E. coli poly(A) polymerase and 2.6 kU ml⁻¹ of RNaseOUT (Applied Biosystems). The merging was applied by pre-filling the electrode section, numbers 6 and 7 (Extended Data Fig. 1d), of the device with 5 M NaCl, as previously described⁶⁰. The flow rates used were 200 µl per hour for the droplet emulsion, number 1 (Extended Data Fig. 1d); 60 µl per hour for the poly(A) mix, number 4 (Extended Data Fig. 1d); 50 µl per hour for the emulsion-diluting oil, number 2 (Extended Data Fig. 1c); and 400 µl per hour for the droplet-spacing oil, number 3 (Extended Data Fig. 1d). This generated ~0.8 nl of droplets at 70 Hz. The droplets were collected in a 1-ml collection tip (Greiner) pre-filled with mineral oil and inserted to the outlet, number 5 (Extended Data Fig. 1d). At the end of the picoinjection, the collection tip was closed by connecting a 1-ml glass syringe (SGE) pre-filled with mineral oil (Sigma-Aldrich) to the tubing and connecting the narrower end of the tip to the glass-bonded PDMS plug. The tip container was then incubated for 25 minutes at room temperature (23 °C) and 8 minutes at 37 °C in a water bath (Grant JB) and then submerged in an ice-cold water bath for 2 minutes. The droplets were then processed for the second picoinjection.

VASA-drop: second picoinjection for reverse transcription

The droplets were re-injected in the second picoinjector (Extended Data Fig. 1e) similarly to the previous step, although this time the droplets were collected in fractions of ~1,000 cells (~27 µl of loaded droplets) in 1-ml LoBind tubes (Eppendorf) pre-filled with 200 µl of mineral oil. The droplets were injected with a reverse transcription mix constituted of 25 mM Tris-HCl, pH 8 (Invitrogen), 8 mM DTT (Invitrogen), 0.75× First-Strand Buffer (Invitrogen), 1 mM dNTPs, 20 kU ml⁻¹ of SuperScript III (Invitrogen) and 1.2 kU ml⁻¹ of RNAseOUT (Applied Biosystems). The flow rates for the second picoinjection were as follows: 70 µl per hour for the emulsion-diluting oil, number 2 (Extended Data Fig. 1e); 700 µl per hour for the droplet-spacing oil, number 3 (Extended Data Fig. 1e); 300 µl per hour for the re-injected droplets, number 1 (Extended Data Fig. 1e); and 255 µl per hour for the reverse transcription mix, number 4 (Extended Data Fig. 1e). The collected fractions were incubated at 50 °C for 2 hours and then heat-inactivated at 70 °C for 20 minutes. For de-emulsification of the droplets, the mineral oil and the excessive fluorocarbon oil phase were aspirated and discarded. Then, 200 µl of filtered HFE-7500 was added to the emulsions, followed by 200 µl of 100% 1H,1H,2H,2H-perfluoro-1-octanol. The tubes were centrifuged for 5 seconds on a tabletop centrifuge, and then 300 µl of the oil phase was removed and 100 µl of fresh HFE-7500 oil was added, as well as 50 µl of TE buffer (Zymo). At this point, the fractions were stored at −80 °C. The protocol, up to and including the IVT step, was the same as for inDrop⁵⁹.

VASA-plate and VASA-drop: downstream library preparation and sequencing

For VASA-plate: after IVT, 2 µl of ExoSAP-IT (Applied Biosystems) was added, and each sample was incubated at 37 °C for 15 minutes. For both VASA-plate and VASA-drop: a 1.8× volumetric ratio AMPure XP clean-up was then performed, and the amplified RNA (aRNA) was eluted in 10 µl of nuclease-free water. The purified aRNA concentration was measured using a Qubit (Invitrogen), and the concentration was adjusted to a maximum of 100 ng µl⁻¹. Next, 6 µl per sample was mixed with 2 µl of rRNA depletion probes (25 µM) (reverse complement of published probes⁶²) and 2 µl of hybridization buffer (pH 7.5, 500 mM Tris-HCl, 1 M NaCl). Samples were incubated at 95 °C for 2 minutes and brought to 45 °C with a gradient of 0.1 °C per second. Once the probes were hybridized, 2 µl of Thermostable RNAseH (Epicentre) and 8 µl of RNAseH buffer (pH 7.5, 125 mM Tris-HCl, 250 mM NaCl, 50 mM MgCl₂) was added. The reaction was incubated at 45 °C for 30 minutes and further kept on ice. Next, 4 µl of RQ DNAse I (Promega), 21 µl of nuclease-free water and 5 µl of CaCl₂ (10 mM) were added to the reaction mixture. The mixture was further incubated at 37 °C for 30 minutes, followed by snap-cooling on ice. A 1.6× volumetric ratio AMPure XP clean-up was then performed, and the aRNA was eluted in 6 µl of nuclease-free water. Next, 1 µl of RA3 ligation oligonucleotide (20 µM; Supplementary Table 12) was added to 5 µl of the aRNA, and the reaction was brought to 70 °C for 2 minutes, followed by snap-cooling on ice. This was followed by the addition of 1 µl of 10× T4 RNA ligase reaction buffer (NEB), 1 µl of NEB T4 RNA Ligase2, truncated (NEB), 1 µl of RNAseOUT (Invitrogen) and 1 µl of nuclease-free water, The reaction was incubated at 25 °C for 1 hour, followed by snap-cooling on ice. The adapter-ligated aRNA was then mixed with 1 µl of dNTPs (10 mM each) (Promega) and 2 µl of RTP oligonucleotide (20 µM; Supplementary Table 12). The mixture was incubated at 65 °C for 5 minutes, followed by snap-cooling on ice. Next, 4 µl of 5× First-Strand Synthesis Buffer (Invitrogen), 1 µl of nuclease-free water, 1 µl of 0.1 M DTT (Invitrogen), 1 µl of RNAseOUT and 1 µl of SuperScript III were added to the sample. The reaction was incubated at 50 °C for 1 hour, followed by 70 °C for 15 minutes and then snap-cooled on ice. To reduce excess RNA material, 1 µl of RNAseA (Thermo Fisher Scientific) was further added to each tube, and the cDNA was incubated at 37 °C for 30 minutes, followed by a 1× volumetric AMPure XP clean-up. The cDNA was eluted in 20 µl of nuclease-free water. Half the material was used for the final PCR (10 µl). Each sample was mixed with 25 µl of NEBNext High-Fidelity 2× PCR Master Mix (VASA-plate) or Kapa HiFi HotStart PCR Mix (VASA-drop), 4 µl of PE1/PE2 primer mix (5 μM each)^1,27 (VASA-plate) or 5 µl PE1/PE2 primer mix (5 μM each) (Supplementary Table 12) (VASA-drop) and 11 µl (VASA-plate) or 10 µl (VASA-drop) of nuclease-free water. The samples were amplified with the following PCR programs. VASA-plate: initial heat denaturation for 30 seconds at 98 °C, 7–8 cycles for 10 seconds at 98 °C, 30 seconds at 60 °C, 30 seconds at 72 °C and final extension for 10 minutes at 72 °C. VASA-drop: initial heat denaturation for 2 minutes at 98 °C, two cycles for 20 seconds at 98 °C, 30 seconds at 55 °C, 40 seconds at 72 °C, 5–6 cycles for 20 seconds at 98 °C, 30 seconds at 65 °C, 40 seconds at 72 °C and final extension for 5 minutes at 72 °C. Each amplified and indexed sample was purified twice using a 0.8× volumetric ratio of AMPure XP beads and eluted in 10 µl. Final libraries were checked for proper length on a Bioanalyzer (Agilent), and concentration was measured with a Qubit (Invitrogen). A detailed catalog of reagents and instrumentation is provided in Supplementary Table 14.

The VASA-drop samples were sequenced on a NovaSeq 6000 S2, 300 cycles flow cell (Illumina), with the following parameters: Read1 247 cycles, Index1 31 cycles, Index2 8 cycles, Read2 14 cycles. VASA-plate samples were sequenced on a NextSeq 500, high-output 150 cycles flow cell (Illumina), with the following parameters: Read1 26 cycles, Index 8 cycles, Read2 135 cycles.

FASTQ file pre-processing in VASA-drop and 10x Chromium

Raw reads for VASA-drop were pre-processed with a Python script to have a favorable format for the pipeline (four reads were demultiplexed and rearranged into two reads). For each Read1, the UMI (6 nucleotides (nt) long in VASA-seq, 10 nt long in 10x Chromium) and the cell-specific barcode (16-nt long in VASA-seq, 14-nt long in 10x Chromium) were extracted. To determine the number of cells in each sample, first the total number of raw reads was determined for each possible barcode. Next, we plotted the histogram of log₁₀(read number) for each possible barcode, which we fitted to a polynomial function that shows two or three minima. We used the position of the minimum with the highest value of log₁₀(reads) as the threshold: only barcodes with reads above this threshold were used for downstream analysis. We merged sequenced barcodes that can be uniquely assigned to an accepted barcode with a Hamming distance of 2 nt or less.

FASTQ file pre-processing in VASA-plate

Read1 starts with a 6-nt-long UFI/UMI, followed by an 8-nt-long cell-specific barcode. There are only 384 cell-specific barcodes, each one corresponding to a well in a 384-well plate (available in GSE176588). We merged sequenced barcodes that can be uniquely assigned to an accepted barcode with a Hamming distance of 1 nt or less.

Mapping data (VASA-seq, 10x Chromium and Smart-seq v3)

Read2 was assigned to accepted barcodes (extracted from Read1) and trimmed with TrimGalore (version 0.4.3) with default parameters. Next, homopolymers at the end of the read were removed with cutadapt (version 2.10)⁶³.

In silico ribosomal depletion was performed by mapping the trimmed reads to mouse or human rRNA (National Center for Biotechnology Information) using bwa mem and bwa aln (version 0.7.10)⁶⁴. Multi-mappers and single-mappers were filtered out. The remaining reads were mapped to the mouse GRCm38 genome (Ensembl 99) or to the human GRCh38 genome (Ensembl 99) using STAR⁶⁵ with default parameters. Assignment of reads to gene biotypes was performed according to the following hierarchy:

• All mappings falling in TEC transcripts were discarded.

• Reads fully falling inside a region annotated as miscRNA, mtRNA, mttRNA, TrJGene, miRNA, rRNA, ribozymes, sRNA, scaRNA, snRNA or snoRNA (for example, biotypes that do not have annotated introns) were assigned to such regions.

• When a read maps to multiple genes simultaneously (because of annotation overlap in the reference GTF file), exonic annotations were given preference to introns. In case all references are exonic or intronic, the read is assigned to a gene whose name is the sequence of all the target gene names.

• Reads falling into exon–intron junctions or inside introns are assigned to unspliced transcripts. Reads falling inside exonic regions are assigned to spliced transcripts.

• If at least one UFI of the same cell from the same transcript has been assigned to an unspliced transcript (because it is mapped in an intron or an intron–exon junction), all the other reads with the same UFI of the same cell for the same transcript are automatically assigned to unspliced transcripts even if they mapped to exons exclusively.

Benchmarking against other methods

To determine the number of potential doublets, barcodes with more that 75% of the genes assigned to only one of either mouse or human were considered singlets. Cells with fewer than 7,500 UFIs were filtered out and not assigned to any organism. For gene body coverage, the BAM files for all single cells were used as a bulk. QoRTs⁶⁶ was used to calculate coverage, and only protein-coding genes were kept. For Smart-seq3, both reads containing a UMI (5′ reads) and non-UMI-containing reads were used together. Average coverages were used for the plotting. To determine percentages of different biotypes, all single cells were used as a bulk. UMI/UFI filtering was carried out for reads where this was possible. For Smart-seq3, both reads containing a UMI (5′ reads) and non-UMI-containing reads were used together. For the gene detection assay, only cells that had been sequenced to the highest numbers of reads (reads with proper barcode and quality/homopolymers trimming) were used (75,000 for saturation curve and 750,000 deep sequencing comparison) (Extended Data Fig. 2f). For Smart-seq3, four cells, with much lower reads than the rest, were removed as they were considered failed libraries. Downsampling was carried out with DropletUtils⁶⁷ on the count matrices (non-UMI/UFI filtered), based on the number of input reads and target reads, and only uniquely assigned genes were counted. For the percentage of intronic reads, each cell was used individually. UMI/UFI filtering was carried out for reads where this was possible. For Smart-seq3, both reads containing a UMI (5′ reads) and non-UMI-containing reads were used together. Mean and standard deviation were calculated and plotted.

scRNA seq analysis for mouse VASA-seq libraries and individual timepoints

The Scrublet⁶⁸ and Scanpy⁶⁹ packages were used together with custom-made code. In brief, for VASA-seq, only cells with more than 10⁴ (E6.5), 10^3.5 (E7.5, E8.5) and 10³ (E9.5) reads and fewer than 10⁶ transcripts were kept. Next, only cells in which 85–95% of transcripts belonging to protein-coding genes, 1–3% of transcripts belonging to lncRNA and 5–15% of transcripts belonging to small RNA were kept. Unspliced and spliced protein-coding genes were treated as different entries in our count tables to recover extra granularity in the downstream two-dimensional projection. Potential doublets as detected by Scrublet with default parameters were removed. The resulting count tables were library-size normalized to 10⁴ transcripts, and data were log-transformed with a pseudo-count equal to 1. Cells with a total transcript count to histone genes above 35 were assumed to be in S-phase (Fig. 3). Differential gene expression analysis between cells in S-phase and not S-phase was performed using the t-test to determine cell cycle genes (default scanpy.tl.rank_genes_groups function in Scanpy), for separate timepoints and all data together (Supplementary Table 4). Next, highly variable genes with mean log expression between 0.0125 and 5 were selected, and cell cycle genes were excluded. Number of counts and cell cycle properties were regressed out (Scanpy function scanpy.pp.regress.out), and data were z-transformed (scanpy.pp.scale). For all timepoints, we selected the top 50 principal components (except for E6.5, for which we selected the first 20). For each timepoint, we constructed a directed graph connecting nearest neighbor cells in the reduced principal components analysis (PCA) space, using the Manhattan metric as previously described³². Initially, for each cell, we identified its ten nearest neighbors. An outgoing edge from cell i to cell j was kept if the distance d_ij was less than the mean + 1.5× s.d. among all the distances connecting ten nearest neighbors. Cells that were not connected to any other cell were filtered out. The directed graph was converted to an undirected graph, and a two-dimensional UMAP was obtained as previously described⁷⁰. We clustered the data using the Leiden algorithm (scanpy.tl.leiden, resolution set to 1) and performed differential gene expression between Leiden clusters using the t-test (default scanpy.tl.rank_genes_groups).

scRNA seq analysis for mouse 10x Chromium libraries and individual timepoints

10x data were analyzed similarly to the VASA-seq data. Here, we kept cells with more than 10^3.5 and fewer than 10⁶ uniquely detected transcripts and with 85–97% protein-coding transcripts. Cell cycle genes were not removed from the set of highly variable genes, and cell cycle regression was not performed. The effect of the libraries was regressed out before Z-score scaling.

Comparison between 10x Chromium and VASA-seq embryo data

For the comparison, only reads mapping at the 80% 3′ end of gene bodies were used to generate count tables for both VASA-seq and 10x Chromium. Only genes expressed in both technologies were used for the comparison. The technology and the number of counts were regressed out from the combined VASA–10x Chromium dataset, and dimensionality reduction was performed by PCA. Manhattan-based distances between cells were calculated in the combined PCA space. Equivalent clusters were defined by fist clustering each dataset for each timepoint independently. Second, for a given cluster and a reference technology (for example, VASA-seq), a background histogram of the distances between cells in that cluster and their corresponding first nearest neighbor in the target technology (for example, 10x Chromium) was obtained. Finally, each cell in the target technology was assigned to the cluster of its nearest neighbor in the reference technology. Cells with low transfer scores were excluded, and equivalent clusters with low numbers of cells in any technology were excluded from the downstream analysis. Equivalent clusters between VASA-seq and 10x Chromium were defined as groups of cells with identical 10x Chromium and VASA cluster assignments. To assign a germ layer to each equivalent cluster, published annotations for the 10x Chromium data²⁴ were used (epiblast: epiblast, primitive streak, anterior primitive streak, caudal epiblast and NMP; ectoderm: ExE ectoderm, caudal neurectoderm, rostral neurectoderm, surface ectoderm, forebrain/midbrain/hindbrain, neural crest and spinal cord; mesoderm: nascent mesoderm, caudal mesoderm, ExE mesoderm, intermediate mesoderm, mesenchyme, mixed mesoderm, paraxial mesoderm, pharyngeal mesoderm, somitic mesoderm and cardiomyocytes; endoderm: allantois, def. endoderm, ExE endoderm, gut, parietal endoderm and visceral endoderm; blood: blood progenitors 1, blood progenitors 2, erythroid1, hematoendothelial progenitors, endothelium, erythroid2 and erythroid3; and PGC: PGC). The prevalent annotation for each equivalent cluster was used.

Master UMAP for VASA-drop mouse embryo data

The master UMAP, where all cells for all timepoints are integrated together, was obtained as previously described³². In brief, we first built a directed graph. For each cell in each timepoint, we found the top 30 nearest neighbors in the subset of cells from the same timepoint and the previous timepoint (cells from E6.5 are only connected to cells from E6.5). To do so, all the cells in the subset are projected to the PCA space of the latest timepoint, and distances are calculated using the Manhattan metric. Next, the undirected graph was extracted and used to project the data to the two-dimensional UMAP.

Expanding the transcriptome annotation

A total of 33,662 demultiplexed and ribo-depleted FASTQ files for each cell were used to reconstruct the transcriptome and quantify AS events. To this end, we implemented a custom computational workflow using Snakemake⁷¹ based on Hisat2/StringTie2 (ref. ⁷²) and additional custom scripts. First, PCR duplicates were removed through a custom Python script that calculates pairwise identity across UMIs for each sequenced read within single cells. Then, reads were grouped by previously obtained Leiden clusters and mapped to the reference mouse genome assembly, version GRCm38, using HISAT2 (ref. ⁷³). We performed the alignments implementing the recommended configuration for HISAT2 and genome indexing to ensure an optimal performance during later steps of the transcriptome assembly⁷⁴.

The alignments for each cluster were assembled and then merged using StringTie2 (ref. ⁷²). The resulting GTF file was then compared to the input transcriptome annotation using gtfcompare⁷², which assigns a classification code to each assembled transcript, which is subsequently used to filter transcripts with codes that indicate additional portions of annotated transcripts or novel genes. Novel transcripts spanning three or more exons that were classified under code k, m, n, j, x, i or y were appended to the input transcriptome annotation, expanding the original set of annotated transcripts. Finally, to further improve the quality of potentially novel transcripts, additional custom filtering steps were implemented to avoid novel transcripts due to false-positive novel exons. This filter is particularly important for transcripts assembled from reads that are mapped to repetitive sequences or exons that are ≤30 nt, which can arise from HISAT2 misalignments. To annotate potentially novel microexons, we used MicroExonator, a specialized computational workflow for discovering and quantifying microexons³⁵. After running MicroExonator’s discovery module, we obtained a transcriptome annotation, which was later processed with custom scripts to limit the number of alternative transcription start and end sites.

Quantification of AS events across cell types

The final GTF from the expanded transcriptome annotation was used to quantify isoforms and AS events using Whippet³⁶. We ran Whippet through MicroExonator’s downstream module to profile AS events using scRNA-seq data, which enabled randomized aggregations of cells into pseudo-bulks and pairwise comparisons of AS profiles across cell types. To determine relevant pairwise comparison of AS profiles across cell types, we used PAGA³⁷ to calculate connectivities between cell clusters based on gene expression. We then compared the 72 pairs of clusters that have a connectivity ≥0.05. For each comparison, cells from each cluster were randomly pooled to form at least three different pseudo-bulks of 200 or fewer cells. To detect reproducible changes of splicing node inclusion across cell types, random pseudo-bulk pooling and differential inclusion steps were repeated 50 times for each pairwise comparison, avoiding the detection of spurious splicing events. As part of MicroExonator’s workflow, the obtained probabilities of each splicing node to be differentially included were used to fit a beta distribution model and calculate CDF-beta values for each event. DISNs were defined as events with CDF-beta values equal to or lower than 0.05. To identify SNMs, we calculated the average ψ values for each splicing node across three randomly defined pseudo-bulk samples for each cell cluster. For splicing nodes where ψ values could be quantified based on at least ten reads across at least 50 pseudo-bulks, we calculated the Z-score by comparing to all other pseudo-bulks. We considered a splicing node as an SNM for a given cell type if at least two pseudo-bulks had significant Z-scores (P ≤ 0.05) and an absolute difference of at least 0.3 from the mean across all pseudo-bulks. To show some functional consequences of detected AS events for protein function, we used the drawProteins package⁷⁵ to draw scaled diagrams of protein domains and other features annotated in UniProt⁷⁶.

Reporting summary

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