Femtosecond laser microdissection for isolation of regenerating C. elegans neurons for single-cell RNA sequencing

Our understanding of nerve regeneration can be enhanced by delineating its underlying molecular activities at single-neuron resolution in model organisms such as Caenorhabditis elegans. Existing cell isolation techniques cannot isolate neurons with specific regeneration phenotypes from C. elegans. We present femtosecond laser microdissection (fs-LM), a single-cell isolation method that dissects specific cells directly from living tissue by leveraging the micrometer-scale precision of fs-laser ablation. We show that fs-LM facilitates sensitive and specific gene expression profiling by single-cell RNA sequencing (scRNA-seq), while mitigating the stress-related transcriptional artifacts induced by tissue dissociation. scRNA-seq of fs-LM isolated regenerating neurons revealed transcriptional programs that are correlated with either successful or failed regeneration in wild-type and dlk-1 (0) animals, respectively. This method also allowed studying heterogeneity displayed by the same type of neuron and found gene modules with expression patterns correlated with axon regrowth rate. Our results establish fs-LM as a spatially resolved single-cell isolation method for phenotype-to-genotype mapping. Femtosecond laser microdissection enables transcriptomic analyses of single neurons based on their phenotype.

Our understanding of nerve regeneration can be enhanced by delineating its underlying molecular activities at single-neuron resolution in model organisms such as Caenorhabditis elegans. Existing cell isolation techniques cannot isolate neurons with specific regeneration phenotypes from C. elegans. We present femtosecond laser microdissection (fs-LM), a single-cell isolation method that dissects specific cells directly from living tissue by leveraging the micrometer-scale precision of fs-laser ablation. We show that fs-LM facilitates sensitive and specific gene expression profiling by single-cell RNA sequencing (scRNA-seq), while mitigating the stress-related transcriptional artifacts induced by tissue dissociation. scRNA-seq of fs-LM isolated regenerating neurons revealed transcriptional programs that are correlated with either successful or failed regeneration in wild-type and dlk-1 (0) animals, respectively. This method also allowed studying heterogeneity displayed by the same type of neuron and found gene modules with expression patterns correlated with axon regrowth rate. Our results establish fs-LM as a spatially resolved single-cell isolation method for phenotype-to-genotype mapping.
Spinal cord injuries result in permanent functional deficits as axons in the adult central nervous system fail to regenerate after trauma 1 . State-of-the-art interventions promote neuron survival and axon regrowth by engaging neuron-intrinsic and proregenerative pathways 2 . However, the clinical outcomes of such interventions are still poor: regrowth can only be stimulated in a small number of neurons, while a seemingly homogeneous neuron population can show diverse or even opposite responses to the same intervention 3,4 . A deeper understanding of nerve regeneration can be attained by dissecting the process at single-neuron resolution in a well-defined nervous system. For such studies, Caenorhabditis elegans has been established as a valuable model organism in conjunction with femtosecond laser axotomy 5,6 , which can reproducibly induce a variety of relevant regeneration phenotypes in vivo. The transparent body of C. elegans further allows analysis of such phenotypes, including axon regrowth rate, guidance and fusion at single-neuron resolution in vivo with fluorescence microscopy. Past nerve regeneration studies in C. elegans have led to the discovery of conserved nerve regeneration pathways 7-26 , notably Article https://doi.org/10.1038/s41592-023-01804-3 precise ablation process can mitigate both collateral and out-of-focus tissue damage as shown in our previous work with fs-laser axotomy in C. elegans 5,6,[42][43][44] .
To isolate a single neuron from C. elegans, fs-LM resects the target neuron by ablating its surrounding tissue at a series of spots arranged in a three-dimensional hemispherical shell encasing the neuron (Fig. 1a,b and Extended Data Fig. 1). The ablation spots are spaced 2-4 µm apart from each other and the target neuron. At each spot, we deliver a train of low-energy laser pulses to ablate a small volume (diameter <2 µm) of tissue (Supplementary Note 3). However, cutting the neuronal processes so close to the neuron might cause extensive loss of cellular content. Therefore, we first perform axotomy at a safe distance of 10-20 µm away from the soma and wait a couple of minutes until the axotomized processes are retracted and sealed 6,45 . This step allows the target neuronal soma to keep its membrane fully intact during the fs-LM isolation process. To release the resected neuron, we then ablate the cuticle of the nematode at approximately 5 µm away from the neuron to create a 2 µm wide incision. The neuron is subsequently released through the incision by the driving force of the internal pressure of C. elegans (Fig. 1a,b and Extended Data Fig. 2a-d).
We visually inspect each released neuron for adherent debris, which can be removed by laser-induced bubbles as needed. We then collect the released neuron with a micropipette and wash it in cleaning buffer twice to remove any remaining contaminant in the aspirate. Finally, we deposit the cleaned neuron into a lysis buffer and freeze it on dry ice.
To demonstrate the feasibility of fs-LM, we isolated PLM neurons from living larval stage (L4) wild-type and dlk-1 mutant (dlk-1(0)) C. elegans and studied their gene expression profiles using scRNA-seq. A typical fs-LM took 4-5 minutes per neuron with a yield of 32% (n = 384, Extended Data Fig. 1a and Supplementary Note 4). We used cytoplasmic green fluorescent protein (GFP) as an indicator of neuron structural integrity and only collected neurons that retained their GFP signal level after being released ( Fig. 1b and Extended Data Fig. 2a-d). In the case of neuronal damage, which primarily occurred as neurons became torn or fragmented during the release process due to incomplete resection, the GFP signal of the neuron dropped drastically (Extended Data Figs. 2e-l). We constructed scRNA-seq libraries from 94.0% of the fs-LM isolated PLM neurons (n = 123) that passed the complementary DNA (cDNA) quality check. The failed samples were mostly a result of neuron lysis while still in the micropipette and were excluded (Extended Data Fig. 3). For comparison, we also collected pooled and single TRNs with the currently prevalent isolation method, dissociation-FACS (Extended Data Fig. 4) 32 . We attained sufficient sequencing depth for 90% of the sequenced single neurons (Extended Data Fig. 5a-c).
The detected average number of genes (2,101 ± 83, mean ± s.e.m.) in fs-LM isolated PLM neurons (n = 46 uninjured neurons, Supplementary Note 5) were similar to those detected (2,150 ± 46, mean ± s.e.m.) in single TRNs (n = 23) isolated by the dissociation-FACS method with a threshold of 1 count per million (CPM) reads ( Fig. 1c and Extended Data Fig. 5d). The average number of genes detected in three FACS-isolated bulk samples (4,763 genes, CPM ≥ 1) agreed well with the previous literature 39 . Pairwise comparison between bulk samples and single PLM neurons identified an average of 1,682 overlapping unique genes, accounting for 35.3% of all expressed genes in bulk (Fig. 1d), which includes all TRNs' genes and thus can be considered as a conservative, upper limit.
Tissue enrichment analysis of the genes detected in fs-LM isolated PLM neurons identified TRNs as the top enriched term (−log 10 (tissue enrichment score) = 21.9 ± 1.39 for n = 5, Extended Data Fig. 5e) and found several TRN-specific marker genes, such as mec-17, mec-12, mec-7, mec-10 and mec-18 (Extended Data Fig. 6). The dissociation-FACS-isolated neurons showed ring interneuron, nerve ring and AVK as the top three terms (Extended Data Fig. 5f), likely due to collection of all six TRNs. The results showed that a minimum of 40 fs-LM isolated cells is needed to recapitulate the number of genes the dual leucine zipper kinase (DLK-1/p38) signaling cascade, whose role in neural development, regeneration and degeneration has been confirmed in numerous vertebrate and invertebrate species 7,[27][28][29][30][31] .
Nerve regeneration research in C. elegans has relied on costly and time-consuming mutant screens that test individual genes for function in axon regeneration (Supplementary Note 1). The rapidly evolving high-throughput sequencing methods offer an opportunity to accelerate progress. Nevertheless, isolating regenerating neurons from C. elegans has so far remained infeasible. As in other models, the prevalent cell isolation method for C. elegans entails fluorescence-activated cell sorting (FACS) preceded by lengthy and disruptive tissue dissociation (dissociation-FACS method) 32 . Several drawbacks, notably low yield of cells (<0.5%) 33 , induction of transcriptional artifacts 34,35 and an inability to isolate specific neurons based on regeneration phenotypes, make this method inapplicable (Supplementary Note 2). Other cell isolation methods such as laser-capture microdissection (LCM) 36 , cutting needles [37][38][39] and Patch-seq 40 can bypass tissue dissociation, but may cause intolerable contamination 41 and degradation of cellular content, especially considering the small size of C. elegans neurons (<3 µm in diameter with <0.1 pg total RNA).
The above challenges motivated us to develop femtosecond laser microdissection (fs-LM) as a method that enables the resection of intact single cells directly from living tissue. The micrometer-scale precision of fs-laser ablation allows us to isolate specific neuronal somas from living C. elegans without inducing any artifacts to preserve cellular RNA. The method also allows us to precisely cut the associated neuronal processes and wait until they are fully sealed before isolating the neuron. fs-LM was used by multiple users to successfully isolate various types of neuron from the ventral nerve cord, head and tail regions with different levels of tissue complexity. The cell yield of fs-LM was 32%, an improvement over the dissociation-FACS method having a cell yield of roughly 0.5%. fs-LM mitigated the stress-related transcriptional artifacts induced by the tissue dissociation protocol. We used fs-LM to isolate one of the six touch receptor neurons (TRNs) of C. elegans, posterior lateral microtubule (PLM, either left or right) neuron, exhibiting successful or failed regeneration in wild-type or DLK-1 knockout animals following fs-laser axotomy. Single-cell RNA sequencing (scRNA-seq) of the isolated neurons revealed transcriptional programs that might be correlated with axon regeneration and their potential upstream regulators. To investigate the molecular basis of regeneration heterogeneity displayed by neurons of the same type, we also isolated single regenerating PLM neurons with different axon regrowth rates. We identified genes whose expression patterns correlated with axon regrowth rate at single-neuron resolution. Weighted gene coexpression network analysis further revealed that accelerated axon regrowth correlates with elevated gene expression profiles associated with functions such as neurogenesis, DNA repair and RNA transcription, as well as reduced gene expression profiles involved in immune response and RNA metabolism. Our results establish that fs-LM can isolate intact cell bodies from living tissue on basis of the specific phenotype of interest, such as level of regeneration success, regrowth length and so on. Such capability uniquely positions fs-LM as a precise spatially resolved single-cell isolation method that can facilitate phenotype-driven studies for deciphering the complexity of tissue.

