Single-cell RNA-seq with Waterfall Reveals Molecular Cascades underlying Adult Neurogenesis

Waterfall: Analyzing Single-Cell Datasets from Continuous in vivo Processes

We next examined the whole-transcriptome dataset of individual CFP^nuc+ NPCs. Unsupervised hierarchical clustering analysis resulted in two super-groups with six sub-groups (Figure 2A). Notably, these six CFP^nuc+ groups were not clearly segregated on the plot of principal component analysis (PCA; Figure 2B), which was consistent over different batches of sequencing runs with multiple biological replicates (Figure S2A; Table S1). The continuous trajectory called for a new approach not relying on segmentation into a few groups of cell clusters. We could not use currently available single-cell analysis software, such as Monocle (Trapnell et al., 2010) or Wanderlust (Bendall et al., 2014), for our system due to the lack of sufficient prior information, such as temporal delineators or a robust set of specific markers (See Supplementary Methods). We thus developed a more generally applicable pipeline of algorithms to perform unbiased statistical analyses of multidimensional single-cell datasets from continuous biological processes. We collectively named the suite of algorithms “Waterfall”, which involves three steps: pre-processing, pseudotime reconstruction, and gene expression analysis (Figure 2C; See Experimental Procedures and Supplementary Methods).

Open in a separate window

Figure 2.

Waterfall for analyzing single-cell data from continuous in vivo process.

(A) Unsupervised clustering analysis of CFP^nuc+ NPCs resulting in two super-groups with six subgroups. Different groups are color coded in the same fashion in this figure and across all other figures.

(B) Principal component analysis (PCA) plot shows one of the possible linear trajectories of different groups with the exception of SA.

(C) A schematic diagram of multiple components and workflow of Waterfall. Waterfall is a full range of algorithms for processing multi-dimensional single-cell datasets derived from continuous biological processes. Please see Supplementary Methods for more information and Waterfall analyses of other biological systems.

(D) Representative expression profiles of marker genes of adult neurogenesis. Each data point represents the gene expression level of a single cell with color scheme following Figure 2B. Data points are fitted with local polynomial regression fitting (red lines) with 95% confidence interval (gray area). HMM-predicted underlying states are represented as black and yellow squares on the bottom of the graphs.

See also Figures S2 and S3

Pre-processing defined the trajectory of interest following dimensionality reduction of the data. Unsupervised learning identified six clusters of cells (Figure 2A), which were then labeled S1–5 and SA based on their relative location in a PCA plot (Figure 2B). SA was recognized as a branch by the minimum spanning tree (MST) algorithm (See Supplementary Methods). Although characterizing SA would be interesting (See Supplementary Methods), we focused on the major neurogenic pathway in the current study. The expression profiles of a few known developmental genes were used to orient the most probable trajectory of interest (Figures 2B and S2B).

To reconstruct the chronology, we first determined the most probable route of transcriptomic progression. We performed k-means clustering of single-cell transcriptomes on the PCA plot after excluding SA, followed by constructing a MST trajectory to connect cluster centers (Figure S3A). We then introduced “pseudotime” (Trapnell et al., 2014) to define the relative location of each cell on the MST trajectory (Figure S3A). Taking the Euclidian distance defined by the whole transcriptomic difference from each cell to the next in pseudotime, we found that the total path length reconstructed by Waterfall was significantly shorter than would result from random ordering of cells (See Supplementary Methods). The pseudotime algorithm reconstructs molecular state transitions of a continuous process by quantifying the gradual divergence of single-cell transcriptomes individually, rather than as members of pre-classified groups.

For gene expression analysis, we developed an algorithm to determine the binary on/high or off/low expression state of each gene along pseudotime in an unbiased fashion using a hidden Markov model (HMM; Figure S3B). Gene expression at the single-cell level is highly stochastic and binary (Muramoto et al., 2012; Novick and Weiner, 1957; Raj et al., 2006). When analyzing continuous processes represented by single-cell transcriptomes in the absence of discrete groups, conventional statistical methods, such as arithmetic mean or t-test, are not appropriate. We adopted HMM to statistically convert stochastic expression patterns of individual genes into binary on/high or off/low states (Figure S3B). The binary gene expression states were then shown as heat maps, which quantified the molecular cascade over time (Figures 2C-​-D).D). To discover novel developmentally regulated genes, we correlated gene expression levels with pseudotime and subjected identified genes to gene ontology (GO) analyses (Figure 2C).

Validation for the Reconstructed Adult Neurogenesis Process

