RNA Sequencing Reveals Novel Transcripts from Sympathetic Stellate Ganglia During Cardiac Sympathetic Hyperactivity

Cardiovascular disease is the most prevalent age-related illness worldwide, causing approximately 15 million deaths every year. Hypertension is central in determining cardiovascular risk and is a strong predictive indicator of morbidity and mortality; however, there remains an unmet clinical need for disease-modifying and prophylactic interventions. Enhanced sympathetic activity is a well-established contributor to the pathophysiology of hypertension, however the cellular and molecular changes that increase sympathetic neurotransmission are not known. The aim of this study was to identify key changes in the transcriptome in normotensive and spontaneously hypertensive rats. We validated 15 of our top-scoring genes using qRT-PCR, and network and enrichment analyses suggest that glutamatergic signalling plays a key role in modulating Ca2+ balance within these ganglia. Additionally, phosphodiesterase activity was found to be altered in stellates obtained from the hypertensive rat, suggesting that impaired cyclic nucleotide signalling may contribute to disturbed Ca2+ homeostasis and sympathetic hyperactivity in hypertension. We have also confirmed the presence of these transcripts in human donor stellate samples, suggesting that key genes coupled to neurotransmission are conserved. The data described here may provide novel targets for future interventions aimed at treating sympathetic hyperactivity associated with cardiovascular disease and other dysautonomias.

Here, we carried out an RNA sequencing transcriptomic analysis to identify differential gene expression in sympathetic neurons from normotensive and hypertensive rats. Importantly, our RNAseq results demonstrated significant differential expression in transcripts that regulate [Ca 2+ ] i signalling in addition to genes that regulate cN signalling in disease. We validated our findings using quantitative real-time polymerase chain reaction (qRT-PCR) and subsequently assessed whether these genes were present in sympathetic stellate ganglia from human donors, in order to provide a framework for future validation and targeting.

Results
Study Design. To study the rat stellate ganglia transcriptome, we obtained RNA from four male 16-week-old SHR with hypertension and four male age-matched Wistar control rats. RNA extracts from these ganglia contained ~8.5 ng/μL RNA. RNA integrity (RIN) was validated using an Agilent 2100 Bioanalyzer. To ensure good quality RNA we only used samples with RINs >8.5 for sequencing. Using the HiSeq 4000 sequencing platform, we generated ~345 million reads of 75-bp paired-end RNAseq data, i.e. ~22 million reads per sample. We performed a comparative analysis of the stellate ganglia transcriptome between Wistar and SHR phenotypes. The number of animals per group (n = 4) and the sequencing parameters established were based on the recommendations of the Wellcome Trust Centre for Human Genomics, in addition to those published by Conesa et al. 18 . Figure 1 illustrates our RNAseq pipeline.
Quasi-mapping and Differential Expression Analysis. Transcript level quantification was performed using the quasi-mapping function of the Salmon package, where inbuilt functions corrected for GC and positional bias 19 . These data were converted into quantitative gene level expression values, and the differential gene expression between the strains was calculated and normalised using DESeq2 20 . We obtained an average number of mapped reads (74.75%) that was not significantly different between groups (Table 1). DESeq2 was used to generate an MA plot for visual representation of our genomic data ( Supplementary Fig. 1a) and a Principal Component Analysis (PCA) plot with log normalisation was produced to demonstrate variation within our samples ( Supplementary Fig. 1b) 20 . As confirmation of the expected neuronal phenotype, we identified that markers consistent with sympathetic neurons were among the most highly expressed genes in both strains (e.g. Dopamine-β hydroxylase, Dbh; and Neuropeptide Y, Npy; see Table 2). The differential expression analysis identified 776 differentially expressed genes in the stellate ganglia, where approximately 55% of genes were upregulated and 45% downregulated in the SHR. The top 20 differentially expressed genes are displayed (Fig. 2).

