Comprehensive analysis of long noncoding RNA expression in dorsal root ganglion reveals cell-type specificity and dysregulation after nerve injury

Supplemental Digital Content is Available in the Text. Novel and annotated LncRNAs showed marked transcriptional change during neuron differentiation and after nerve injury. Identified LncRNAs were antisense or adjacent to pain genes and ion channels.


Introduction
Dorsal root ganglion (DRG) neurons provide connectivity between peripheral targets and the spinal cord. These neurons show significant heterogeneity in relation to their morphology, functional properties, growth factor dependence, and transcriptional profile. 1,4,63 This reflects the highly specialised nature of these neurons subserving distinct sensory modalities, including temperature, pain, itch, touch, and proprioception. Recent single-cell RNA-seq studies have provided a means to classify these neurons and have identified multiple DRG neuron subgroups. 38,67 Pathologies of DRG neurons-for instance, in the form of acquired or inherited peripheral neuropathies-have a significant impact on human health as a consequence of sensory loss and neuropathic pain 12 and are destined to become more common with an ageing population and increased prevalence of type II diabetes.
A wide variety of injuries applied to sensory neurons, whether traumatic or metabolic, result in marked alterations in transcription of protein-coding genes. 13,34,53 Such changes can have either beneficial or maladaptive outcomes, including increased expression of regeneration-associated genes 9,51 and altered expression of ion channels resulting in enhanced DRG neuronal excitability and neuropathic pain. The DRG has therefore become a model system to study the transcriptional changes after injury. This focus on RNAs' encoding proteins is understandable, given their obvious link with function. However, there are other types of RNA and one of these, long noncoding RNA (LncRNA) has been relatively neglected and little studied in the context of sensory neurones.
LncRNAs are usually multiexonic transcripts of more than 200 base pairs that can modulate gene expression through cis and trans signalling, and have important functional effects. 6,31,58,69,76 The mechanisms by which LncRNAs may alter gene expression are very diverse, including: complementary binding of antisense LncRNAs, transcriptional interference at promoter sites, altered chromatin structure, competing for miRNA binding, and binding to transcription factors. 26,44,81 Ion channels are key determinants of the excitability and hence the functional properties of sensory neurons. Antisense LncRNAs have previously been identified to the voltage-gated potassium channel, KCNA2, 86 and the voltage-gated sodium channel, SCN9A. 32 In the former case, induction of the antisense LncRNA after nerve injury was shown to result in reduced expression of KCNA2 (which acts as an excitability break), leading to sensory neuronal hyperexcitability and the development of neuropathic pain. These are selected examples of functionally relevant LncRNAs that illustrate the important role they can play. However, there has, to date, been no comprehensive analysis of LncRNA expression within the DRG partly because LncRNAs are typically expressed at low levels and are known to vary by species, tissue, and developmental stage. 31,68 Our aim was to use high-coverage RNA-seq combined with a dedicated bioinformatics platform to identify as comprehensively as possible LncRNAs expressed in the DRG. We compared rats and 2 different mouse strains (which show differing degrees of mechanical hypersensitivity after nerve injury 67 ). We wished to determine whether LncRNAs were expressed in a cell-type-specific manner and also to assess the effect of nerve injury on LncRNA expression. In addition, we also assessed LncRNA expression in human Induced Pluripotent Stem Cell -derived sensory neurons.

