Detection of transcriptome-wide microRNA–target interactions in single cells with agoTRIBE

Ago2-tagging DNA construct and plasmids

DNA-encoding ADAR2DD_E488Q (human ADAR2 adenosine deaminase domain 316–701-amino acid (aa) hyperactive mutant E488Q, the ADAR-only control) and ADAR2DD_E488Q in-frame with 55-aa flexible linker and Ago2 (agoTRIBE) were chemically synthetized and inserted in pcDNA3.1(+) vector (GeneArt, Invitrogen/Thermofisher). The synthetic coiled-coil (EIAALEK)3, E3- and (KIAALKE)3, K3-tags were described previously³⁷. The E3-GGSG linker-Ago2 and K3-GGSG linker-ADAR2DD_E488Q constructs were synthetized and inserted in pcDNA3.1(+) vector (GeneArt, Invitrogen/Thermofisher). Plasmids encoding EGFP-Ago2 and FLAG/HA-TNRC6B protein were described previously^46,47. eGFP_ADAR2DD_E488Q was synthetized and cloned in pcDNA3.1(+) vector (GeneArt, Invitrogen/Thermofisher). The plasmid encoding GFP nanobody, pGNb-mCherry³⁶ was used as pGNb-vector backbone. GNb-Ago2 and GNb-ADAR2DD_E488Q constructs were made by PCR from the corresponding complementary DNA, followed by insertion into the pGNb-vector lacking mCherry using Gibson assembly cloning (NEB). The following primers were used: Ago2-forward (GGGGGATCTGGATCCATGTACTCGGGAGCCGGCCCCGC);

Ago2-reverse (GCCCTCTAGACTCGAGTCAAGCAAAGTACATGGTGCGC);

ADAR2DD-forward (GGGGGATCTGGATCCATGCAGCTGCATTTACCGCAGG);

ADAR2DD–reverse (GCCCTCTAGACTCGAGTCAGGGCGTGAGTGAGAAC).

The DNA-encoding human T6B peptide (TRNC6B 599–683 aa³⁰) and T6B fusion with GSG linker and mCherry were synthetized and cloned in pcDNA3.1(+) vector (GeneArt, Invitrogen/Thermofisher). The main constructs have been submitted to the Addgene repository with ID nos. 205598 (pcDNA3.1_ADAR2DD_E488Q), 205599 (pcDNA3.1_ADAR2DD_E488Q_Ago2) and 205600 (pcDNA3.1_T6B-mCherry).

Cell culture and transfection

HEK-293T cells were grown in DMEM medium (Sigma) supplemented with 10% FBS (ThermoFisher) in the presence of penicillin/streptomycin (Sigma) under standard conditions at +37 °C in 5% CO₂. Before transfection, cells were seeded in 12-well plates (Corning) at 60–70% confluency. Cells were transfected with 0.5–2.0 μg of the recombinant plasmid DNA in two biological replicates 18–24 h after seeding using Lipofectamine 3000 (Invitrogen) according to the manufacturer’s instructions.

Immunofluorescence staining and microscopy

HEK-293T cells were seeded on four-well glass slides (Nunc Lab-TekII) at 60–70% confluency and transfected with 0.5–1.0 μg of the recombinant plasmid DNA using Lipofectamine 3000 (Invitrogen). The 48-h posttransfected cells were fixed with 3.7% formaldehyde (Sigma), permeabilized with PBS Tween containing 0.1% Triton X-100 (Sigma) and blocked in PBS (Gibco) with 1% bovine serum albumin (Sigma) and 0.1% Tween 20 (VWR). The primary antibodies against Ago2 and ADAR2 were added for 1 h at room temperature. After washing, secondary antibodies were added for 1 h in the dark. Cells were mounted using Vectashield mounting medium (Vector Laboratories) with 0.1 μg ml^–1 DAPI (ThermoScientific) for DNA counterstaining. Cells were imaged with a ZOE fluorescent cell imager (Bio-Rad). The following antibodies were used: rabbit polyclonal human ADAR2 (no. GTX54916, GeneTex), mouse monoclonal Ago2 (no. ab57113, abcam), mouse monoclonal anti-Flag (no. F3165, Sigma), goat anti-rabbit IgG (H + L) Alexa 594 (no. A11012, ThermoFisher) and donkey anti-mouse IgG (H + L) Alexa Fluor Plus 488 (no. A32766, ThermoFisher).

Bulk and single-cell RNA-seq