fs-LM isolates single intact neuronal somas from living C. elegans
Conventional LCM of single cells from living tissue is challenging as out-of-focus and collateral damage from laser ablation can damage the target cells (Supplementary Note 3). The challenge is even greater when isolating single neurons from C. elegans due to the small size of its neurons (<3 µm in diameter). Here, we leveraged the precision of fs-laser microsurgery, in which high peak intensity of tightly focused fs-laser pulses allows for submicrometer-scale tissue ablation. Such a Article https://doi.org/10.1038/s41592-023-01804-3 detected in a bulk sample (Fig. 1e), consistent with previous studies 46,47 . Moreover, gene expression profiles exhibited a good correlation between the ensemble of single neurons and the bulk samples (Pearson correlation coefficient of 0.96, Fig. 1f) as well as single TRNs isolated by the dissociation-FACS method (Pearson correlation coefficient of 0.99, Extended Data Fig. 5g). Our results confirm that fs-LM yields intact single neuronal somas for gene expression profiling including those with low expression levels.
To show that the fs-LM method can isolate neurons located within complex and crowded tissues, we isolated GFP-labeled neurons from head ganglion (cephalic sensilla (CEP) and RME neurons) and ventral nerve cord (VD and ALM neurons) using multiple C. elegans strains. The isolated neuronal somas looked fully intact as they could retain the full strength of their fluorescence signal during the entire process (Extended Data Fig. 7 and Supplementary Videos 1 and 2). Furthermore, we confirmed their cell-specific gene expression levels by single-cell PCR reverse transcription (RT-PCR).