We validated molecular dynamics revealed by Waterfall at multiple levels. First, known NSC markers Gfap and Apoe, and eIPC markers Sox11 and Tbr2, showed non-overlapping expression states over developmental pseudotime (Figure 2D; Table S2). Second, we evaluated in vivo expression of Aldoc and Stmn1, which have not been previously studied in adult hippocampal neurogenesis. Over pseudotime, Aldoc was initially highly expressed (on/high state), then downregulated (off/low state), whereas Stmn1 was initially off, then upregulated (Figures 3A and ​and3C).3C). In the hippocampal dentate subgranular zone (SGZ) of adult Nes-GFP^cyto mice (Mignone et al., 2004), Aldoc⁺GFP⁺ precursors were almost exclusively PCNA^– qNSCs, with very few PCNA⁺ active NSCs (aNSCs) or IPCs (Figure 3B). In contrast, Stmn1⁺GFP⁺ precursors were mostly eIPCs and PCNA⁺ aNSCs, but not PCNA^– qNSCs (Figure 3D). Thus, Waterfall accurately predicted the in vivo expression dynamics of both known and unknown genes. Third, for functional validation, we explored the possibility of genetic labeling of a specific developmental stage based on Waterfall results. Hopx (Hop homeobox) was highly expressed in qNSCs, but downregulated around the transition point from qNPC to aNSC (Figure 3E). Upon a single low-dose tamoxifen injection into a Hopx-CreER^T2 mouse line (Takeda et al., 2011) crossed with a mT/mG^f/f reporter line for clonal lineage-tracing (Bonaguidi et al., 2011), almost all labeled precursors at three days post injection were nestin⁺GFAP⁺ qNSCs in the adult SGZ (Figure 3F). By 7 days, we were able to observe GFP-labeled clones that contained both NSCs and their progeny (Figure 3G), indicating self-renewal and differentiation, two hallmarks of NSCs.

Open in a separate window

Figure 3.

Validation for Waterfall predictions for early adult neurogenesis.

(A-D) Validation of gene expression patterns and on/off binary states of Aldolase C (Aldoc, A) and Stmn1 (C) over pseudotime by immunohistology. Also shown are sample confocal images of GFP, cell proliferation marker PCNA, Aldoc or Stmn1 in the dentate gyrus of adult Nestin-GFP^cyco mice (left panels) and quantifications (right panels). Values represent mean ± SEM (n = 3 animals). The pie chart represents the proportion of Aldoc⁺ or Stmn1⁺ cells among each category of the GFP⁺ progenitor population. Scale bars: 20 µm (left) and 10 µm (right).

(E-G) Validation of gene expression patterns and on/off binary states of Hopx over pseudotime (E) by genetic labeling and lineage-tracing. Adult Hopx-CreER^T2::mT/mG mice were injected with a single dose of tamoxifen and examined 3 (F) or 7 days (G) later. Shown in (F) are sample confocal images of GFP, GFAP, Nestin and DAPI. Also shown is quantification of percentages of GFP⁺ cells as NSCs or IPCs. Values represent mean ± SEM (n = 5 dentate gyri). Shown in (G) is an example of a labeled clone containing an NSC and multiple Tbr2⁺ neuronal progeny. Scale bars: 20 µm (left) and 10 µm (right).

Transcription Factor Expression during Adult NSC Activation and Neurogenesis

Waterfall allows an unbiased prediction of the relative chronological position of each individual cell and distribution of binary gene expression over the developmental trajectory. To delineate molecular cascades underlying adult qNSC activation and neurogenesis, we generated a list of the top 1000 negatively correlated genes with pseudotime (DOWN¹⁰⁰⁰ genes; Spearman correlation coefficient < –0.13; Table S4), which represent qNSC-enriched genes down-regulated during activation and neurogenesis, as well as the top 1000 positively correlated genes with pseudotime (UP¹⁰⁰⁰ genes; Spearman correlation coefficient > 0.20; Table S4), which represent newly activated genes during qNSC activation and early neurogenesis.

Out of these 2000 genes, we initially focused on transcription factors (TFs). Most known TFs involved in adult neurogenesis were discovered by extrapolating findings from embryonic studies. In contrast, our database provides unbiased genome-wide profiles of TF expression. Systematic analyses of our dataset revealed a total of 41 down-regulated TFs and 42 up-regulated TFs during adult hippocampal neurogenesis (Figures 4A and S4; Table S5).

Open in a separate window

Figure 4.