Gene Ontology Analysis.
To assess the relevance of the observed differentially expressed genes, we grouped these genes into Gene Ontology (GO) families 21 using the ClusterProfiler R package 22 . We identified significant differences between the SHR and Wistar transcriptome in the following GO categories: Molecular Function (MF)   Tables 1, 2). Importantly, the differentially expressed transcripts within these GO groups may be linked to sympathetic dysfunction observed in the SHR model. There were no significantly different GO groups in the Biological Process (BP) category between Wistar and SHR.

Differential Expression of Extracellular Ligand-Gated Ion Channel Genes: Confirmation by
qRT-PCR. We obtained four independent Wistar and SHR RNA stellate ganglia to investigate whether we could replicate the RNAseq differential gene expression using qRT-PCR. In the 15 selected genes, we obtained similar trends using qRT-PCR as that observed in the RNAseq data; although due to variation within samples, only 4:15 genes reached significance (Fig. 4

Extracellular Ligand-Gated Ion Channel Genes Are Expressed in Human Sympathetic Stellate
Ganglia. To assess whether these genes are conserved in other species and have potential relevance for human sympathetic activity, we validated the expression of these genes in ganglia obtained from human donors. Given the limited availability of human tissue from patients, it was not possible to assess the role of these genes in human disease itself. Importantly, however, we detected the presence of glutamate (AMPA, NMDA and Kainate) and GABA receptor subunits which have not been previously identified within the human sympathetic stellate ganglia (Fig. 4). We also identified the presence of purinergic, glycinergic and serotonergic receptors.
Differential Expression of Phosphoric Ester Hydrolase Activity Genes. The 'phosphoric ester hydrolase activity' (GO:0042578) GO group was found to contain the highest number of differentially expressed  Table 3). Importantly, many of the same genes were represented in each of these groups, therefore we combined the KEGG signalling pathways (Circadian Entrainment, Dopaminergic Synapse and Glutamatergic Synapse) to identify the possible involvement of these transcriptomic changes on intracellular signalling (Fig. 6b,c). The KEGG enrichment pathway observations indicate that RNA transcripts associated with glutamate signalling are differentially expressed within SHR cardiac PGSNs. Within these KEGG groups, genes associated with glutamate signalling are coupled to [Ca 2+ ] i directly, where differences in subunit composition of glutamate-activated ion channels affect Ca 2+ permeability. Moreover, cyclic nucleotide generation and activation of their effectors are where green and blue represent up-and downregulated genes respectively (a). The presence of these genes was validated using qRT-PCR and the ΔΔCT comparative method was applied for analysis (SHR/Wistar); with β2 microglobulin (B2m) selected as the housekeeping gene for comparison as there was no difference in B2m expression between strains in the RNAseq dataset. Log 2 fold change of the differentially expressed genes in SHR stellates as identified by RNAseq are represented in red (b). Log 2 fold change ± SEM of qRT-PCR SHR data are also depicted (blue; b) relative to Wistar. Significance is indicated above the relevant bars. The RNAseq trends for up-or downregulation of each gene was confirmed with qRT-PCR; however, only Chrnb4 (p < 0.05), Gria1 (p < 0.001), Grik1 (p < 0.05) and Grin2b (p < 0.05) were significantly different to the expression levels measured in Wistar samples. The presence of these genes was confirmed in human stellate ganglia samples (c). To quantify mRNA expression in human stellate samples, data were normalised to the housekeeping gene B2m and the ΔCT method was applied for analysis.