fs-LM mitigates transcriptional artifacts in scRNA-seq
Past studies have found induction of stress response genes in cells isolated by the dissociation-FACS method 34,35 . Such artifacts are a result of the tissue dissociation protocol, which involves prolonged enzymatic digestion at elevated temperatures, and can affect data interpretation, obscuring identification of lowly expressed genes and/or study-specific stress genes. fs-LM can mitigate such artifacts by minimizing stress to the target cell thanks to the high precision of fs-laser ablation. In addition, as early transcriptional responses peak at 25 min following stimuli 48,49 , the relatively short fs-LM protocol The dissociation-FACS-isolated neurons were labeled as ALM, AVM, PLM and PVM, as identified using the neuron-specific CeNGEN data 50 . UNK represents an unknown neuron. h, Volcano plot showing DEGs with blue highlighting the substantially down-and up-regulated genes (FDR < 0.05) in dissociation-FACSisolated neurons. The list shows representative GO terms associated with the upregulated genes, including key functions related to stress response (highlighted in red). i, Correlation among gene expression profiles of single neurons isolated by fs-LM or dissociation-FACS methods or predicted-PLM neurons isolated by the dissociation-FACS method, showing a statistically higher correlation among fs-LM isolated neurons (two-sided Wilcoxon rank sum test, P = 2.87 × 10 −5 between fs-LM and predicted-PLM neurons for the top ten genes). Data are presented as mean ± 95% CI.
Article https://doi.org/10.1038/s41592-023-01804-3 (<5 min) can be completed before the artifacts become substantial. We compared gene expression profiles in PLM neurons isolated by fs-LM and TRNs isolated by dissociation-FACS methods using principal component analysis (Fig. 1g). As dissociation-FACS data encompassed all six TRNs, we classified the data into four groups of neurons (anterior, ALM and AVM, and posterior, PLM and PVM) based on their cell-specific gene expression profiles from the CeNGEN (C. elegans neuronal gene expression map and network) database 50 (Extended Data Fig. 8a). PLM neurons isolated by the different methods formed two distinct clusters ( Fig. 1g and Extended Data Fig. 8b), indicating that the transcriptome analysis was dependent on the isolation method.
Next, we investigated differentially expressed genes (DEGs) between neurons isolated by these two methods ( Fig. 1h and Supplementary Fig. 1a). Using gene expression profiles of fs-LM isolated neurons as a reference, we found 121 up-regulated and four down-regulated genes with a false discovery rate (FDR) <0.05, (Supplementary Table 1). Gene Ontology (GO) analysis of the up-regulated genes identified stress-related terms such as response to heat (FDR = 8.3 × 10 −7 ), cellular response to unfolded proteins (FDR = 3.5 × 10 −5 ) and mitochondrial electron transport chain (FDR = 1.2 × 10 −4 ). Among the top up-regulated genes, we identified genes that were found to be up-regulated during heat shock and mitochondrial unfolded responses (Supplementary Fig. 1b-d). The down-regulated genes included mec-7, mec-12 and mec-17, with known functions in microtubule formation and stability in C. elegans 51 . These genes were found to be down-regulated when the microtubules were disrupted by colchicine 28,52 . The results above indicate that fs-LM can alleviate stress-related transcriptional artifacts observed in the dissociation-FACS method.

fs-LM improves reproducibility of scRNA-seq
To examine whether fs-LM can improve the reproducibility of scRNA-seq compared to the dissociation-FACS method, we compared the PLM data obtained by the two methods. The fs-LM isolated neurons showed a significantly higher correlation in expression levels (two-sided Wilcoxon rank sum test, P = 5.38 × 10 −79 ), especially among the top ten highly expressed genes ( Fig. 1i and Extended Data Fig. 8c-e). We expected this result as fs-LM can avoid technical variability introduced by stress-related transcriptional artifacts. Moreover, fs-LM can isolate a specific neuron of interest, while the dissociation-FACS method isolates all labeled neurons and relies on gene expression profiles to identify specific neurons of interest. Together, the results above indicate that fs-LM can mitigate transcriptional artifacts induced by the dissociation-FACS method and improve the reproducibility of scRNA-seq.