Molecular cascade underlying adult quiescent neural stem cell activation and neurogenesis initiation

(A) Lists of DOWN and UP TFs and their Spearman correlation coefficient with pseudotime.

(B) ON/HIGH (yellow) or OFF/LOW (black) states of top 150 DOWN (left) and UP (right) genes sorted by the timing of transition points. Shown on the top are histograms of the numbers of individual cells examined along the pseudotime progression. The colors on the histogram follow the color scheme of Figure 2B.

See also Figure S4.

First, the set of dynamic TFs we identified included known regulators of adult NSCs and neurogenesis, which provided additional validation of our approach. Among DOWN TFs, Sox2, Sox9, Id3, nuclear receptor Nr2e1/Tlx, and Hes1 have been shown to regulate adult NSC maintenance and function. Among UP TFs, SoxC (Sox4 and Sox11), Foxg1, Tbr2, Insm1, Tcf12 and Nfib are critical in proliferative adult NPCs (Table S2).

Second, we identified multiple dynamic TFs that are regulators of embryonic neurogenesis but have not yet been studied in adult neurogenesis. DOWN TFs include homeobox protein Dbx2, nuclear glucocorticoid receptor Nr3c1 and Id4. UP TFs included chromatin protein Hmgb1 and proto-oncogene N-myc. During embryonic neurogenesis, these DOWN TFs inhibit cell cycle (Nr3c1) or prevent premature differentiation (Id4), whereas UP TFs regulate progenitor proliferation (Hmgb1, N-myc), suggesting conserved functions during embryonic and adult neurogenesis (Table S2).

Third, more than half of UP TFs and DOWN TFs are largely uncharacterized in the context of neurogenesis, but many of them are close paralogs or binding partners to other neurogenesis-related genes, or have been implicated in other somatic stem cell systems. Examples include Hmgb1 paralogs (Hmgb2, Hmgb3, Hmga1-rs1), SWI/SNF-related Brg1/Smarca4-associated factors (Smarcc1/Baf155, Smarce1/Baf57), and Nfib paralogs (Nfia, Nfix). In addition, Mxd3, Zeb2 and ZT3/Zfp, regulate adipocyte, melanocyte and myogenic differentiation, respectively, whereas Tsc22d3/Gilz inhibits myogenic differentiation (Table S2). Hopx, which we found to mark qNSCs in the adult dentate gyrus (Figures 3E-​-G),G), is expressed in intestinal stem cells and a subset of multipotent hair follicle stem cells (Table S2).

Together, our single-cell transcriptome datasets are a rich resource of genome-wide dynamic expression profiles of TFs during adult neurogenesis. Many TFs that were previously uncharacterized in adult neurogenesis are known to be involved in embryonic neurogenesis or regulation of other somatic stem cells, suggesting shared biology among different stem cell systems and the potential utility of our resource for the general stem cell field.

Molecular Cascades underlying Adult qNSC Activation and Neurogenesis Initiation

The vast majority of UP¹⁰⁰⁰ and DOWN¹⁰⁰⁰ genes in our dataset were not TFs (Table S4). We investigated their characteristics from three perspectives: transition patterns along the developmental trajectory, cellular location of gene products, and biological function.

For transition patterns, plots of the top 150 genes each from UP¹⁰⁰⁰ and DOWN¹⁰⁰⁰ lists showed a wave of molecular activation or inactivation events over time, highlighting the sequential transition of gene expression during qNSC activation and neurogenesis (Figure 4B). To obtain biological insight into these transition patterns, we performed multiple gene ontology analyses. Strikingly, the predicted cellular localizations of protein products of UP¹⁰⁰⁰ genes and DOWN¹⁰⁰⁰ genes were drastically different. 51% of DOWN¹⁰⁰⁰ genes, as opposed to 20% of UP¹⁰⁰⁰ genes, encode proteins associated with the membrane (Figure 5A, p = 2.6 × 10⁻²⁹). On the other hand, 58% of UP¹⁰⁰⁰ genes, as opposed to 20% of DOWN¹⁰⁰⁰ genes, encode proteins associated with the nucleus (Figure 5A, p = 1.2 × 10⁻³⁶). Similar results were obtained using different thresholds for generating lists of UP and DOWN genes (Figure S5A).

Open in a separate window

Figure 5.

Functional characterization of UP¹⁰⁰⁰ genes and DOWN¹⁰⁰⁰ genes.