Discussion
We have used RNAseq to identify alterations in the transcriptome of the stellate ganglion between SHR and Wistar rats, in order to establish lead genetic candidates that may be linked to sympathetic hyperactivity in this hypertensive model [26][27][28][29][30] .
Here we report several novel findings. First, we observed differential expression in the GO group encoding extracellular ligand-gated ion channels that are linked to intracellular Ca 2+ regulation. These changes were confirmed by qRT-PCR. Secondly, we provide evidence for the impaired expression for a family of genes that encode PDEs and phosphatase activity in hypertension. Altered transcription of PDEs and associated phosphatases that regulate PDE activity may be coupled to the dysregulation of cyclic nucleotide signalling, which is impaired in neurogenic hypertension 2 . Finally, to provide physiological context, we show that key differentially expressed genes are conserved and present in human stellate ganglia.
We observed high levels of transcripts involved in pathways governing adrenergic synthesis (e.g. Dbh, Th), neurotransmission (e.g. Snap25) and reuptake (e.g. Slc6a2), corroborating similar recently published profiles of stellate neurons [31][32][33] . We found that the highest number of differentially expressed genes between Wistar and SHR stellate ganglia within the molecular function (MF) gene ontology (GO) category, were associated with the 'extracellular ligand-gated ion channel' (GO:0005230) and the 'phosphoric ester hydrolase activity' (GO:0042578) superfamilies. In particular, our data highlighted significant differences in the mRNA transcript expression levels of several glutamate-activated ion channel subunits including: upregulation of the NMDA receptor (R) NR2 subunit (Grin2b) and AMPA GluA3 subunits (Gria3); as well as downregulation of AMPA GluA1 and GluA2 subunits (Gria1, 2) and Kainate 1 and 2 R subunits (Grik1, 2 respectively), supporting previous evidence for the role of amino acids in regulating cardiac-sympathetic activity [34][35][36][37] . Within the 'extracellular ligand-gated ion channel' GO Group, we also identified significant changes in the transcript expression levels of purinergic Rs (P2rx3, 4, 6) and GABA A R subunits (Gabra1, Gabra2). We also observed a significant downregulation of the ionotropic serotonin type 3 R subunits (Htr3a, Htr3b), nicotinic acetylcholine R β4 subunit (Chrnb4) and the chloride (Cl − )-permeable glycine channel β subunit (Glrb). Our RNAseq data provides an accurate prediction of the differential expression of these genes, as the trend for up-or downregulation was confirmed by qRT-PCR (Fig. 4). Moreover, using qRT-PCR we confirmed the presence of these genes within human male donor PGSNs  Fig. 4). Importantly, differential expression of ligand-gated ion channel subunits may have implications for ion permeability. The expression of these as-yet unreported ion channels on human sympathetic nerves has clear clinical relevance, as it may offer new targets for interventions aimed at modulating neurotransmission as well as considerations regarding the side effect profiles of current therapeutics.
We also analysed the relevance of differentially expressed genes by performing a functional enrichment KEGG analysis. This identified 'Circadian Entrainment' , 'Dopaminergic Synapse' , 'Retrograde Endocannabinoid Signalling' , 'Glutamatergic Synapse' and 'Neuroactive ligand-receptor Interaction' as the most over-represented pathways in the hypertensive strain (Fig. 6). Within these KEGG enrichment pathways, the genes associated with glutamate signalling were the most heavily over-represented genes, supporting previous evidence that amino acids may play a role in regulating sympathetic activity in healthy physiological conditions 34,36,37 . The NMDAR NR2 subunit (Grin2b) that was found to be upregulated in SHR ganglia, is the binding site for glutamate and plays a direct role in regulating the electrophysiological properties of the channel [38][39][40] . Moreover, decreases in the AMPA subunit GluA2 (Gria2) mRNA expression may be associated with increased Ca 2+ permeability. Indeed, in CA1 pyramidal neurons, increased cAMP and PKA signalling shift GluA2/3-containing AMPARs to a more conductive state 41,42 . Whether impaired cN signalling within SHR PGSNs affects AMPA R conductance and contributes to the observed Ca 2+ dysfunction within sympathetic stellate neurons remains to be established.
Within the over-represented KEGG groups we identified several differentially expressed transcripts indicating that cGMP-PKG signalling may be downregulated. We have previously identified an impairment in cyclic nucleotide (cN) signalling linked to dysregulated [Ca 2+ ] i in hypertensive states 2,3 and have demonstrated that enhancing cGMP signalling rectifies this feature 2 . We observed that the transcripts encoding the soluble guanylyl cyclase 1 soluble subunit α (Gucy1a3) and the cGMP-dependent Protein Kinase (PKG) type II subunit (Prkg2) are significantly downregulated in SHR ganglia. Importantly, we identified that the 'phosphoric ester hydrolase activity' (GO:0042578) was also significantly over-represented within the MF category. Within this GO group, we identified 33 significantly differentially expressed transcripts, including downregulation of Pde11a and Pde2a transcripts and upregulation of Pde6b (Fig. 5) supporting previous evidence that PDE dysregulation may play a role in cardiovascular diseases 13,16,[43][44][45][46] . Recently, it was demonstrated that 1% of the Swedish population carry a loss-of-function mutation in PDE11a that is directly linked to hypertension, obesity and ischaemic stroke 47 . Moreover, a single nucleotide polymorphism (rs197163010) has been identified in PDE2a in the SHR/ola strain and reported in the NCBI dbSNP database 48 . Importantly, a number of other differentially expressed genes within this group also relates to the regulation of cell function through protein phosphorylation, including protein tyrosine and phosphoserine phosphatases that regulate protein kinase effector activity. A full overview of the differentially expressed transcripts are displayed in Fig. 5; however, it is pertinent to note that cN activity is dependent on PDE activity that is closely localised within specific subcellular microdomains 49 .
Several genes that we identified within these pathways have been previously associated with hypertension and CVDs. These include the class III metabotropic glutamate receptor 7 (Grm7) 50,51 , G-protein β3 subunit (Gnb3) [52][53][54][55] and the voltage gated L-type Ca 2+ Cannel α subunit (Cacna1c) 56,57 . Of particular relevance to this study is the finding that Protein Kinase A regulatory subunit 1a (Prkar1a) is downregulated in these cell types. Prkar1a transcription is integral for normal development of peripheral neurons including axonal sorting, myelination and proliferation of neurons 58 . Additional support for the importance of Prkar1a in peripheral neurons is found in patients with the inherited Carney Complex (CNC) syndrome 59,60 . Notably, patients with global loss-of-function mutations in Prkar1a develop peripheral nerve tumours and display a wide range of cardiovascular comorbidities that are often fatal, including dysregulated maintenance of blood pressure, cardiomyopathies and arrhythmia 60 . Indeed, in one study approximately 60% of CNC patients died of cardiovascular complications 59 . Limitations. In this study, we carried out a hypothesis neutral, non-biased approach to sequencing the transcriptome of sympathetic stellate neurons of SHR and normotensive Wistar rats. There are several limitations to this approach. First, the stellate ganglion comprises a heterogeneous population of cell types. This was apparent in the analysis, as we identified several markers of fibroblasts and astrocytes including vimentin and glial fibrillary acidic protein (GFAP) respectively, however; a high number of identified transcripts were of neuronal phenotype. Further studies to fully uncover the transcriptome of PGSNs require single cell RNAseq or confirmation of RNA transcripts in single PGSNs, using techniques such as in situ hybridisation for nucleic acids. Secondly, we used stellates obtained from male rats. While there is a plethora of data regarding the sex differences in hypertension and CVD incidence 61 , in this study we focussed on investigating the transcriptome of the male rat stellate ganglia, given that the prevalence for hypertension and CVD is significantly higher in males compared with premenopausal women 61 . Thirdly, given the low concentration of RNA in our stellate samples, we chose to amplify our samples during the library preparation using the SMARTer amplification protocol. We are aware that low level transcripts may be under-represented in the amplification process and in our dataset. Moreover, since we sequenced the transcriptome from rats with established hypertension, the transcriptomic changes that we identified may occur as compensatory rather than causative factors. Finally, our gene profiling and ontology analyses are enabled, yet limited, by the GO groups and KEGG enrichment pathways currently available.
Notwithstanding these limitations, RNAseq transcriptomics has revealed differential gene expression profiles in PGSNs obtained from SHR with established hypertension and normotensive Wistar rats. Leading transcripts were also present in human stellate ganglia. Taken together, these genes provide a starting framework to further investigate their role in sympathetic hyperactivity and other diseases associated with autonomic dysfunction.