fs-LM facilitates profiling of transcriptional dynamics
With fs-LM, we were able to isolate regenerating C. elegans neurons, which has not been possible with existing methods. In wild-type animals, axotomized PLM neurons start regrowing from their proximal stump at 4-6 hours post-axotomy, and result in distinct regeneration outcomes at a 24-hour time point (Fig. 2a,b and Extended Data Fig. 9). To avoid confounding the gene expression induced by regeneration with the immediate injury response to axotomy 53 , we isolated regenerating neurons 24 hours post-axotomy. In wild-type animals, we only collected neurons that regrew but displayed no reconnection to the distal fragment (roughly 56% of the regrowing PLM neurons, Extended Data Fig. 9a), since potential axonal fusion can trigger restorative processes altering the neurons' transcriptional state 54 . We also isolated axotomized PLM neurons from dlk-1 (0) mutants 24 hours post-axotomy The red dotted box shows the region of interest in the tail region where the axotomy was performed on one of the PLM neuronal processes. b, Microscope images of axotomized PLM neurons, located close to the C. elegans tail, immediately after axotomy and 24 hours post-axotomy, with regrowth and no reconnection (n = 18), regrowth and reconnection (n = 16) and regrowth and fusion (n = 9) phenotypes in wild-type (WT) animals and no-regrowth phenotype in dlk-1 (0) mutants (n = 13), observed in at least three independent batches. Yellow asterisks indicate sites of axotomy, and the yellow arrowheads indicate sites of reconnection of the proximal fragment (P, neuronal process connected to the cell body) to the distal fragment (D, neuronal process disconnected from the cell body). Scale bar, 10 µm. c, Unsupervised hierarchical clustering (top) and pseudotemporal ordering (bottom) of three types of neuron based on 1,000 genes with the highest variance in expression levels; uninjured neurons from both wild-type and dlk-1 (0) animals, axotomized neurons from wild-type animals with regrowth and no reconnection, and axotomized neurons from dlk-1 (0) animals with no-regrowth. We combined the uninjured neurons from both wild-type and dlk-1 (0) animals since they showed few DEGs. d, Overlap between DEGs in successful or failed regeneration. Blue circles (regrowth) represent DEGs between regrowing and uninjured neurons from wild-type animals. Green circles (no-regrowth) represent DEGs between no-regrowing and uninjured neurons from dlk-1 (0) animals. Article https://doi.org/10.1038/s41592-023-01804-3 from animals that failed to initiate regrowth by forming a growth cone (84% of the axotomized neurons 43 , Extended Data Fig. 9). We, therefore, refer to data from regrowing but not reconnecting neurons in wild-type animals as 'regrowth' and data from dlk-1 (0) animals of not regrowing neurons as 'no-regrowth' in the following. As controls, we isolated uninjured PLM neurons from both wild-type and dlk-1 (0) animals that underwent mock axotomy.
To understand the transcriptional dynamics leading to successful regeneration in wild-type animals and regeneration failure in dlk-1 (0) animals, we performed scRNA-seq on the collected neurons. Unsupervised hierarchical clustering showed that the uninjured and axotomized neurons formed two clusters independent of their genetic background, indicating that the observed transcriptional activities were dominated by axon regrowth instead of DLK-1 knockdown (Fig. 2c).
Further pseudotemporal ordering by Monocle2 (ref. 55 ) positioned the neurons into three branches according to their regeneration states (Fig. 2c). The S1 branch comprised mostly uninjured neurons from both wild-type and dlk-1 (0) animals, the S2 branch comprised mostly regrowing neurons from wild-type animals and the S3 branch included most of the nonregrowing neurons from dlk-1 (0) animals. The single branching point suggested a transition from an uninjured state into two distinct transcriptional states of successful or failed regeneration.
To identify genetic players involved in the regeneration process, we profiled the DEGs between axotomized and uninjured neurons from wild-type and dlk-1 (0) animals (Extended Data Fig. 10a). Between the regrowing and uninjured neurons from wild-type animals, we identified 1,007 DEGs (fold change >2, FDR < 0.05), including 427 up-regulated and 580 down-regulated genes (Fig. 2d, Extended Data Fig. 10b and  Table 2). Between the nonregrowing and uninjured neurons from dlk-1 (0) animals, we identified 510 DEGs, 399 of which were down-regulated (Extended Data Fig. 10c). Because both populations underwent a similar axotomy process, the nonoverlapping DEGs might be candidates relevant to the regeneration processes whereas the overlapping DEGs are expected to be involved in the stress response to the regeneration process. The 152 overlapping DEGs between these two sets exhibited a significant correlation in fold change of expression levels (Pearson correlation coefficient of 0.94, P = 2 × 10 −16 ), suggesting that these DEGs were likely induced by the stress response to axonal injury or loss of neuronal function instead of promoting the regeneration itself ( Supplementary Fig. 2a and Supplementary Note 5). The nonoverlapping DEGs were likely involved in axon regeneration, with their selective induction possibly leading to successful or failed axon regeneration. GO analysis followed by network visualization of these DEGs indeed revealed major network hubs associated with relevant terms such as the cytoskeleton, synapse signaling and development (Fig. 3a), which were also identified in regenerating mammalian neurons 56 .