(A) Quantification of predicted cellular location of gene products of UP¹⁰⁰⁰ genes and DOWN¹⁰⁰⁰ genes (two middle panels). Also shown are numbers of genes with indicated functions for membrane-associated DOWN genes (left panel) and those for nucleus-associated UP genes (right panel).

(B) Functional GO analysis for DOWN¹⁰⁰⁰ and UP¹⁰⁰⁰ genes. Color of each functional entity represents proportion of UP genes (blue) and DOWN genes (red). Connections between each pair of data points represent sharing more than 5 genes between the pair. Functionally similar entities are grouped with same background colors. Broken lines represent two categories of DOWN genes: genes encoding proteins involved in the intra- or extra-cellular communication and genes encoding proteins defining intrinsic stem cell properties. The P values are from hypergeometric tests, and corrected by Holm–Bonferroni method.

See also Figures S5, S6 and S7.

Molecular Signatures of Adult qNSCs Revealed by DOWN¹⁰⁰⁰ Genes

Functional annotation of DOWN¹⁰⁰⁰ membrane genes revealed enrichment for ion or protein transport, cell communication and cell adhesion (Figure 5B). Further classification identified proteins specific to the plasma membrane, endoplasmic reticulum, Golgi apparatus, and cytoplasmic vesicles (Figure S5B). KEGG pathway analysis revealed diverse functional entities involved in intra- and inter-cellular communication (Figure 5B). Specifically, Notch signaling, GABAergic synapses, glutamatergic synapses, BMP pathways, MAPK pathway, calcium and cell adhesion-related genes were down-regulated upon qNSC exit from quiescence (Figures 5B, ​,6A6A-​-BB and S6A). Electrophysiological recordings of Nestin-GFP^cyto+ NSCs in acute slices from adult animals showed responses to both AMPA and NMDA, suggesting expression of functional receptors (Figure 6B). Each functional signaling pathway entity contained key genes that encode receptors, subunits or downstream mediators (Figures 6A-​-B).B). Notably, many ligands for these receptors, including glutamate, GABA, Wnts, BDNF/neurotrophin, Jagged1, BMPs, FGFs and Insulin/IGF2, are known to be present in the adult SGZ niche (Table S2), suggesting active signaling in qNSCs. While previous studies have examined each of these ligands and receptors in regulation of adult neurogenesis in isolation, our systematic genome-wide analyses unified disparate information and suggested a novel model that quiescent adult NSCs are not passive or dormant, but instead are actively integrating various niche signals. More surprisingly, qNSC activation was associated with decreased expression of genes involved in transducing local environment cues and pervasive down-regulation of various signaling pathway-related genes (Figure 6A). These results suggest that, once activated, adult NSCs shunt their capacity to respond to external regulation.

Open in a separate window

Figure 6.

Sequential molecular dynamics during adult neural stem cell activation and neurogenesis.

(A) Gradual down-regulation of averaged binary states of each functional entity defining intrinsic stem cell properties (top) and those defining intra- or extracellular communication (bottom). On/off binary states of expressed genes within each functional entity were averaged and normalized to show the timing of the transition.

(B) Number of genes in different DOWN gene ontology groups and representative example genes. Expression patterns of representative genes over developmental pseudotime are shown in Figure S6A. Electrophysiological recording of Nestin-GFP^cyto+ NSCs in acute hippocampal slices showed responses to both AMPA and NMDA (bottom).

(C) Gradual up-regulation of averaged binary states of each functional entity defining cell cycle checkpoints and cytoplasmic ribosomal subunits. On/off binary states of up-regulated genes within each functional entity are averaged and normalized to exemplify the timing of the transition. Cell cycle checkpoint genes were up-regulated in the sequence of the cell cycle progression. The up-regulation of cytoplasmic ribosomal subunits preceded up-regulation of the earliest cell cycle checkpoint gene.

See also Figure S6.

KEGG analysis of DOWN¹⁰⁰⁰ genes also revealed a shift in energy source and metabolism. First, multiple lipid metabolism-related functional entities, including fatty acid degradation and sphingolipid metabolism, were enriched in qNSCs, but down-regulated upon activation (Figure 5B). As previously reported (Knobloch et al., 2013), qNSCs exhibited the highest level of Spot14 (Figure S6B), which regulates lipid metabolism. qNSCs also maintained an active fatty acid degradation pathway (Acsl3, Acsl6, and Acsbg1; Figure S6B; Table S2). Second, pathway analysis consistently indicated glutathione metabolism and glycolysis as an adult qNSC characteristic, which was lost upon activation. Among glycolysis genes, aldolase A, aldolase C, and Ldhb decreased significantly, whereas most other glycolysis genes including Gapdh did not change during initiation of neurogenesis (Figure S6C; Table S2).