Methods
Animals. 16-week old male SHRs with established hypertension and age-matched normotensive Wistar control rats were obtained from Envigo, UK. At 16-weeks, it is well-established that SHR display hypertension and sympathetic hyperactivity 4,30,62-66 . In this study, we used the Wistar rat strain as the normotensive control, given that Wistar rats are the progenitor strain for WKY and the two strains display similar hemodynamic profiles at all ages as measured using multiple methods 4,66-72 . We have previously shown that neither Wistar nor WKY display a sympathetic Ca 2+ phenotype (unlike the SHR) 13 making the Wistar a suitable control in this study 73 . All rats were housed in standard plastic cages and artificial lighting was fixed to a natural 12-hour light/dark cycle. Food and water were available ad libitum. All experiments were performed in accordance with the UK Home Office Animal Scientific Procedures Act 1986 (ASPA) and approved by the University of Oxford (PPL 30/3131; David J. Paterson).

RNA Extraction from Sympathetic Cardiac Ganglia. Right and left postganglionic sympathetic
stellate neurons were dissected from 16-week-old male SHR and age-matched male Wistar rats. Right stellate ganglia were removed for RNA sequencing and the left stellate ganglia were used for matched qRT-PCR experiments. Rats were anaesthetised in an induction chamber (3-5% isoflurane) and humanely killed by a Home Office approved Schedule 1 method: overdose of pentobarbital (Euthatal, 200 mg/mL) and exsanguination. The stellate ganglia were placed in Hanks Buffered Saline Solution without Ca 2+ and Mg 2+ . For clinical samples, human ganglia (n = 3 patients, two left and two right stellates) were extracted at UCLA and shipped on dry ice in RNAlater ® RNA Stabilization Solution (ThermoFisher). Rat and human ganglia were cleaned and de-sheathed and each sample contained RNA from one stellate. Briefly, stellate ganglia were transferred immediately to RLT lysis buffer (Qiagen) with β-mercaptoethanol (1%). Stellates were finely chopped and carefully triturated with fire-blown Pasteur pipettes until adequately digested. An equal volume of ethanol (70%) was added to each sample and vortexed. Rat RNA was extracted using an RNeasy Mini RNA Extraction Kit (Qiagen) and human RNA was extracted using an RNeasy Maxi RNA Extraction Kit (Qiagen) in accordance with the manufacturer's instructions. RNA samples were divided and aliquoted for qRT-PCR and quality control experiments before snap freezing in liquid nitrogen and storage at −80 °C. The RNA quality and integrity from each sample, was confirmed using a 2100 Bioanalyzer Instrument with an RNA 6000 pico kit (Agilent). Rat samples with an RNA integrity number (RIN) less than 8.5 were discarded. Due to difficulties with fast shipping of human samples we accepted human RNA samples with RINs above 5.7. RNA concentrations were determined using a Qubit RNA High Sensitivity Assay Kit (Molecular Probes, Life Technologies) and a Qubit ® 2.0 Fluorometer (Invitrogen, Life Technologies). Regarding the extraction of human ganglia; all methods were carried out in accordance with relevant guidelines and regulations. The human study was approved by an institutional committee, the UCLA Institutional review board, approval # 12-000701, and informed consent was obtained from all subjects. cDNA library preparation for RNA Sequencing. 4 samples containing total RNA extracted from the right sympathetic stellate ganglion of 16-week-old male SHR (n = 4) and age-matched Wistar (n = 4) were sent to the High-Throughput Genomics Group at the Wellcome Trust Centre for Human Genetics (WTCHG) for RNAseq library construction and sequencing using an Illumina HiSeq 4000 (Illumina, Inc., San Diego, USA). The sequencing libraries were amplified using a SMARTer (first strand synthesis) amplification protocol, due to the low initial RNA concentrations obtained from a single stellate, and prepared for paired-end sequencing (2 × 75 bp). Each sample was sequenced on two separate lanes to minimise technical error and to increase the sequencing depth (~15-25 million reads per lane). Samples were randomized and blinded to the experimenter. The sequencing parameters established, were based on recommendations from WTCHG and those published by Conesa et al. 18 .
Data files were assigned alternative names to blind the experimenter during the relevant stages of the quasi-mapping analysis. Differential Expression Analysis. Once each sample was quantified, the data were imported into R and summarised to a gene-level using the 'Tximport' function (v1.4.0) 74 . Gene level differential expression analysis between Wistar and SHR samples was performed using the 'DESeq2' command in DESeq2 (version 1.16.1) 20 . Significance for differential expression was accepted at the Benjamini-Hochberg adjusted p < 0.05 level. The 'LFCshrink' function was used to shrink log fold changes after analysis as per the DESeq2 vignette 20 .
Gene ontology and KEGG analysis. Gene ontology (GO) and KEGG over-representation tests were performed using ClusterProfiler (Version 3.4.4) with 'EnrichGO' and 'EnrichKEGG' functions respectively 22 . The differentially expressed gene input was taken at the Benjamini-Hochberg adjusted p < 0.05 level. Significant groups were subsequently determined using Benjamini-Hochberg correction and cut-off levels of p < 0.01 and q < 0.05 were enforced. Redundancy amongst GO terms was removed using the ClusterProfiler simplify command, where the adjusted p-value was used as the deciding variable, with a cut-off value of 0.7. The ClusterProfiler 'Enrichmap' function was used to generate KEGG enrichment maps [23][24][25] . We subsequently generated an integrated KEGG enrichment diagram by combining signalling pathways from the over-represented KEGG pathways 'Circadian entrainment' , 'Dopaminergic synapse' and 'Glutamatergic synapse' .
cDNA Library Preparation for qRT-PCR. We  (GStorm) at 25 °C (10 mins), 42 °C (60 mins), 80 °C (5 mins). For human PGSN RNA, the SuperScript TM IV VILO TM cDNA synthesis protocol was followed according to manufacturer's instructions (ThermoFisher). Briefly, human RNA samples (50 ng) were added to ezDNAse buffer, exDNAse enzyme and nuclease-free H 2 O. The samples were mixed and incubated (37 °C, 2 min). SuperScript TM IV VILO TM master mix (or for control experiments SuperScript TM IV VILO TM no reverse transcriptase) was added to each sample with nuclease-free H 2 O. The samples were mixed and incubated in a thermocycler (GStorm, GS4822) at 25 °C (10 mins), 50 °C (10 mins), 85 °C (5 mins). Concentrations of cDNA in each sample and the 260/280 ratios were calculated (NanoDrop Lite) to detect the presence of contaminants. Samples with an abnormal 260/280 ratio (<1.7 and >1.95) were discarded. Samples were aliquoted and frozen at −80 °C for long-term storage or kept at 4 °C for immediate use.
Statistical Analysis. qRT-PCR data were analysed using GraphPad Prism software (v7). When the data passed normality tests, an unpaired Student's t-Test was used. When the data were not normally distributed, nonparametric tests were used with the specific statistical test reported in the figure legend. All qRT-PCR data are expressed as Log 2 fold change ± SEM. Statistical significance was accepted at p < 0.05, unless otherwise described. For details of quality control and further methodology, please refer to 76 . Data availability. Our RNA sequencing (RNAseq) raw FastQ files are deposited in the NCBI short reads archive under SRA number SRP132271 (NCBI Short Read Archive) and our quasi-mapped data will be available under Gene Expression Omnibus (GEO) accession number (GSE110197).