Data availability
Data are in GSE107182 super series that consists of GSE107180 (rodents' DRG) and GSE107181 (human IPSC). Splicing junctions (SJ), differentially expressed (DE) data analysis results, and GTF files with annotations used are included as downloadable supplementary files of the GEO  All data have been integrated to PainNetworks. In http:// www.painnetworks.org, the user can examine a gene (or set of genes) of interest alongside known interaction partners on the protein level. This information is displayed by the resource in the form of a network. Moreover, the user can access all expression data (log2 fold change and false discovery rate-adjusted P values) and download these in the form of spreadsheets. A tutorial on how to use PainNetworks can be accessed following this link http://www.painnetworks.org/tutorials/RefMan.pdf.
All intergenic and antisense LncRNAs' profiling data are accessible in PainNetworks (http://www.painnetworks.org) → ExpressionData → Mouse centric/Rat centric/Human centric. Experiment names are GB-BALBC-LNCRNAS for the BALB/c mouse, GB-SNI-B10D2-LNCRNAS for the B10.D2 mouse, GB-RAT-LNCRNAS for the rat, and IPSC_HS_AD2-LNCRNAS for IPSC-derived neurons. Naming is as follows: Closest {gene or sense gene}_LNCRNA_{IG or nothing}_chr:start-end(strand All procedures on rats were performed in accordance with UK home office regulations and in line with the Animals Scientific Procedures Act 1986 at a licensed facility at King's College London. Animals were group housed in temperature-and humidity-controlled rooms where food and water was available ad libitum, with a 12-hour light-dark cycle. The welfare of all animals was continually assessed throughout all procedures. In total, 24 rats were used. Rats were humanely culled. L5 DRG tissue from male Wistar rats was collected 21 days after the spinal nerve transection (SNT) surgery, placed into sterile tubes, frozen on dry ice, and stored at 280˚C. Each sample comprises 3 pooled animals, and we had 4 samples of each condition (SNT vs sham).

Mouse
All procedures in mice were performed in McGill University, Montreal, Canada, were approved by the McGill University Animal Care Committee and are fully consistent with Canadian Council on Animal Care guidelines. All mice strains were procured from Jackson Laboratories (Bar Harbor, ME) at 4 to 8 weeks of age. All animals of the same sex were group housed in a vivarium at 21˚C in standard shoebox cages, 2 to 4 per cage, with access to food (Harlan Teklad 8604 each strain, 6 spared nerve injury [SNI]-6 Sham stratified for sex). All operations were performed on adult mice. Brain and DRG tissue has also been dissected from 3 wild-type mice (C57/bl6) and was used to determine relative expression of mRNA using quantitative real-time polymerase chain reaction (qPCR) (see below). In total, 27 mice were used.
The tip of the iliac bone, "the first articular process more than 1 mm rostral to the iliac crest," 61 was used as the landmark for identifying the L5 DRG in all samples. L3 and L4 DRG were dissected from all mice 28 days after peripheral nerve injury. Each sample represents 1 animal and consists of both L3 and L4 DRG. Twelve BALB/c mice and 12 B10.D2 mice stratified for condition and sex were used. All dissections were performed on dry ice, and RNase Decontamination Solution was used to prevent RNA degradation. Tissue was placed into sterile Eppendorf tubes and initially stored on dry ice. For long-term storage, samples were stored in a 280˚C freezer.

IPS-derived human neurons
Human fibroblast-derived IPSC was generated as described previously. 11 Neural differentiation was performed using 8 the protocol with modifications.

Mouse spared nerve injury
The surgical procedure for SNI followed a published protocol developed for rats 15 and adapted for mice. 64 Under general anaesthesia (isoflurane and oxygen), the common peroneal and the sural branch of the sciatic nerve were cut and the tibial branch spared. For sham surgery, the same surgical and anaesthetization procedures were followed, but the nerve branches were simply exposed and not damaged. We assessed mechanical hypersensitivity after SNI surgery on the ipsilateral mouse paw.

Rat spinal nerve transection
The left L5 spinal nerve was ligated and transected, and the L4 and L6 branches were left intact. In sham animals, the spinal nerve was exposed but not ligated.

Behavioural tests
The behavioural test was conducted in a specially allocated room in the animal facility unit at McGill University, performed at a consistent time of the day and by the same experimenter. Mice habituated to the vivarium for at least 1 week before testing. Mechanical pain-related hypersensitivity in mice was assessed using von Frey filaments and the up-down method of Dixon 10 to determine the 50% withdrawal threshold. Mice were first acclimatised to behaviour equipment and baseline behaviour performed 3 times and an average was calculated before surgery. Baseline paw withdrawal threshold was 1.27 g (SD 5 0.22) for BALB/c strain and 1.36 g (SD 5 0.23) for B10.D2 strain (N 5 12). Mice were assigned to the sham or SNI group randomly and postinjury mechanical sensitivity tested at day 1, 7, 14, 21, and 28 (N 5 12 per strain stratified for sex and condition, 6 SNI-6 Sham mice per strain). Assuming an effect of 30% and an SD 5 20%, we need an N 5 6 to achieve power 5 80 at an a 5 0.05 twosided 1-way analysis of variance (ANOVA). Mice from both strains and surgery groups were tested on the same day. The experimenter was not informed about the condition (injury vs sham) of animals but could not be blinded because of the coat colour of the different strains.

RNA isolation and library preparation
RNA was extracted using a hybrid method of phenol extraction (TriPure; Roche, Welwyn Garden City, United Kingdom) and combined with column purification (High Pure RNA tissue Kit; Roche). 14 Dorsal root ganglion samples were first homogenised in TriPure using a handheld homogeniser (Cole-Palmer, Saint Neots, United Kingdom). For IPS cells, TriPure was added directly to the well after removal of media. The concentration of RNA in the samples was measured using a nanodrop. Total RNA was provided to the sequencing centre, and the ribodepleted fraction was selected for further sequencing. In rats, this was the polyadenylated fraction. It was then converted to cDNA using the strand-specific deoxy-UTP strand-marking protocol.

Sequencing and mapping
All samples were sequenced at the Oxford Genomics Centre. Sequencing was performed using the Illumina HiSeq4000 pairedend protocol with 100 bp reads for the mouse DRG, 75 bp for human IPSC/neurons, and Illumina HiSeq2000-100 bp reads for the rat DRG.
The DRG from 24 mice (12 per strain stratified for sex and condition) was sent for sequencing. During library preparation, 2 samples (sample 72 BALB/c SNI Male and sample 68 B10.D2 SNI Female) were accidentally mixed together and destroyed. From the 22 samples sent for sequencing, 2 were excluded (sample 59 B10.D2 SHAM Male and sample 66 BALB/c SNI Female) because of having more ambiguously mapped reads, lower percentage of mapped reads, and higher Cook's distance than all the other samples.
Mapping to the genome was done using STAR aligner. 18 Reads were mapped on the mm10 mouse genome, rn6 rat genome, and Hg38 human genome, all downloaded from ENSEMBL. Conditions and strains were multiplexed in lanes and library batches. Lanes were merged as BAM files after mapping. 39

Differentially expressed and counting features
Differentially expressed analysis was performed using DESeq2 42 default settings. Significant cutoff in all cases was false discovery rate-adjusted P value ,0.05. Counting of features was done using HTSeq 3 and the intersection and not empty strategy to resolve ambiguously counted reads.
All visualisations used regularised log2-transformed counts. 42 Principal component analysis was always performed on regularized log-transformed counts using the top 10,000 genes in mice and humans, and 5000 genes in rats ranked by their SD. Hierarchical clustering was done on regularised log2 counts of the whole gene set, using Euclidean distances and complete linkage.
Gene Ontology (GO) enrichment for DE genes was carried in R using top GO and GSEA. 2,49 Enrichment of neuron subtypespecific genes was calculated with the Fisher exact test, and enrichment of Biological Process (BP) in network modules was calculated using hypergeometric distribution.

Identification of novel LncRNAs
We used a customised reference-based transcript assembly pipeline that requires a reference genome and gene set Only properly paired and uniquely mapped reads were selected. We selected SJs covered with .2 reads and with lengths .20 and ,100,000. We discarded all reads overlapping annotated gene models. We then used the remaining subset of RNA-seq reads to identify I.o.E outside known gene models using a coverage window approach. Gene models were extended by 1000 bp in each direction to ensure that elongated untranslated regions (UTRs) or not yet annotated exons would not be considered putative novel genes. We selected continuous regions above the coverage threshold of more than a read-mate length to ensure that overlapping read mates would not artificially increase coverage. For I.o.E length $100 and depth .2, I.o.E were identified using the function "BAM_to_IOE." Islands of expression were collapsed and clustered as co-overlapping features connected by SJs. A connectivity matrix was created holding all interconnecting I.o.E. In each cluster, consensus introns were calculated by the relative frequency of each discrete segment of a set of SJs. We then subtracted the genomic intervals of these consensus introns from the genomic intervals of the grouped (I.o.E) to reconstruct full-length putative LncRNAs. For novel I.o.E with no overlapping SJs, we first selected only the intersect of the respective genomic regions across all samples. Then, monoexonic putative LncRNAs were kept for further analysis only if the length-normalised coverage had Pr (.jZj) , 0.1. Coverage across I.o.E was fed into a smoothed z-score signal processing algorithm. Z-score thresholding was used to identify introns not identified by the aligner and sudden coverage drops indicating end of transcription activity. 5 Rolling coverage was calculated over a smoothing window of 31 bp, and the minimum coverage drop threshold was set to 5 and the minimum intron length to 20 bp. We only kept novel intronic genes if they were supported by evidence of novel splicing junction and did not contain retained introns.
We included putative LncRNAs in this novel annotation only if they were present in all replicates of a biological condition or strain. Annotations were exported in the Gene Transfer Format (GTF). Subsequently, we filtered out transcripts with length ,200 bp, and we used CPAT 72 to assess coding potential. An average expression cutoff threshold similar to Ref. 56 of .0.5 fpkm for at least one condition was applied to novel LncRNAs performed for downstream analysis. The pipeline was scripted in R 59 using bioconductor 21 packages and custom scripts. All iterative processes were executed in parallel to optimise run times using parallel and BiocParallel. 48,50,56 More details in supplementary methods are (available at http:// links.lww.com/PAIN/A676). All scripts of the workflow are available in github: http://github.com/gbaskozos/ Scripts_LncRNAs.

Transcription start sites mapping to mm10
Transcription start sites (TSS) data were downloaded from FANTOM 5 database. 20,40 We downloaded TSS data that have been classified as "True TSS" by the "TSS classifier." The UCSC Lift Over tool 46 has been used to translate genomic coordinates from the mm9 genome to the mm10. Fifty-one percent of the true TSS were unambiguously mapped to mm10.

Tissue specificity
Tissue specificity was calculated using the tau metric 77 :

Primer design
Primers for the detection of LncRNA and reference gene expression were designed using Primer-BLAST (https:// www.ncbi.nlm.nih.gov/tools/primer-blast/). Primers were designed not to overlap any other annotated gene. Thus, primers designed for antisense LncRNAs were not able to detect regions complementary to sense gene exons. Primer efficiency and specificity were validated before experimental use.

Reverse transcription polymerase chain reaction
For qPCR analysis, RNA samples were converted into cDNA using Evoscript Universal cDNA Master kit (Roche) and by following the manufacturer's instructions. This kit uses random primers. For LncRNA2754, the primers designed also detected a putative UTR, and therefore, a strand-specific real-time reaction was used. Strand-specific RT primers (see table 18, available at http:// links.lww.com/PAIN/A676) were used for LncRNA2754 and for HPRT1 (final concentration 0.5 mM), and 200 ng of RNA was used for each reaction. Strand-specific reverse translation into cDNA synthesis was performed using Transcriptor reverse transcriptase (Roche) and dNTPs (Promega, Southampton, United Kingdom).

Quantitative real-time polymerase chain reaction
For qPCR analysis, cDNA (5 ng) and primer pairs (1 mM, see table 18, available at http://links.lww.com/PAIN/A676) were mixed with LightCycler 480 SYBR Green Master (Roche) in a 1:1 ratio and added to white 384-well plates (Roche). Plates were run on a 45cycle protocol using the LC 480 II system (Roche). Gene expression for each mouse target primer was normalized against the reference gene HPRT1, and the relative expression (delta CT) was calculated. For human IPS cells, transcript expression was normalised against the average CT of GAPDH and HPRT1. For each target, transcript expression is shown relative to a control group (eg, Sham). Significance was calculated using ANOVA with a design of ; sex 1 condition for each mouse strain and ; condition for human IPSC averaging over all cell lines. N 5 10 per strain (6 Sham-4 SNI for BALB/c, 5 Sham-5 SNI for B10.D2; for strandspecific qPCR, B10.D2 strain N 58, 5 Sham-3 SNI). Quantitative real-time polymerase chain reaction was also used to assess relative mRNA expression in brain vs DRG, N 5 3. Significance was calculated using 1-way ANOVA.

In situ hybridisation
In situ hybridization was performed as previously shown. 14 Once cut, sections were air-dried onto superfrost slides in the cryostat for 0.5 hours and then stored in the 280˚C freezer. In situ hybridization was performed using the RNAScope 2.5 RED chromogenic assay kit and by following the manufacturer's instructions (Advanced Cell Diagnostics, Newark, CA). Briefly, slides containing tissue sections were removed from the 280˚C freezer, allowed to equilibrate to RT and rehydrated in phosphatebuffered saline (PBS). Pretreatment required a hydrogen peroxide step at RT, followed by an antigen retrieval step and protease treatment in a hybridization oven at 40˚C. Slides were then incubated with the target or control probes at 40˚C for 2 hours. For HAGLR mRNA, the probes were designed to target position 1153 to 2443 of NR_110445.1. After probe incubation, slides were subjected to 6 rounds of amplification, and the probe signal was developed via a reaction with fast red. To combine with immunohistochemistry (IHC), tissue sections were then washed with PBS-Tx (0.3%) and treated with either the isolectin B4 (IB4) conjugated to biotin (Sigma, Gillingham, United Kingdom, 1:100) or primary antibodies against NF200 (mouse anti-NF200, Sigma, 1: 250) and CGRP (rabbit anti-CGRP; Peninsula Labs, San Carlos, CA, 1:1000) overnight at room temperature. IB4 and primary antibodies were diluted in PBS-Tx (0.3%). Slides were then washed with PBS-Tx (0.3%) and then incubated with the appropriate secondary antibody (anti-mouse Pacific Blue or ant-rabbit Alexa 488; Thermo Fisher Scientific, Abingdon, United Kingdom) or a Pacific bluestreptavidin (Thermo Fisher Scientific) at a concentration of 1:500 for 3 to 4 hours at room temperature. Slides were then washed, mounted using Vectashield, and imaged using a confocal microscope (Zeiss, Cambridge, United Kingdom).

Novel long noncoding RNAs expressed in the rodent dorsal root ganglion
To identify unknown transcribed regions outside known gene models (ie, representation of RNA transcripts produced by a gene) and then group them into transcriptional units representing putative novel LncRNAs on the gene level, we performed RNA-sequencing (RNA-Seq) of DRG tissue harvested from animal models of nerve injury vs sham controls. We profiled novel LncRNAs alongside known annotated protein-coding genes and LncRNAs from the ENSEMBL genome database. 79 To obtain computational predictions of novel LncRNAs, we used a reference-based customised pipeline 45 that relies on a reference genome and gene set annotations. It uses the output of the STAR 18 aligner, and the quality of predictions is dependent on the aligner accuracy, the quality of the reference genome, and the completeness of the gene set annotation. A coverage threshold is used 7 to identify nonannotated continuously transcribed regions, 22 ie, I.o.E., and then clustered together and trimmed these regions using de novo identified SJs from mapping the RNA-seq reads to the organism's genome. We also applied a signal processing, smoothed z-score thresholding algorithm to further identify coverage drops (putative introns and transcription ends) and peaks (putative exons). To identify nonannotated I.o.E and to differentiate them from UTRs or not-yet-annotated exons belonging to known gene models, we filtered out reads overlapping ENSEMBL and RefSeq annotations as well as genomic regions of 1000 bp from the 59 and 39 end of known gene models. Doing so, we have excluded any predictions overlapping with a region of 1000 bp from both ends of annotated gene models. The mean length of 39 UTRs is 424 bp and 524 bp in rodents and humans, respectively. For the 59 UTRs, mean lengths are 127 bp in rodents and 146 bp in humans. 55 We should note that we used both major annotation consortia, so an incomplete gene annotation (ie, missing exons, UTRs, and isoforms) in just one of the annotated gene sets would not influence results. Given the fact that we discarded all predictions overlapping annotated genes, we could not also identify extracoding RNAs, as these overlap protein-coding genes on the sense strand. 62 32.6% of read pairs on average were overlapping regions outside gene models (as we used paired-end data, the units of evidence for gene expression are always pairs of reads). This finding is consistent with Ref. 22.
To acquire a complete representation of the nonannotated transcribed regions, we intentionally applied a low coverage threshold (.2 sequence coverage for the region) for at least the length of an RNA-seq read ($100 bp). The cutoff threshold is similar to the one applied in Ref. 22. Then, we clustered and trimmed these regions using splicing information from novel SJs identified by at least 2 RNA-seq read-pairs and a smoothed zscore over a rolling window to identify coverage drops. To predict and analyse LncRNAs on the gene level, we merged together all identified transcripts from individual samples to create a unified set of nonredundant, novel, putative LncRNAs in the form of a GTF file. The bioinformatic workflow is illustrated in Figure 1; for more details, see methods.
We then calculated the coding potential and performed transcriptional profiling to identify DE novel and known gene models between animal models of peripheral neuropathy and control (sham surgery) samples. To perform DE analysis without overestimating fold changes for lowly expressed transcripts, we used the analysis software DESeq2 42 and filtered DE results according to expression levels and consistency. Significant cutoff for DE was adjusted P value ,0.05 in all cases. We filtered out novel LncRNAs that were not expressed in at least all replicates of a condition or strain or were below an average expression threshold of .0.5 fpkm in at least one condition. We particularly focused on antisense LncRNAs, ie, overlapping exonic regions of the gene on the opposite strand, and intergenic LncRNAs, lying on the intergenic space between known gene models. All expression data including antisense and intergenic LncRNAs are available in http://www.painnetworks.org, 54 and all RNA-seq data (raw data and the whole gene set, ie, novel LncRNAs and annotated genes) are available in Gene Expression Omnibus (GEO) (GSE107182, GSE107180, and GSE107181). Supplemental spreadsheets of the complete data set are available at http://doi.org/10.6084/m9.figshare.6508205.
LncRNAs are known to have relatively poor conservation between species, so we aimed to compare one rat (Wistar) and 2 mouse strains (BALB/c and B10.D2 strains). RNA-seq quality was good for both experiments (supplementary tables 1-4, available at http://links.lww.com/PAIN/A676). Quality was assessed based on: the percentage of uniquely mapped reads (89% for mice and 72.51% for rats), the number of properly and concordantly paired reads, the average Phred score for read quality (32.2 for mice and 34.3 for rats), the base calling at the extremities of reads (0.08 for mice and 0 for rats, low Phred scores at the end of reads), the median of the insert between paired read mates (192.1 for mice and 153.4 for rats), and the GC content of reads for all sequencing lanes and samples (48.2% for mice and 51.4% for rats). Two mouse samples were excluded because they had much lower mapped reads (sample 59: 73.3%), a very high percentage of unmapped reads because they were too short (sample 59: 18.1%), much more reads mapped to multiple loci (sample 66 and sample 59), and higher Cook's distance (supplementary figure 1, available at http://links.lww.com/PAIN/A676). Uniquely mapped read pairs were used for the downstream analysis. Raw and processed gene expression data are available in GEO.
In total, we had on average 64 million uniquely aligned pairs of reads per sample for mice and 41 million pairs of reads for rats, more than enough both for identification of de novo LncRNAs and then to ask whether they were DE in injury vs sham conditions. 19 Using this approach, we initially identified and reconstructed thousands of nonannotated transcribed loci in the mouse and rat DRG with length of more than 200 bp. Then, we evaluated whether these novel transcripts were protein coding or noncoding. After coding potential calculation using the coding potential calculator CPAT, 72 we obtained 6657 long consistently expressed transcripts classified as noncoding in the rat DRG and 4729 in the mouse DRG. Four thousand three hundred fifteen of these passed the expression threshold in rats and 2693 in mice and were retained for the downstream analysis. These long novel transcripts with no coding potential were considered putative novel LncRNAs. The majority of novel LncRNAs were intergenic in both species, with 21% and 13% antisense (overlapping exonic regions of protein-coding genes on the opposite strand) in the mouse and rat DRG, respectively ( Figs. 2A and B). Most novel LncRNAs were multiexonic, with a distribution of exons heavily skewed towards biexonic transcripts (Figs. 2C and D). This exon distribution is very similar to GENCODE findings. 27 The usage of ribodepleted RNA for the library preparation allowed us to completely sample and interrogate the noncoding transcriptome and led to a relatively high proportion of intronic noncoding transcripts being identified. 29 To increase confidence and to gather more evidence regarding the completeness and expression of predicted novel LncRNAs, we examined the relationship between their genomic loci and annotated TSS. To do this, we used 59 CAGE (cap analysis gene expression) TSS data 40 that are available in the mouse. CAGE is a technique for high-throughput identification of sequence tags corresponding to 59 ends of RNA at the cap sites and the identification of the TSS. 65 As TSS data were not available for the mouse reference genome mm10, we translated and mapped coordinates from the mm9 genome to mm10. We then calculated the kernel density of the distance between the TSS that were mapped to mouse genome mm10 (51% of the TSS available for February 2019 · Volume 160 · Number 2 www.painjournalonline.com 469 mm9) and a region within 100 bp of the 59 end of the putative LncRNA similarly to Ref. 27. For both known and novel LncRNAs, the kernel density was highest at 0 but with more spread for the novel LncRNAs (Fig. 2E, supplementary Figure 2; available at http://links.lww.com/PAIN/A676). 33.7% of the antisense and intergenic LncRNAs had a predicted TSS on their 59 end on the same strand of the genome. When we removed the outlying distant TSS-ie, more than 1.5 times the interquantile range of the distribution-the mean distance of TSS and novel LncRNAs was 119 bp upstream of the predicted transcript. In a GENCODE study, 27 15% of identified LncRNAs in humans had a TSS on their 59 end. These results highlight how close the 59 end of novel LncRNAs was to experimentally determined TSS. Due to the fact that only a fraction of TSS were mapped to the current mouse genome, these data are likely to be an underestimate. These findings are in concordance with those of Ref. 17. We also note that a significant fraction of novel LncRNAs are either incomplete, not yet annotated extended UTRs or that LncRNAs are indeed  In both species, we found that LncRNAs were significantly and consistently expressed at a lower level than protein-coding genes (Fig. 3A). This ratio of about 10-fold lower expression of LncRNAs to protein-coding genes is similar to previous studies. 69,76 Nociceptors are a major component of the DRG, and a number of pain genes are selectively expressed by these neurons. We therefore studied the genomic context of identified LncRNAs, and in particular, whether they were antisense or in close genomic proximity with known pain genes downloaded from the Pain Genes Database. 35 The Pain Genes Database catalogues genes that have been shown to have an impact on pain-related behaviour in rodents based on transgenic knockout experiments. Of the 449 genes in the database, we have found 13 novel LncRNAs antisense to pain genes in mice and 19 in rats (supplementary tables 5 and 6, available at http://links.lww.com/ PAIN/A676). Twenty-three intergenic LncRNAs had a pain gene as their closest genomic neighbour in mice and 57 in rats (supplementary tables 7 and 8 available at http://links.lww.com/ PAIN/A676). Ion channels are key determinants of sensory neuron excitability. We identified LncRNAs antisense to voltagegated sodium channels, potassium channels, calcium channels, chloride channels, and TRP channels, within the mouse and rat DRG, as shown in supplementary tables 9 and 10 (available at http://links.lww.com/PAIN/A676).
In general, we had modest syntenic conservation (ie, in equivalent genomic positions) between species. We used synteny portal 37 to retrieve synteny blocks conserved between the human (GRCh38.88), mouse (mm10), and rat (rn5) with a resolution of 150,000 bp. We lifted genomic coordinates from rn5 to rn6 genome and found in total 912 conserved synteny blocks in humans and mice and 443 uniquely mapped in the current rat genome. Eight hundred (18.5%) novel LncRNAs in rats and 1271 (47%) in mice were in these syntenic blocks conserved between the 3 species. Moreover, 649 LncRNAs in mice and 782 in rats were in 200 common syntenic blocks between the 2 species (supplementary data available at http://doi.org/10.6084/ m9.figshare.6508205, S. Data 1). Five hundred nine LncRNAs in rats and 397 in mice were antisense of the same orthologous genes in mice and rats (S. Data 2 and 3, available at http:// doi.org/10.6084/m9.figshare.6508205).

Different types of dorsal root ganglion neurons selectively express LncRNAs
Dorsal root ganglion neurons are heterogenous both in terms of physiology and molecular profile; to identify whether the expression of LncRNA may be DRG subtype-specific, we reanalysed a published data set of single-cell RNA sequencing data 38 derived from C57BL/6 mouse DRG neurons for expression of the novel LncRNAs that we had identified. The authors had previously generated 10 different DRG neuron subtypes from their analysis of 197 neurons. We found that we could effectively identify transcriptome-and size-based neuronal types based on the selective expression of both ENSEMBL annotated genes and our novel LncRNAs.
Principal component analyses of the expression set of novel LncRNAs showed that samples belonging to most of the sensory neuron subtypes were clustered together (Fig. 3B) but with a higher spread than for annotated genes (supplementary figure 3A, available at http://links.lww.com/PAIN/A676), suggesting that these can be important transcriptional units the expression of which is highly related to the different subtypes of DRG neurons.
We then used the Tau tissue specificity metric 77 to identify highly DRG subtype-specific LncRNAs expressed in the mouse DRG. The density plot of the Tau index distribution (supplementary figure 3B, available at http://links.lww.com/PAIN/A676; and S. Data 4, available at http://doi.org/10.6084/m9.figshare.6508205) suggested that the Tau .0.8 was an appropriate threshold for separating neuron subtype-specific from ubiquitous annotated genes and novel LncRNAs. To consider a gene or LncRNA neuron subtype-specific, we also imposed a cutoff mean expression threshold of .3 log2-normalised read pairs across all samples of a neuronal subtype. Eighty-seven novel LncRNAs, 119 annotated LncRNAs, and 597 protein-coding genes were neuron subtype-specific. This reveals a significant (Fisher exact test P value ,0.001) enrichment of novel and annotated LncRNAs among the neuron subtype-specific transcripts (Fig. 3C). This quantification of the DRG subtype specificity confirms that LncRNAs' expression pattern is more tissue specific than protein-coding genes'. In studying the different DRG subtypes (as defined by Li et al., 2016), we noted that there was an uneven distribution of subtype-specific LncRNAs, Figure 3D.

IPSC-derived sensory neurons differentially express known and novel LncRNAs compared with their respective IPS cells
We assessed LncRNA expression in human IPSC-derived sensory neurons generated from healthy individuals. We followed a differentiation protocol known to produce neurons with a gene expression profile and functional characteristics that are very similar to rodent sensory neurons. 8,11,80 Virtually, all the resulting neurons express the sensory neuron marker Brn3a, project extensively arborized neurites, a subset can be myelinated and exhibit mature electrophysiological characteristics of sensory afferents. 8,73 At least 84% of the ion channel genes known to be expressed in the human DRG were shown to be expressed by sensory neurons generated using this protocol. 80 We compared gene expression in IPSC-derived sensory neurons with expression in IPSC parent lines to focus on LncRNAs enriched in differentiated sensory neurons (GEO GSE107181, S. data 5-7, available at http://doi.org/10.6084/m9.figshare.6508205). We first interrogated genes annotated in the ENSEMBL GRCh38.88 gene set, and we also identified and profiled 2948 novel LncRNAs in human IPSC and IPSC-derived sensory neurons, most of them were intergenic (Fig. 4A) and biexonic (Fig. 4B). Again, we found that LncRNAs were significantly lower expressed than proteincoding RNAs (Fig. 4C). IPSC-derived sensory neurons from 3 different cell lines (identified as AD2, AD4, and NHDF) were very similar to each other and considerably different to the IPSC parent lines (supplementary figure 4A, available at http://links.lww.com/ PAIN/A676). Forty-two percent of all expressed genes were DE in IPSC-derived sensory neurons vs IPSCs, which is a remarkable transcriptional change. A total of 6103 annotated LncRNAs (ENSEMBL GRCh38.88) were consistently expressed in at least all samples of either IPSCs or IPSC-derived sensory neurons. A total of 1830 (30%) of them were significantly DE in all 3 cell lines between IPSC-derived sensory neurons and IPSCs; the majority of these were intergenic (47%). Seventy-seven percent of the expressed novel LncRNAs were significantly DE between IPSCderived sensory neurons and IPS cells. All 3 cell lines had in common 371 novel LncRNAs significantly DE.
In both annotated and novel intergenic LncRNAs, distance was modestly but significantly anticorrelated with expression correlation with their closest genomic neighbour (supplementary figure 4B, C; February 2019 · Volume 160 · Number 2 www.painjournalonline.com  PAIN ® available at http://links.lww.com/PAIN/A676). Two thousand seven hundred forty-three intergenic LncRNAs (novel and annotated) were significantly DE between IPSC-derived neurons and IPSC (Fig. 4D). Eighty annotated and 16 novel LncRNAs were antisense to pain genes (S. Data 8, available at http://doi.org/10.6084/m9.figshare.6508205), and 15 annotated and 18 novel LncRNAs had a pain gene as its closest genomic neighbour (S. Data 9, available at http://doi.org/ 10.6084/m9.figshare.6508205). Eight of the LncRNAs antisense to pain genes were significantly DE between IPSCderived sensory neurons and IPSCs with a significantly DE gene on the opposite strand and anticorrelated expression changes, 1 of them was novel (supplementary table 11, available at http:// links.lww.com/PAIN/A676). IPSC-derived sensory neurons expressed LncRNAs (including both annotated and novel LncRNAs) antisense to TRP, voltage-gated sodium, potassium, chloride, and calcium ion channels; in some cases, these showed the opposite expression changes relative to the sense gene on the opposite strand (Figs. 4E and F).
One example of the antisense LncRNAs, the HOXD Cluster Antisense RNA 1 (HAGLR) was further investigated because it was very highly upregulated in neurons vs IPSC, and HoxD genes have been implicated in nociceptor specification. 25 HAGLR was ranked first among the significantly DE annotated antisense LncRNAs by its log2 fold change (8.74) and base mean counts (261.5). Differential expression was validated by qPCR (Fig. 5,  supplementary table 12, available at http://links.lww.com/PAIN/ A676). In situ hybridization also revealed that it was expressed by mouse DRG cells in vivo, and it was found to be significantly downregulated in the mouse DRG after nerve injury (see below as well as Fig. 5).
We also observed similar extent of syntenic conservation between the human, mouse, and rat in novel LncRNAs identified in human IPSC and IPSC-derived sensory neurons. One  Human DRG eQTLs (expression quantitative trait loci) have recently been identified and associated with pain phenotypes. 52 We interrogated this data set for overlaps with novel and annotated LncRNAs expressed in IPSC-derived sensory neurons, as this may provide an underlying mechanism through which a single nucleotide polymorphism (SNP) at this site may impact on gene expression. We only considered overlaps valid if they were found on exons of the respective transcript. We found that 4 annotated LncRNAs (3 antisense and 1 intergenic) that had DRG eQTLS were expressed in IPSC-derived sensory neurons, and also that 5 novel LncRNAs were overlapping DRG eQTLS ( Table 1).

Expression profiling of LncRNAs after nerve injury
We then assessed differential expression of LncRNAs within the DRG after nerve injury both in rats and mice. We used the SNI model in mice and the SNT model in rats. We used 2 different mouse strains: BALB/c and B10.D2, which have previously been shown to develop different levels of mechanical hypersensitivity after nerve injury. 67

Transcriptional changes of LncRNAs in the rat dorsal root ganglion
We confirmed that the expression patterns of the top 5000 ENSEMBL annotated genes and novel LncRNAs ranked by their SD could separate samples according to condition (Fig. 6A). We also observed a significant transcriptional response after nerve injury in the rat DRG, which amounts to 25.5% (4215 ENSEMBL annotated genes 1 novel LncRNAs) of all expressed genes (  (supplementary table  13, available at http://links.lww.com/PAIN/A676). There were novel LncRNAs antisense of pain genes, voltage-gated potassium and sodium channels (Fig. 6C).
Six hundred fifty-four intergenic LncRNAs were found to be significantly DE in SNT vs sham, 575 of them were novel (Fig. 6D share.6508205). Five of them (4 novel and 1 annotated) were adjacent to and highly correlated with pain genes (supplementary Table 1 Annotated and novel LncRNAs overlapped by eQTLs.

Transcriptional changes of LncRNAs in the mouse dorsal root ganglion
Withdrawal thresholds to von Frey filament stimulation confirmed previous findings that the B10.D2 demonstrates less mechanical pain-related hypersensitivity after SNI than BALB/c 67 (Fig. 7A). There were no significant behavioural differences between male and female mice for either strain. Principal components analysis of the expression of the top 10,000 ENSEMBL annotated genes and novel LncRNAs (ranked by SD) showed that they were able to optimally separate mouse samples according to sex, strain, and within them, condition (Fig. 7B).
Volcano plots show the range of transcriptional changes of ENSEMBL annotated genes and novel LncRNAs in mice after peripheral nerve injury (Figs. 7C and D, S. Data 15, available at http://doi.org/10.6084/m9.figshare.6508205). Novel LncRNAs were a substantial component of the DE genes following the SNI pain model (highlighted in Figs. 7C and D).
Regarding the BALB/c mouse strain, which showed significantly more pain-related hypersensitivity from day 7 onwards than the B10.D2 strain, we found significantly more DE ENSEMBL annotated genes and novel LncRNAs, 2750 compared with 1441. In comparison to rats, we found less DE genes in the mouse after nerve injury. Spinal nerve transection is a more severe injury model than SNI (in which less sensory neurons are axotomised and the injury is more distal), and so, this may reflect model severity rather than species differences. In those ENSEMBL annotated genes that were DE, we found significant GO enrichment in terms related to the nervous system, regulation of excitability, signal transduction, and response to injury (supplementary figure 7, available at http://links.lww.com/ PAIN/A676). As novel LncRNAs were not assigned with GO terms, we used their expression profiling context, in an unbiased way, to gather more insights regarding their possible functional importance. Under the assumption that they may regulate gene expression either in-cis or in-trans, we further studied them in the context of modules of closely coexpressed genes and associated them with enriched GO BP terms. We first created a weighted gene coexpression network (WGCNA) 36 and then identified modules of highly coexpressed ENSEMBL annotated genes and novel LncRNAs (scale independence and dendrogram of merged/unmerged modules in supplementary figure 8, available at http://links.lww.com/PAIN/A676). The absolute weighted bicorrelation across all samples (n 5 20) was used to construct the network. An n . 15 is needed to robustly calculate expression correlations. We also calculated the representative gene (ie, eigengene) for each module, defined as the first principal component of the expression of all member genes. We then quantified the novel LncRNAs' module membership by calculating the robust weighted bicorrelation of the LncRNAs with these eigengenes. Next, we performed a GO enrichment analysis and annotated each module with its top enriched BP term, based on an overrepresentation analysis of the annotated genes. The strength of module membership for novel LncRNAs (Fig. 8A) shows highly correlated modules of LncRNAs associated with Figure 7. Dorsal root ganglion transcriptional response in mice after peripheral neuropathy. (A) Hind paw withdrawal thresholds to von Frey filament stimulation 1 SEM in grams. We calculated the area over the curve (AOC) for each strain and obtained the percentage of maximum induced hypersensitivity. Two-way ANOVA showed a significant effect of surgery and a significant interaction of strain:surgery (P 5 0.001). One-way ANOVA between SNI groups on D28 showed significant difference in % of maximum hypersensitivity (P 5 0.002). Black bar shows comparison between SNI groups, P , 0.01**.    RNA-processing and some related to the nervous system, signalling, and regeneration. The vast majority of novel LncRNAs were in the module associated with RNA-processing (Fig. 8B).
Of the 7990 annotated ENSEMBL (GRCm38.87) LncRNAs, 2406 were expressed in the mouse DRG. A total of 296 LncRNAs were significantly DE in BALB/c strain SNI vs sham, 193 of them were novel. In B10.D2 strain, 146 LncRNAs, 97 of which were novel, were significantly DE (PainNetworks and GEO GSE107180). Although the absolute numbers are much smaller in comparison with annotated genes, percentages of significantly dysregulated transcripts are similar. Most of them were intergenic (S. Data  ). The majority of LncRNAs antisense to potassium channels were downregulated after nerve injury and between IPSC-derived neurons and IPSC. On the other hand, although the LncRNA antisense to Scn9a was downregulated after nerve injury, all other LncRNAs antisense to sodium channels were upregulated in IPSC-derived neurons vs IPSC. In some cases, these showed an opposite expression pattern to their sense gene.
In total, 2365 intergenic LncRNAs were expressed in the mouse DRG, 1282 of them were novel, and 126 were significantly DE in mouse SNI vs Sham (Fig. 9E). Similarly to rats and humans, intergenic LncRNAs showed positive correlation with their adjacent gene (supplementary figure 9, available at http:// links.lww.com/PAIN/A676). Twenty-three novel and 18 known intergenic LncRNAs were adjacent to pain genes in mice. Two of these novel intergenic LncRNAs were significantly DE and highly correlated with their adjacent pain gene (supplementary table 17, available at http://links.lww.com/PAIN/A676). In some cases, multiple novel LncRNAs related to coding genes in the same pain-related signalling pathway could be identified and were DE in mice and rats. For instance, 2 intergenic LncRNAs upstream of the opioid receptor genes, Oprl1 and Oprd1, were significantly DE and highly correlated with their adjacent significantly DE gene in both species. Oprl1 and Oprd1 form heterodimers and appear in the same network of highly correlated genes. 54 We also performed unsupervised clustering of samples based on the expression of all novel and annotated LncRNAs. We observed a pattern of highly sex-specific expression changes within strains and then, within sex, we had separation according to condition (Fig. 9F). This indicated that LncRNAs are dysregulated after peripheral nerve injury in a sex-and straindependent manner.
We selected and validated 7 representative novel LncRNAs based on the expression strength, the significance and size of the effect (log2 fold change), and the genomic context. These novel LncRNAs were all found to be significantly DE in the mouse DRG after nerve injury and among these are antisense and intergenic, Quantitative real-time polymerase chain reaction confirmed the expression of these novel LncRNAs in the mouse DRG and validated their dysregulation after nerve injury, Figure 10. With the exception of 1 LncRNA, these all demonstrated higher expression in the DRG compared with brain (Fig. 10E). This comprehensive study 47 showed that upstream of genes where these LncRNAs lie, there were no previously unannotated lengthened 39 UTRs. With the exception of LncRNA4714, upstream of Oprd1, where the multiexonic transcript we have identified is much longer than any predicted elongated UTR. When studying the changes in LncRNA expression evoked by SNI in B10.D2 and BALB/c mouse strains, there was in general good agreement between RNA-seq and qPCR findings. We also found that in all cases, where RNA-sequencing showed significant dysregulation, qPCR confirmed the result (Fig. 10).

Discussion
We used high-depth RNA sequencing and a dedicated bioinformatic pipeline to identify thousands of known and putative novel LncRNAs expressed in the mouse and rat DRG and human IPSC-derived sensory neurons. Many of these LncRNAs were February 2019 · Volume 160 · Number 2 www.painjournalonline.com  antisense or adjacent to known pain and ion-channel genes. A significant proportion demonstrated selective expression in DRG neuron subtypes. Novel LncRNAs were DE after peripheral nerve injury in a species-and strain-dependent manner, including novel antisense LncRNAs with opposite transcriptional changes to their sense genes and intergenic LncRNAs highly correlated with their adjacent gene. Thus, LncRNAs that have been a relatively unexplored part of the DRG transcriptome demonstrate remarkable complexity in terms of their relationship to genes known to impact on sensory function, DRG cell-type-specific expression, and the transcriptional response to nerve injury. LncRNAs have been shown to be: smaller than protein-coding genes, spliced, biexonic, transcribed as independent transcription units, expressed at a lower level than protein-coding genes, and show highly spatially constrained correlation with their antisense and adjacent genes. 17,69,81 These characteristics are highly consistent with the antisense and intergenic LncRNAs that we identified in this study as being expressed in the DRG.
No previous attempts have comprehensively determined LncRNA expression within the DRG. Previous studies have undertaken a candidate gene approach, 32,86 whereas others have used a high-throughput approach such as RNA-seq but have only interrogated the expression of previously annotated LncRNAs. 24,30,75 As we have shown, this strategy will miss thousands of unannotated and potentially important LncRNAs. In fact, 8% of genes within the Pain Genes Database in the mouse, 16% in the rats, and 7% in human IPSC were found to have an antisense or neighbouring intergenic LncRNA. Furthermore, LncRNAs antisense to ion channels were not limited to KCNA2 and SCN9A, but could also be identified to other sodium, potassium, calcium, and also TRP channels, all of which have key roles in modulating sensory neuron function. As with all computational predictions, this study has the limitation of the possible inclusion of false positives in our list of novel LncRNAs; however, we have validated the expression of a number of novel LncRNAs using the independent technique of qPCR.
It is known that LncRNAs can be highly tissue-and celltype-specific. 57,68 Dorsal root ganglion cells are heterogeneous in terms of their physiology, anatomy, and connectivity; recent single-cell analysis shows how gene expression relates to such specialised features but also reveals even greater complexity in DRG subtypes based on molecular profiling. 38,70 We showed that both novel and annotated LncRNAs demonstrate selective expression within sensory neuron subtypes. As such, identification of LncRNAs may help in the molecular classification of DRG cells. Moreover, both novel and annotated LncRNAs were significantly enriched in neuron-subtype-specific genes vs protein-coding genes. The functional specialisation of different neuron subtypes will principally arise due to the neuron subtypedependent expression of protein-coding genes; however, such expression may be sculpted by LncRNAs.
We have also defined LncRNAs expressed in human IPSCderived sensory neurons. Advantages of studying LncRNA expression in this model (rather than, eg, cadaveric human DRG) is the ready source of high-quality RNA from a pure sensory neuronal population (with very few contaminating Schwann cells) and the ability to compare expression to the parent IPSC line. Disadvantages include the fact that although subpopulations of sensory neurons exist in these cultures, these are not as diverse as native DRG cells and are less mature given the lack of target interactions.
We identified almost 2000 previously annotated and over 2000 novel putative LncRNAs, which were DE when comparing IPSCderived sensory neurons and parent IPSC lines. These demonstrated many characteristics of the LncRNAs expressed in the rodent DRG. Twenty-nine percent of genes in the Pain Genes Database had a potentially relevant relationship to a novel or annotated LncRNA either antisense LncRNAs on the opposite strand or adjacent to an intergenic LncRNA. Interestingly, some of these LncRNAs may have a role in shaping the expression of the sense pain gene during differentiation, as both novel and annotated, DE antisense LncRNAs demonstrated anticorrelated expression changes to their DE sense pain gene. We also identified human LncRNAs that were antisense to voltage-gated ion channel genes. It has been shown that many LncRNAs overlap genomewide association studies' traits, and LncRNAs overlapping traitassociated SNPs are specifically expressed in cell types relevant to the traits. 28 The identification of human LncRNAs expressed in IPSC-derived sensory neurons may therefore enable investigation of the genetic basis of chronic pain states in humans, especially under conditions of neuropathic pain wherein maladaptive plasticity in the DRG is an important pathophysiological driver. A recent study has mapped eQTLs in the human DRG, 52 and a number of SNPs that impacted on gene expression were found to be within annotated LncRNAs. In our data, 9 LncRNAs expressed in IPSCderived sensory neurons were overlapped by eQTLs. Three of the novel and 3 of the known were also DE in IPSC-derived neurons vs IPSC. Dorsal root ganglion eQTLs identified by Ref. 52 were found among hits in numerous genome-wide association studies. Such studies highlight the utility in identifying LncRNAs, in this case by providing explanatory power as to how an SNP linked to complex disorders may impact on gene expression and phenotype. As more tissue becomes available for RNA-seq, it will also be possible to extend LncRNA discovery to the postmortem human DRG. 60 Nerve injury is known to elicit remarkable transcriptional changes within the DRG, which has deleterious consequences such as neuropathic pain, but in certain contexts can also be adaptive; eg, by promoting nerve repair. 51 We investigated the expression of LncRNAs in Wistar rats after SNT and in 2 different mouse strains after SNI. The BALB/c strain was previously shown to develop greater levels of mechanical hypersensitivity after SNI vs the B10.D2 strain. 67 Hundreds of LncRNAs were DE in both rats and mice after nerve injury. We found that more proteincoding genes, annotated LncRNAs, and novel LncRNAs were DE in the high pain BALB/c strain compared with the low pain B10.D2 strain. Most of the variance of the LncRNAs' expression in our data set was between strains and not conditions; namely, LncRNAs demonstrated strain-dependent expression plasticity after nerve injury. Twenty-seven percent of LncRNAs were significantly DE between mouse strains compared with 12.2% for ENSEMBL protein-coding genes. Considering the rapidly increasing literature on sex differences in pain processing 48,66 that there were sex-dependent effects on LncRNA expression.
Weighted gene coexpression network analysis revealed that novel LncRNAs in the mouse DRG were in network modules related to RNA-processing and some of them in modules related to myelination, development, and regeneration of the nervous system, immune response, and signalling. This could indicate that these LncRNAs function together with genes acting as transcriptional regulators associated with posttranscriptional modifications. This enrichment is different from that observed in the ENSEMBL annotated gene set, where terms related to axon guidance, ion channel transport, regulation of synapse assembly and of neuron apoptotic process, nervous system development, and memory were significantly enriched. These findings are similar to known biological processes enriched after nerve injury 83,84 and in pain genes. 41 Enrichment in biological processes of splicing, mRNA processing, and polyadenylation (ie, parent, child, and related terms to RNA-processing) has however been reported. 82,83 Our finding that the expression of novel LncRNAs changes together with annotated genes regulating RNA-processing is consistent with known mechanisms by which LncRNAs alter gene expression (discussed below).
LncRNAs can regulate expression in cis and in trans, 33 but the genomic context of LncRNAs is important for both antisense transcripts regulating the gene on the sense strand or intergenic LncRNAs. Intergenic and antisense LncRNAs, which tend to lie in genomic regions populated with genes, 31 have correlated expression patterns with their adjacent genes and may regulate gene expression (Fig. 4). 58 We found that the shorter the distance, the stronger the correlation in all species, both for known and annotated LncRNAs. However, the network analysis we performed allowed us to identify in an unbiased way modules of highly correlated genes and LncRNAs across all samples. Thus, these LncRNAs are putative both in-cis and in-trans regulators. Regarding the relation of LncRNAs with their genomic context, we found that more than 45% of antisense LncRNAs had anticorrelated-opposite expression changes in respect to the sense gene, whereas only 10% of intergenic LncRNAs had negative correlation with the expression of their closest genomic neighbour. Twelve pairs of antisense LncRNA/sense gene with opposite LFCs reached significance between pain models and control animals. These antisense LncRNAs fit into the paradigm of Kcna2 antisense, 86 which silences the gene on the opposite strand. LncRNAs are also known for regulating clusters of imprinted genes or close genes such as the Hoxd cluster. 43,78 One of these LncRNAs, HAGLR, was the most upregulated and stronger expressed LncRNA in human IPSCderived neurons and was also found by in situ hybridization to be expressed in the mouse DRG and significantly downregulated in both mouse strains after nerve injury (supplementary table 12, available at http://links.lww.com/PAIN/A676).
We used qPCR to validate the expression of a number of novel LncRNAs in close genomic relationship to protein-coding genes of relevance to sensory function: A novel LncRNA antisense of Nefl gene implicated in the Charcot-Marie-Tooth disease, 85 and an intergenic LncRNA close to and highly correlated with the pain gene Oprd1 was found DE and validated with qPCR. We also describe a novel intergenic LncRNA in close proximity to the Lrrc4 gene that relates to axon guidance and synapse organisation, 16,74 and finally, a significantly DE intergenic LncRNA was validated close to the voltage-gated sodium channel subunit Scn4b. The majority of LncRNAs validated by qPCR were found to be more highly expressed in the DRG than brain. LncRNAs are putative therapeutic agents that could regulate the expression of target genes related to disease. 71 Although application of such therapeutics to pain would need to overcome the hurdle of delivering therapeutics to DRG cell bodies.
In summary, we have provided a resource, in which we have defined LncRNA expression within the DRG across species. We show that marked changes in LncRNA expression occur after nerve injury and during sensory neuron differentiation. LncRNAs' expression is DRG subtype-specific, and there is often highly spatially constrained correlation/anticorrelation with their antisense and adjacent genes.