To validate results from analyses of the top 1000 significantly down-regulated genes, which contained a limited number of genes in each particular pathway, we performed analysis using all expressed genes and an independent functional annotation database wikipathway (Pico et al., 2008). Virtually identical results were obtained (Figure S7; Table S7). Together, analyses of down-regulated genes provided novel insight into molecular signatures of adult qNSCs, including both intrinsic properties and regulation of intra- or inter-cellular signaling pathways.

Sequential Molecular Dynamics during Adult Neurogenesis Revealed by UP¹⁰⁰⁰ Genes

We next analyzed UP¹⁰⁰⁰ genes, which were nucleus-associated and/or related to cell cycle, DNA/RNA metabolism and chromosome organization (Figure 5A). Detailed analysis revealed pervasive activation of cell cycle-related genes, ranging from cell cycle supporting genes, such as nucleotide synthesis, protein/RNA synthesis, and DNA fidelity controls (DNA repair and p53 signaling pathways), to genes directly involved in cell cycle, such as DNA replication, kinetochore complex, cyclin/cyclin-dependent kinases, or cytosolic mitotic spindle genes (Figure 5B).

As opposed to down-regulation of glycolysis-related genes, oxidative phosphorylation-related genes were up-regulated (Figure S6C). Specifically, in contrast to stable expression of earlier mitochondrial respiratory chain complexes (complex I, II, III and IV), expression of subsequent complexes (complex V) increased over pseudotime, implying a gradual completion of the full electron transport chain during neurogenesis (Figure S6C).

The high resolution of Waterfall analyses revealed temporal relationships among genes in different functional groups. Cell cycle checkpoint genes were sequentially activated following the known biological sequence of cell cycles: G1 to S transition, followed by G2 to M transition and then chromosomal segregation, indicating that our pseudotime accurately reconstructs sequential biological events (Figure S6D). Initiation of cell cycle preceded the major transcriptomic shift (Figure 6C). Notably, up-regulation of genes encoding ribosomal subunits preceded the appearance of any cell cycle checkpoint genes (Figures 6C and S6D), suggesting that priming of protein synthesis machinery may mark the G0 to G1 transition ahead of cell cycle entry during adult qNSC activation.

Together, analyses of UP¹⁰⁰⁰ genes suggested that molecular dynamics of qNSC activation and initiation of neurogenesis are largely defined by priming of protein synthesis machinery, cell cycle entry, activation of RNA and protein biogenesis, and a shift in energy metabolism from glycolysis to oxidative phosphorylation. An independent approach using a different functional annotation database showed similar results (Figure S7; Table S7).

A Resource of in vivo Single-cell Transcriptomes of qNSCs and their Immediate Progeny

Our study provides a single-cell RNA-seq dataset and the first comprehensive view of transcriptome dynamics underlying adult qNSC behavior in vivo. Currently, there is no published dataset for transcriptome dynamics during stem cell development in any somatic system in vivo. We performed multiple levels of validation of our dataset and approach, including comparison with an in situ database, confirmation with known and unknown marker expression during adult neurogenesis in vivo, and functional validation via clonal lineage tracing and electrophysiology (Figures 3 and 6B; Table S3).

Our resource of holistic molecular profiles during early neurogenesishas three unique features that increase its versatility. First, our whole-transcriptome information includes unannotated transcripts, isoforms, and retrotransposon-derived transcripts, as opposed to multiplexed qPCR or microarray-based studies which can only provide limited annotated transcripts (Hoppe et al., 2014). Second, we animated static single-cell transcriptomes over the in vivo neurogenesis trajectory, allowing queries of molecular dynamics for each gene over the continuous biological process and generation of novel hypotheses during specific phases of adult qNSC maintenance, activation and neurogenesis initiation. For example, expression of nuclear glucocorticoid receptor Nr3c1 in adult qNSCs but not in eIPCs (Figure 4A) suggests a cellular target of glucocorticoids during adult hippocampal neurogenesis and a means to manipulate qNSCs in vivo. Our resource also reveals potential prospective markers of adult neurogenesis, such as using Hopx-CreER^T2 to target adult qNSCs in vivo (Figures 3E-​-G).G). Third, each transcriptome in our dataset represents a biological state within a single cell. This modular construction allows for flexible reorganization of the dataset to probe different questions as our understanding of the biology evolves. For example, instead of generating reporter lines for individual genes or using different surface markers for physical sorting of specific cell populations, investigators can perform unlimited in silico cell sorting with any individual gene, or multiple genes in combination, to obtain a selected cell population to probe their gene expression characteristics at the genome-wide level using our single-cell datasets.