Discovery of genetic programs driving axon regrowth
To further understand the regulatory mechanisms involved in axon regeneration, we investigated the links between the identified DEGs and the known regeneration mechanisms in C. elegans. Among roughly 300 known regeneration regulators in C. elegans (Supplementary Table 3), we identified 48 regulator genes overlapping with the DEGs (P = 0.0015, chi-squared test with Yates correction; Supplementary Table 4); 39 required and nine inhibitory regulators for axon regeneration. Most of the required regulators were up-regulated in regrowing neurons (30 out of 39, P < 0.0001, chi-square test) and down-regulated in nonregrowing neurons (23 out of 39, P = 0.0024, chi-square test). These regulators also exhibited previously unknown expression patterns (Fig. 3b, Supplementary Fig. 2b and Supplementary Notes 5 and 6).
To explore the upstream regulators of the identified DEGs, we screened for transcription factors with strong binding affinities to the DEGs using the chromatin immunoprecipitation with sequencing (ChIP-seq) data from modENCODE and modERN repositories 57,58 ( Supplementary Fig. 3). We identified 37 transcription factors (P < 0.01, Kolmogorov-Smirnov-like test) among the 204 screened C. elegans transcription factors (Supplementary Table 5 Apart from the known regulators, we also identified transcription factors with high binding affinity upstream of the promoter region and potential regulatory roles in axon regeneration ( Fig. 3d and Supplementary Fig. 3e). These results demonstrate that fs-LM can facilitate the discovery of insights on the regulatory mechanisms of nerve regeneration in C. elegans. Further functional studies are required for the validation of these candidate genes.

Regrowth heterogeneity displayed by the same type of neurons
In both C. elegans and mammalian models, neurons display heterogeneity in regeneration capacity depending on the neuronal type, age, preconditioning lesion, or injury condition 59 . Previous investigation into the molecular basis of such heterogeneity has largely focused on differential gene expression displayed by different neuronal types with varying regeneration capacities [60][61][62] . Nevertheless, regeneration heterogeneity also manifests among neurons of the same type 63,64 . In C. elegans, the PLM neurons display regeneration heterogeneity hallmarked by a variable axon regrowth length (Fig. 4a,b), even when animal age and injury conditions are tightly controlled. Such regeneration heterogeneity is likely induced by random fluctuations in the intrinsic genetic machinery of axon regeneration. For example, previous work has shown that the amplitude of calcium influx and release from internal storage following axon injury can modulate the axon regrowth rate 8 . Unfortunately, this aspect of nerve regeneration has remained understudied due to a lack of methods. With fs-LM, we isolated regenerating PLM neurons based on the regrowth length of their axon at 24 hours post-axotomy. We combined parametric and nonparametric tests to identify 226 candidate genes whose expression patterns correlated with axon regrowth rate (Supplementary Fig. 4), with most of them exhibiting a negative correlation (144 out of 226, Supplementary Table 6). We further analyzed the coexpression patterns among the identified genes using weighted gene coexpression analysis 65 . We identified six gene modules with distinct expression patterns and biological functions (Fig. 4c). GO analysis of each module revealed that accelerated axon regrowth is associated with elevated expressions of genes located in modules 1-3 with functions related to neural development (Fig. 5a-c and Supplementary Note 6). Compared to module 1, module 2 displayed lower expression at a regrowth rate <5 µm h −1 , whereas module 3 showed higher overall expression yet smaller change with increasing axon regrowth rate.
By contrast, genes in modules 4-6 displayed reduced expression levels in fast regrowing neurons (Fig. 5d-f and Supplementary Note 6). While genes in module 4 showed gradually decreasing expression levels with increasing regrowth rates, the expression of genes in modules 5 and 6 showed a steep decrease at 2 and 6 µm h −1 , respectively. Genes in these modules had minimal overlap with the DEGs identified when comparing regrowing and uninjured neurons in the previous section (Fig. 5g). GO analysis of all modules indicated a role of the calcium signaling pathway, corroborating previous findings 8 ( Supplementary  Fig. 4c). This result highlights the power of fs-LM to isolate specific neurons based on a phenotype of interest. Correlation network analysis allowed us to study the molecular basis of regeneration heterogeneity   displayed by the isolated neurons of the same type and identify genetic players with roles in regeneration biology. Beyond C. elegans, fs-LM can also be readily applied to other tissues such as acute mice brain slices ( Supplementary Fig. 5 and Supplementary Note 5).

Discussion
In this work, we present fs-LM as a method for isolating specific single cells with specific phenotypes and from a living organism with complex surrounding tissue architecture. The high precision of fs-laser ablation minimizes collateral and out-of-focus damage, allowing us to resect intact single neuronal somas with their content unchanged. This study shows that fs-LM can overcome technical limitations inherent to the state-of-the-art cell isolation methods. Specifically, when compared to the prevalent dissociation-FACS method, fs-LM can offer a few advantages: (1) fs-LM can isolate a specific neuron among all labeled neurons in a given strain, (2) it can also isolate unlabeled cells, (3) it can associate cells with their specific phenotype, (4) it can reduce artifacts (such as stressed-induced gene expressions) and (5) it allows collecting rare cells thanks to its 35% efficiency. Dissociation-FACS with its 0.2% efficiency requires dissociating thousands of worms to isolate only a few neurons, which would have been infeasible for nerve regeneration studies.
Among other existing methods, LCM can retain contextual and phenotypic information of the isolated cells 36 similar to fs-LM, but it requires a series of tissue processing steps (freezing, thin slicing, fixing, staining, dehydration and so on) that pose a high risk of degradation 40 or contamination 38 . As for purification-based isolation methods, such as transcriptome in vivo analysis 66 and translating ribosome affinity purification 67 , only part of cellular content can be captured, which will likely suffer additional losses during column purification. Recent years have seen the rise of in situ sequencing methods 68,69 . fs-LM is complementary to such methods as they are still limited to messenger RNA, while intact cells isolated by fs-LM are amenable to all kinds of assays.
In terms of limitations, fs-LM is relatively labor-intensive and has low throughput (roughly five neurons per hour) compared to the dissociation-FACS method, making it more suitable for precision and phenotype-driven studies. Automation can potentially increase this throughput while affordable laser sources can facilitate a broad adaption of fs-LM method.
With a complete blueprint of anatomy, wiring diagram and function, the C. elegans nervous system has long been one of the most well-defined models in neuroscience. Most recently, the CeNGEN consortium has also produced gene expression profiles of C. elegans neuron classes 50,70 . Combining such a comprehensive knowledge with fs-LM's capabilities will allow studying mechanisms of development, regeneration and degeneration in relation to neuronal genetic composition and physiology in a complete nervous system. In summary, our work establishes fs-LM as a versatile spatially resolved, single-cell isolation method with unique capabilities to facilitate future studies.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41592-023-01804-3.

Methods
Experiments presented in this paper complies with all relevant ethical regulations, approved by Institutional Biosafety Committee and Institutional Animal Care and Use Committee of The University of Texas at Austin for C. elegans and mice studies, respectively.

C. elegans strains and maintenance
We used the following strains in this work: CZ10175 zdIs5 [Pmec-4::GFP+lin-15(+)], CZ11327 dlk-1(ju476); zdIs5, CZ1200 [Punc-25::GFP], BY200 [Pdat-1::GFP] and N2. Both CZ10175 (wild-type animal) and CZ11327 (dlk-1(0)) strains express GFP in the six TRNs 27,71 . The strain CZ1200 expresses GFP in the 26 GABAergic motor neurons 72,73 . The BY200 strain expresses GFP in the six dopaminergic neurons 74 . C elegans were grown and maintained on nematode growth media (NGM) plates with bacteria of strain HB101 at 20 °C as previously described 5 . To obtain age-synchronized populations for FACS experiments, we lysed gravid adults with alkaline hypochlorite solution and collected the eggs in M9 buffer. The eggs were suspended in M9 buffer and allowed to hatch overnight on a rotor for aeration. The hatched larvae were arrested in larval stage 1 (L1) until the food was introduced. The synchronized L1 population was plated on an NGM plate and allowed to grow for 50-52 h until ready for harvest at middle to late larval stage 4.

Laser axotomy for nerve regeneration studies
Femtosecond laser axotomies as an injury mechanism for nerve regeneration studies were performed on C. elegans hermaphrodites anesthetized on agar pads on an upright microscope using a ×63, 1.4 numerical aperture (NA) oil immersion objective (Plan-Apochromat, Zeiss). Agar pads were prepared by sandwiching 0.25 ml of melted 2% agar between two microscope slides. For each batch of axotomy, we used 15 µl of 10 mM levamisole solution to immobilize 8-12 L4 stage animals under a coverslip on the agar pad. We ablated the axon of the PLM neuron that was closer to the coverslip at 50 µm from the soma. For axotomy, we used a train of 300 laser pulses with 7.5 nJ pulse energies that were generated from an amplified fs-laser system (Spectra Physics 'Spitfire', 805 nm center wavelength, 250 fs pulse width, 1 kHz repetition rate). We recovered the axotomized animals in a small droplet of M9 buffer and placed them on the bacteria lawn of a fresh NGM plate. The axotomized animals were allowed to recover at 20 °C before imaging or fs-LM. From immobilizing the first animal of the batch to collecting the last axotomized animal from the agar pad, each batch of axotomy was consistently completed within 7-8 min to avoid stressing the animals. For axon regrowth length measurements, we took z stack images (2 µm spaced) of the regenerating axons. The length measurement was taken on a max-projection image formed from the z stack images by ImageJ (v.1.52 u). Only processes regrowing longer than 5 µm were included in length measurements. We excluded reconnected axons from the statistical analysis. For mock axotomies, ablation was performed on the hypodermis instead of the axon (at 2 µm from the position where actual axotomies would be performed).

C. elegans dissociation and FACS
For dissociation-FACS isolation of C. elegans TRNs, we followed an established protocol for dissociation of C. elegans 32 . We collected a synchronized animal population with M9 buffer from NGM plates and washed the animal pellet with egg buffer three times to remove bacteria. The pellet was then incubated in a freshly prepared 0.25% sodium dodecyl sulfate, and 200 mM dithiothreitol solution for 4 min to weaken the cuticle of the animals. We washed the treated pellet five times with egg buffer, then suspended the pellet in freshly prepared 15 mg ml −1 pronase solution to further digest the cuticle for 20 min. We gently pipetted the entire suspension up and down with a syringe and 18G needle to promote digestion and release of cells. During pipetting, bubbles were carefully avoided. At the end of the pronase treatment, we stopped the reaction by adding a prechilled L15 medium supplemented with 10% fetal bovine serum (FBS) and allowed the undigested carcasses to settle at 4 °C for 5 min. The cell suspension remained chilled from this point forward. We collected the supernatant and washed the collected cells twice with L15 containing 2% FBS. The final concentration-adjusted suspension was filtered through a 5 µm syringe filter. Before FACS, we added propidium iodide to the cell suspension for live/dead staining. FACS was performed using a BD FACSAria III sorter and GFP-labeled C. elegans TRNs were collected using FACSDiva (v.8.0.1, BD Biosciences) (Extended Data Fig. 4). The GFP sorting gate was determined using cell suspension collected from N2 animals as a reference, and the propidium iodide sorting gate was determined by test sorting a stained cell population and observing the separation of live or dead cells into two clusters. The GFP + propidium iodide-negative events were sorted directly into a 96-well plate (single-neuron sample, n = 26) or 1.5 ml tubes (bulk sample each containing >5,000 single neurons, n = 3) containing lysis buffer with 2 U µl −1 RNAse inhibitor. The collected GFP + propidium iodide-negative neurons were 0.3% of all sorted cells. We confirmed the purity of the collected neurons by sorting extra GFP + propidium iodide-negative events into an L15 medium and visually inspecting the collected cells under a microscope. The sorted samples were promptly frozen on dry ice for subsequent reverse transcription of cDNA.

fs-LM of single C. elegans neurons
Before fs-LM, we washed each batch of four animals in M9 buffer to remove bacteria and anesthetized them using 20 mM levamisole for up to 10 min until movements were suppressed. Meanwhile, we taped a clean, hydrophobic coverslip to a microscope slide to provide mechanical support. We pipetted three droplets of three different solutions onto the coverslip: 50 µl of dissection buffer (1-20 mg ml −1 pronase in egg buffer), 25 µl of L15 medium with 15-30% FBS and 25 µl of fresh L15 medium (Extended Data Fig. 1 and Supplementary Notes 4 and 5). The anesthetized animals were transferred into the dissection buffer droplet using an eyelash pick. We then mounted the slide onto an upright microscope, equipped with an XYZ motorized stage (ASI PZ-2000FT). fs-LM was subsequently performed using a 2 mm working distance, ×60, 1.0 NA water dipping objective (Olympus LUMPLFL-N60XW). We first ablated the neuronal processes of the target neuron 10-20 µm away from its soma using 300 laser pulses with 16 nJ pulse energy (Supplementary Note 3). Once the processes were sealed, we ablated the surrounding tissue using 4-6 ablation spots distributed in a hemispherical pattern encasing the neuron. Each ablation spot included a train of 300 laser pulses with 16 nJ pulse energy, generated from the same laser as used for laser axotomy. We slightly increased the pulse energy to 20-25 nJ when ablating spots deeper underneath the target neuron to compensate for scattering. For puncturing the cuticle, 150-300 laser pulses with 50 nJ pulse energy were used. The resected neurons would normally be released through the incision driven by the internal pressure of the animal. The release process could also be assisted by delivering 150-300 laser pulses with 50 nJ pulse energy to a spot inside the worm roughly 50 µm away from the neuron and generating laser-induced bubbles that can push the neuron out of the incision. In case the released neuron was stuck to debris, it could be separated by inducing laser-induced bubbles in the buffer at roughly 50 µm from the neuron, which could generate a torrent in the buffer that effectively breaks up the debris without damaging the neuron.
To aspirate the released neuron, we used a glass micropipette controlled pneumatically by a microinjector (Narishige IM-11-2), which was mounted on a three-axis micromanipulator (Narishige MN-153). The micropipettes had a tip inner diameter of 10-15 µm and were polished by a microforge. We coated the inner surface of the micropipette with silicone (Sigmacote, Sigma-Aldrich). Previous to fs-LM, we filled the tip of the micropipette with 15-30% FBS/L15 buffer, positioned the tip at the center of the field of view and then retracted the tip so that during Nature Methods Article https://doi.org/10.1038/s41592-023-01804-3 fs-LM the tip would not be in the path of the laser. When a neuron was successfully released, we brought the tip close to the dissected neuron while applying low positive pressure to generate a slight outflow. We used the outflow to wash and elevate the released neuron away from other cells or debris, after which we aspirated and collected the neuron from the dissection buffer droplet. The collected neuron was next washed once in the 15-30% FBS/L15 buffer droplet and twice in fresh L15 medium droplets. For each wash, we completely flushed the micropipette or up to the position where the cell was stuck on the pipette surface (without creating bubbles) to guarantee cleanliness. After the last wash, we pipetted 6 µl of lysis buffer prepared from Clontech SMART-seq v4 kit with 2 U µl −1 RNAse inhibitor onto the coverslip. The collected neuron was then deposited into the droplet of lysis buffer, which we aspirated and transferred to a PCR tube for freezing at −80 °C. An identical fs-LM protocol was applied to collect uninjured wild-type (n = 23), uninjured dlk-1 (0) (n = 25), axotomized wild-type (n = 47) and axotomized dlk-1 (0) (n = 25) PLM neurons. The axotomized neurons were collected over several days and with appropriate control populations (uninjured animals). To demonstrate the wide applicability of the fs-LM protocol, we isolated single CEP neurons, ALM neurons and RME neurons using BY200, CZ10175 and CZ1200 strains, respectively. For each animal, only one neuron was isolated by fs-LM.

scRNA-sequencing and data analysis
cDNA libraries were constructed with Clontech SMART-seq v.4 3′ DE Kit per the manufacturer's instruction. cDNA libraries were pooled and next-generation sequencing libraries were constructed by Illumina Nextera XT DNA Library Prep Kit. We performed pair-end 150 bp sequencing with Illumina Hi-seq 4000. The raw sequencing data were demultiplexed and filtered to remove low-quality reads by Cutadapt v.2.3 and v2.8. We aligned the filtered reads to C. elegans reference genome WS271 using STAR aligner v.2.7.2. The alignments were summarized using the featureCount function from Subreads package v.1.6.4. Unique alignments were retained. Each neuron received over 3 million summarized counts, which was deep enough to maximize gene discovery according to subsampling analysis. As Clontech SMART-seq v.4 3′ DE Kit sequences the 3′ end of amplified cDNA fragments, we used CPM as a normalized metric of gene expression. For differential expression analysis, we used single-cell differential expression package v.3.12. For a detailed description of the data analysis discussed in the paper, please refer to Supplementary Note 5.

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

Data availability
The ChIP-seq datasets used for identifying potential regulators of nerve regeneration are publicly available at https://www.encodeproject.org/ (ENCODE Project 75,76   The microdissection is performed with a 60× 1.0 NA water dipping objective. The sample manipulation is performed with a 10× 0.3NA air objective with a 10 mm working distance. A glass microneedle is mounted on a 3-axis manipulator and suction pressure is controlled using a precise pneumatic puller. The pi-chart shows the isolation of n = 123 PLM neurons (32.0% out of n = 384 total PLM neurons). The pie chart also depicts the unsuccessful situations and classifies them under different categories based on the step at which we lost the neuron. Most failures (incomplete resection, n = 88 neurons) occurred primarily when the neurons were torn into pieces during the release process through the cuticle (Extended Data Fig. 2e-l). b, Priming steps for the glass micropipette before it can be used for cell manipulation. c, fs-LM isolation steps including C. elegans immobilization, single-cell dissection, cell collection in the glass micropipette, and cell lysis for single-cell transcriptomic analysis. The resected neuron was released through the incision smoothly without damage, as evidenced by consistent volume and cytoplasmic GFP intensity (n = 123 out of total 384 attempts). e-h, Incomplete resection of neuron causing it to remain in place after the incision was made (n = 88 out of total 384 attempts). As a result, the neuron was torn by the outflow of surrounding tissue. i-l, A resected neuron was fragmented when migrating through a small cuticle incision. In this case, the neuron was moving towards the incision, which indicated complete resection. Scale bar, 50 µm. Fig. 3 | Quantification of amplified cDNA libraries from fs-LM isolated single C. elegans neurons. a, Example of libraries of good quality. Inset: the isolated neurons retained GFP intensity until being deposited into the lysis buffer. b, Example of libraries that failed quality check, and were thus excluded from this study. Inset: the isolated neuron lysed before being deposited into the lysis buffer. Although the cellular content of the neuron supposedly remained in the micropipette, degradation, and attachment to the inner wall of the micropipette affected reverse transcription (RT)-PCR, resulting in libraries dominated by short fragments that were mostly primer dimers. c, Example of NTC controls. Article https://doi.org/10.1038/s41592-023-01804-3 Extended Data Fig. 4 | Isolation of C. elegans touch receptor neurons by the dissociation-FACS method. a-d, FACSAria II gate settings for isolating living GFP-labeled C. elegans touch receptor neurons from dissociated animals. Prior to FACS, the cell suspension was filtered through a 5 µm syringe filter to prevent clogging. SSC-A, SSC-W, and FSC gates were set up to exclude debris and clusters of cells. We stained the dead neurons with propidium iodide. GFP + PI-events were sorted into 96 well plates or 1.5 mL conical tubes containing lysis buffer. The GFP threshold was determined in reference to the autofluorescence of cells from dissociated N2 worms. The PI threshold was determined by test sorting a small number of PI-stained cells and observing two distinct populations of live/ dead cells. e, Histogram of GFP levels of the isolated neurons. f, Cell suspension prior to sorting and g, example of a collected neuron. Green: GFP. Red: PI. Multiple locations were imaged to confirm GFP + neurons using the fluorescent microscope (confirmed with 3 independent sample preparation). Scale bar: 5 µm. h, Quantification of cDNA library prepared from one of the GFP + sorted cells using FACS method. Article https://doi.org/10.1038/s41592-023-01804-3 Extended Data Fig. 7 | Examples of successful fs-LM isolation of single neurons from different locations in multiple C. elegans strains. a, Isolation of one of the PLM neurons from the tail region in CZ10175 animal (we resected n = 14 PLM neurons out of 15 attempts over 6 independent batches). Scale bar, 50 µm. b, Isolation of VD7 neuron from the mid-body region (close to the vulva muscle and residing on the ventral cord region) in the CZ1200 animal (resected n = 3 VD7 neurons over 3 independent batches). c, Isolation of one of the CEP neurons from the head region (in the head ganglion) in the BY200 animal (resected n = 24 CEP neurons out of 30 attempts over 10 independent batches). d, Isolation of one of the RME neurons from the head region (in the head ganglion) in the CZ1200 animal (We resected n = 8 RME neurons out of 9 attempts over 4 independent batches). The top panel shows the neurons' schematics. Below are the fluorescence images describing the steps for the axotomy of neuronal processes, neuron dissection, cuticle ablation, neuron extraction, and neuron collection. The lightening arrows indicate the locations of the laser spots for axotomy and fs-LM (light pink, indicating low pulse energies) and cuticle ablation (red, indicating higher pulse energies required to ablate the cuticle). The yellow arrows indicate the isolated neurons outside the animal body and as collected inside the micropipette tip. The bottom panel shows 2% agarose gel with PCR amplification products of single-cell lysates and visualized with ethidium bromide. Lanes are marked with the target genes or ladder used in the experiment. Fig. 8 | Correlation of gene expression profiles between fs-LM and dissociation-FACS methods. a, Normalized expression of the selected 36 genes, obtained from CeNGEN feature lists, utilized to visualize and identify cell clusters that were assigned a predicted neuronal category. With this classification analysis, we obtained PLM (n = 6), PVM (n = 5), ALM (n = 8), and AVM (n = 3) subgroups with one remaining unclassified. The color panels show the four neuron categories and the normalized level of expression. b, Principal component analysis (PCA) on gene expression profiles in single PLM neurons isolated by fs-LM (fs-LM (PLM), n = 40) and dissociation-FACS (dissociation-FACS (PLM), n = 6) methods. The six dissociation-FACS isolated neurons were classified using the neuron-specific gene expression data available in CeNGEN data. c, Pearson correlation coefficients as calculated from 50 groups of 6 neurons isolated using the fs-LM method (n = 40), predicted-PLM method (dissociation-FACS, n = 6), and all TRNs (dissociation-FACS, n = 23). Shaded areas represent a 95% confidence interval around the mean. d, Box and whisker plot for the top 10 expressed genes. e, Box and whisker plot for all 1,000 ranked genes. The box extends from the 25 th to 75 th percentiles, the line is plotted at the median, the whiskers drawn down to the 10 th percentile and up to the 90 th percentile. Points below and above the whiskers are drawn as individual points. The statistical test used is a two-tailed unpaired Wilcoxon rank sum test. P-value = 0.024 (dissociation FACS and predicted PLM for top 10 genes, *) and P-value < 0.0001 (****).