For standard bulk RNA-seq, cells were harvested 24 and 48 h after transfection. For bulk RNA library preparation, total RNA was isolated using the Quick-RNA miniprep kit (Zymo Research). Quality and integrity of RNA samples were assessed using a Bioanalyzer instrument and the RNA 6000 Nano kit (Agilent Technologies). One microgram of total RNA was used for RNA library preparation with the TruSeq Stranded mRNA kit (Illumina). Sequencing was carried out on an Illumina NextSeq 500. For single-cell RNA library preparation using the Smart-seq3 protocol, cells were harvested 48 h after transfection, stained with 1 μg ml^–1 propidium iodide (Invitrogen) and sorted into 384-well plates using a BD Bioscience FACS Aria fusion instrument. Single-cell Smart-seq3 libraries were prepared at the Eukaryotic Single-Cell Genomics facility and sequenced on an Illumina NovaSeq 6000 platform at the NGI facility (SciLifeLab) as previously described⁴⁰.

Bulk and single-cell RNA-seq data preprocessing and mapping

Both bulk and single-cell RNA-seq data were first processed for quality and adapter trimming using cutadapt v.4.0, allowing a minimum length of 36 base pairs (bp) after trimming (with parameters –q 30 -m 36), and then mapped to the hg38 reference genome using STAR v.2.7.2b (ref. ⁴⁸) with parameters detailed on GitHub (see link below). Putative PCR duplicates were removed using Picard using default options. The six samples used for estimation of transcriptome-wide editing levels (Supplementary Table 1) were subsampled to a sequencing depth of 18 million mapped reads using the Sambamba tool⁴⁹ with the following command line options: view -f bam -t 20–subsampling-seed=3. We used the featureCounts tool⁵⁰ for quantification of the expression level of genes based on Ensembl Hg38 genome annotation (https://www.ensembl.org/info/data/ftp/index.html, release v.99). Multimapping reads and/or reads overlapping by more than one feature (genes in our case) were accounted for using a fractional assignment (-M -O–fraction).

Bulk and single-cell RNA editing analysis

Alignment mismatches were called using a custom variant-calling pipeline (available on GitHub: https://github.com/vaishnoviS/agoTRIBE.git) with two main filtering steps: (1) Phred scale-based quality filtering and (2) read depth-based filtering. All mismatches supported by Phred score >30 and at least two reads supporting the nucleotide substitutions were reported. Single-nucleotide substitutions were then intersected with positions in the single-nucleotide polymorphism (SNP) database (NCBI dbSNP v.151)⁵¹, and all known SNPs were excluded from further analyses. The remaining A-to-G and T-to-C (on reverse strand) substitutions were further compared with editing events compiled in the REDIportal database⁵² and known editing sites were excluded. Strand-aware A-to-G substitutions were then reported as editing events. For the purpose of identification of agoTRIBE-specific targets, the number of editing events was computed for each annotated 3′ UTR in protein-coding genes. Nonsense-mediated decay, nonstop decay, processed transcripts, retained introns and unprocessed pseudogene isoforms were removed. Regression analysis was performed with the editing events in agoTRIBE-transfected samples compared with ADAR-only control samples to identify 3′ UTRs having higher editing. Standard residual values were calculated based on linear regression applied on editing events per gene using the lm function from the R stats package. The top 1,000 3′ UTRs with highest residual values were considered as highly edited 3′ UTRs by agoTRIBE. Conversely, the 1,000 3′ UTRs with lowest residual values were considered as background ADAR-only targets.

Small RNA-seq

Total RNA was isolated using the Quick-RNA miniprep kit (no. R1054, Zymo Research). The quality and integrity of RNA samples were assessed using a Bioanalyzer instrument and the RNA 6000 Nano kit (no. 5067-1511, Agilent Technologies). Small RNA libraries were prepared using the NextFlex small RNA library v.3 kit (no. 5132, BioScientific, Perkin Elmer Applied Genomics). The following amounts of total RNA were used as input: 0.25–1.00 μg of imatinib-treated K562 samples, 0.5 μg of agoTRIBE-transfected and untransfected HEK-293T control and 0.1 μg of cell cycle-specific populations (G0/G1, S, G2/M), and both agoTRIBE-transfected and untransfected HEK-293T control.

Small RNA-seq analysis

Bulk small RNA-seq data were analyzed using miRTrace to trim adapter sequences. Processed reads from miRTrace were then quantified using quantify.pl from the miRDeep2 package. To identify differentially expressed miRNAs in control and agoTRIBE-transfected HEK-293T samples (Supplementary Fig. 2), DESeq2 was used on raw counts for miRNA and all miRNAs with P_adj. < 0.01 were considered significant. For cell cycle analysis (Supplementary Fig. 13), reads per million (RPM) values obtained from quantify.pl were used to identify the top ten miRNAs in each stage and were plotted as a line plot.

Overlaps between miRNA targets reported by agoTRIBE, CLIP-seq methods and TargetScan

Those genes with 3′ UTRs that were specifically edited in agoTRIBE samples (above) were compared with the top target genes from publicly available HITS-CLIP²⁵, PAR-CLIP²⁴ and eCLIP²⁶ data. HITS-CLIP (Ago2 CLIP-Seq) and PAR-CLIP (AGO1234 PAR-CLIP) binding information was downloaded in BED format from the Dorina database⁵³ and eCLIP data (GSE140367, GSM4247216, empty vector control rep 1 (AGO-eCLIP-seq)) were downloaded from NCBI GEO⁵⁴. Targeted genes from HITS-CLIP and PAR-CLIP data were sorted based on total peak coverage, while genes from the eCLIP target was sorted based on P values reported in the GEO supplementary file. For each set, the top 1,000 genes were considered for further comparison with our agoTRIBE top 1,000 target genes. For the purpose of identifying top TargetScan targets, the ten most highly expressed miRNAs in HEK-293T cells were identified (miR-10a-5p, miR-92a-3p, miR-10b-5p, miR-191-5p, miR-148a-3p, miR-16-5p, miR-378a-3p, miR-182-5p, miR-222-3p and miR-186-5p, in descending order by expression) using our miRNA-seq data (Supplementary Table 13). The top 100 targets for each of these miRNAs (based on the cumulative context ++ score) were downloaded from the TargetScan database²⁸. Overall, 699 unique high-confidence target genes were compiled for the top ten most highly expressed miRNAs in HEK-293T cells (given that some targets were shared between miRNAs and were redundant).

Identification of target derepression following T6B transfection

Gene expression was estimated as described above. Unexpressed genes and those with <50 counts following normalization in any of the considered samples were removed. Between-sample normalization of expression was performed using the trimmed mean of the M-values method⁵⁵. Fold change values between control and T6B-mCherry samples were calculated with the foldchange function from the gtools R package and visualized with their empirical cumulative distribution function. The top 699 TargetScan targets were defined as above, while the remaining expressed genes were defined as background.

Identification of positional information in target editing

To identify whether agoTRIBE editing reflects miRNA binding positional information, we compared agoTRIBE editing positions with those reported by AGO-eCLIP–seq data (GSM4247216). To do so we focused on genes with monoexonic 3′ UTRs and their largest isoforms that showed one eCLIP editing event. Genes with 3′ UTRs shorter than 40 bp were discarded. We then calculated the distance (in bp) between agoTRIBE and eCLIP editing events if detected in the same gene. As background control we defined a random iteration of uniform distributions (n = 100) to generate dummy editing positions within each targeted 3′ UTR. Background distances were then computed using eCLIP positional data and randomly distributed editing sites.

Single-cell RNA-seq cell filtering

Smart-seq3 data were preprocessed and mapped as described above. As a preliminary quality control step we filtered out cells with <9,000 detected genes or <150,000 deduplicated reads. To assess the efficiency of transfection we mapped our reads to the linker sequence with bowtie2 v.2.3.5.1 (ref. ⁵⁶) using local alignment (–local) while allowing for one mismatch (–N 1), and quantified the aligned reads with the featureCounts v.2.0.0 tool⁵⁰. For HEK-293T cells we then calculated normalized linker counts per cell as raw linker counts divided by the total number of reads, and plotted these against the editing events in 3′ UTRs. A pseudocount of 1 was added to all linker values to avoid division by zero. By visual inspection of the resulting plot we decided to exclude 163 cells with log₁₀-normalized linker below −4.5 and below 1,000 editing events, leaving 566 cells in total. For K562 differentiated and undifferentiated cells, individual cells with a minimum of one linker count were considered positively transfected and used for further analysis. The K562 single-cell trajectory was constructed using transcriptomics counts with Monocle v.2.28.0 in R, with the top 1,000 genes differentially expressed between stages.

Cell cycle assignment for single cells

Cell cycle annotation was based on the Seurat v.4.0 pipeline⁴¹. Cells were count per million normalized using the NormalizeData function from Seurat, with the normalization method of relative counts (normalization.method = ‘RC’) and a scale factor of 10⁶ (scale.factor = 1,000,000). We then annotated cells into cell cycle phases. To do this, we customized the CellCycleScoring function from Seurat to better fit Smart-seq3 data and used it with the default markers provided by Seurat. Specifically, cells with S-score and G2/M-score (as defined by the Seurat software) >50 were assigned to G1. These same marker genes were considered for scaling of data and calculation of the first 50 principal components with the Seurat functions ScaleData and RunPCA and default parameters. Finally, uniform manifold approximation and projection (UMAP) dimensionality reduction on the first ten dimensions was performed with the RunUMAP function. Based on this graphical representation, every cell cycle phase was split into ‘early’ and ‘late’ substages, with both substages containing approximately the same number of cells. For the purpose of identifying changes in gene expression and miRNA targeting across the cell cycle, sequence data from all cells belonging to a given cell cycle substage were pooled and mapped to Thermofisher ERCC92 spike-in sequences using bowtie, with strict parameters (-v 1) and a pseudobulk approach. The number of reads mapping to the spike-in was divided by the corresponding number of cells in the specific subphase to derive a scaling factor. This scaling factor was then used to normalize expression values of the pooled cells mapping to each cell cycle substage.

Transcript half-life analyses

The top 20% of stable and unstable transcripts was obtained from ref. ⁵⁷ (GEO: GSE49831). For the violin plot (Supplementary Fig. 9), the per-transcript sum of expression and sum editing over all agoTRIBE-transfected cells (540) were calculated. For cell phase analysis (Supplementary Fig. 10), cells from each cell cycle phase were pooled and expression and editing were calculated. Normalization was performed as in Fig. 3i.

Erythroid differentiation from K562 cells

K562 cells (1.0–1.3 × 10⁵ ml^–1) were cultured in 12-well plates containing RPMI (no. 21875034, Gibco) with 10% fetal bovine serum and penicillin/streptomycin supplemented with 5 mM imatinib (no. 9084 S, Cell Signaling) for erythroid differentiation. Cells cultured in RPMI supplemented with 10% FBS, penicillin/streptomycin and 5 mM DMSO (no. 8418, Sigma) were used as an undifferentiated control. The Neon transfection system (no. MPK 1096, Invitrogen) and 1 μg of agoTRIBE and ADAR-only constructs were used for K562 electroporation. Cells harvested 24 h after electroporation and corresponding DMSO (no. D8418, Sigma) controls were used for Smart-seq3 preparation.

Cell sorting by cell cycle stage

An asynchronous population of un- and transfected HEK-293T cells was used for cell cycle analysis and cell sampling for small RNA-seq. HEK-293T cells were transfected with 2.5 μg of agoTRIBE- and ADAR-only constructs using Lipofectamine 3000 and harvested 48 h after transfection. Following cell cycle analysis and cell sorting, cells were stained with Vybrant DyeCycle Green (no. V35004, Invitrogen) and propidium iodide. A BD Bioscience FACS Aria fusion cytometer was used for cell sorting in the G1, S and G2/M phases as determined by Vybrant DyeCycle Green staining.

Cell cycle profiling of agoTRIBE protein levels

HEK-293T cells were transfected with 2.5 μg of the agoTRIBE recombinant plasmid DNA using Lipofectamine 3000, then 48-h posttransfected cells were harvested. The suspended cells were washed with PBS (no. 10010-031, Gibco), fixed with 3.7% formaldehyde (no. 47608, Sigma), permeabilized with PBS Tween containing 0.1% Triton X-100 (no. T8787, Sigma) and blocked in PBS with 10% FBS, 0.36 μM IgG (no. 026202, Life Technologies) and 0.1% Tween 20 (no. 0777, VWR). The primary antibody against ADAR2 (no. GTX54916, GeneTex) was added for 30 min at room temperature and the secondary antibody, goat anti-rabbit IgG Alexa 594 (no. A11012, ThermoFisher), was added for 20 min in the dark. After washing with PBS, cells were resuspended in 400 μl of PBS and stained with Vybrant DyeCycle Green and 30 nM DAPI (no. 62248, ThermoScientific). A negative control with IgG and without ADAR2 primary antibody was used for Alexa 594 compensation, and the agoTRIBE Alexa 594-positive cell population was used for cell cycle analysis. G0/G1, S and G2/M phases were determined using Vybrant DyeCycle Green staining and a BD Bioscience FACS Aria fusion cytometer. agoTRIBE protein expression was determined as mean Alexa 594 fluorescence intensities normalized to G0/G1, S and G2/M event counts.

Statistical analyses

Differences in fold change derepression following T6B transfection in the top 1,000 edited genes by agoTRIBE and in miRNA targets according to HITS-CLIP, compared with background unedited and/or nontargeted genes, were tested with the DTS test, which calculates a reweighted integral of the distance between two empirical cumulative distributions, available in the dts_test function from the twosamples R package. Overlap significance between agoTRIBE targets, CLIP targets and TargetScan targets was estimated using binomial statistics. Differences in long linker reads and in editing events comparing control samples and samples following transfection with agoTRIBE were assessed with a nonparametric approach using the Wilcoxon rank-sum test.

AI protein structure predictions

The human Ago2 (Q9UKV8) three-dimensional protein structure was predicted by AlphaFold2 DB v.1 (ref. ¹⁹). K3/E3 peptides, eGFP and GNb and anti-GFP nanobody (LaG-16) structure predictions were generated using the Robetta, RoseTTAFold¹⁸ modeling method available at https://robetta.bakerlab.org/.

Reporting summary

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