Developmental Dynamics of Adult Neurogenesis at the System Level

Recent genome-wide studies have begun to provide a system-level understanding of in vivo adult NSC biology using marker-defined NPC populations (Bracko et al., 2012; Codega et al., 2014; Kriegstein and Alvarez-Buylla, 2009). Yet previous studies have only acquired snapshots of transciptomes, which limits investigation of developmental dynamics among different cellular states. We co-opted the imperfection of the Nestin-CFP^nuc genetic labeling system to collect individual qNSCs and their immediate neuronal progeny concurrently. Aligning cells along the developmental trajectory yielded, for the first time, a molecular continuum with sequential progression of the individual transcriptome from qNSC to aNSC and then eIPC. Importantly, this novel approach does not divide developmental processes into discrete stages that are defined a priori by capturing populations sharing specific markers.

Our resources provide unparalleled temporal resolution to identify new mechanisms underlying adult NSC biology. For example, we showed that Acyl-CoA synthetases (Acsl3, Acsl6 and Acsbg1), the enzymes for the first step of fatty acid β-oxidation, were highly expressed only in qNSCs (Figure S6B), suggesting a novel role for active fatty acid β-oxidation in qNSCs and thereby extending previous findings on the role of fatty acid metabolism in adult neurogenesis (Knobloch et al., 2013) (Table S2). We also found that ribosomal subunits were the first genes up-regulated upon adult qNSC activation and early differentiation (Figures 6C and S6D), suggesting a possible demarcation of G0 to G1 transition and providing the timing for switches in protein synthesis, the regulation of which is important for somatic stem cell function (Signer et al., 2014). The holistic picture we obtained unifies disparate information and illuminates novel biological themes in stem cell biology. For example, our analysis suggests that qNSCs actively respond to local environmental cues through various signaling pathways, but gradually and globally shut off signaling capacity upon activation. These observations support the concept of a niche wherein mammalian somatic stem cells are tightly controlled by a regulatory microenvironment (Schofield, 1978), and predict that eIPCs are less responsive to environmental input (Berg et al., 2015). This novel biological insight may be applicable to many somatic stem systems defined by stochastic behavior (Simons and Clevers, 2011).

Library Preparation and Sequencing

cDNA amplification followed the previously published SMART protocol (Ramskold et al., 2012). Briefly, the DNase I was first inactivated by increasing the temperature (75°C for 10 minutes) and samples were then stored on ice. Custom designed 2A oligo 1 µl (12 µM, Integrated DNA Technologies, sequence shown in Figure S1A) was added and annealed to the polyadenylated RNA by increasing temperature (75°C for 3 minutes) and quenching on ice. A mixture of 2 µl Superscript II First-Strand Buffer (5X, Invitrogen), 1 µl custom designed TS oligo (12 µM, Integrated DNA Technologies, Figure S1A), 0.3 µl MgCl2 (200 mM, Sigma), 0.5 µl RNase inhibitor (Neb), 1 µl dNTP (10 mM each, Thermo), 0.25 μl DTT (100 mM, Invitrogen), and 1 µl Superscript II (200 U/μl, Invitrogen), were incubated at 42°C for 90 minutes, followed by enzyme inactivation at 75°C for 10 minutes. A mixture of 29 µl Water, 5 µl Advantage2 taq polymerase buffer, 2 µl dNTP (10 mM each, Thermo), 2 µl custom designed PCR primer (12 µM, Integrated DNA Technologies, Figure S1A), 2 µl Advantage2 taq polymerase was directly added to the reverse transcription product and the amplification was performed for 19 cycles. The amplification product was purified using Ampure XP beads (Beckman-Coulter). Library preparation was performed using Ovation Ultralow library systems (Nugen inc). Libraries were multiplexed and sequenced using Illumina Hiseq 2500 (Illumina Inc) (Table S1).

Bioinformatic Analyses

We thank S.L. Salzberg, H.H. Kazazian, and B. Langmead for suggestions, members of Song and Ming laboratories for discussion, Y. Cai and L. Liu for technical support. This work was supported by MSCRF to H.S. and G.L.M., fellowship from Samsung to J. Shin, and EMBO long-term postdoctoral fellowship and Swedish Research Council to D.A.B.