Identifying Candidate Genes that Underlie Cellular pH Sensitivity in Serotonin Neurons Using Transcriptomics: A Potential Role for Kir5.1 Channels

Ventilation is continuously adjusted by a neural network to maintain blood gases and pH. Acute CO2 and/or pH regulation requires neural feedback from brainstem cells that encode CO2/pH to modulate ventilation, including but not limited to brainstem serotonin (5-HT) neurons. Brainstem 5-HT neurons modulate ventilation and are stimulated by hypercapnic acidosis, the sensitivity of which increases with increasing postnatal age. The proper function of brainstem 5-HT neurons, particularly during post-natal development is critical given that multiple abnormalities in the 5-HT system have been identified in victims of Sudden Infant Death Syndrome. Here, we tested the hypothesis that there are age-dependent increases in expression of pH-sensitive ion channels in brainstem 5-HT neurons, which may underlie their cellular CO2/pH sensitivity. Midline raphe neurons were acutely dissociated from neonatal and mature transgenic SSePet-eGFP rats [which have enhanced green fluorescent protein (eGFP) expression in all 5-HT neurons] and sorted with fluorescence-activated cell sorting (FACS) into 5-HT-enriched and non-5-HT cell pools for subsequent RNA extraction, cDNA library preparation and RNA sequencing. Overlapping differential expression analyses pointed to age-dependent shifts in multiple ion channels, including but not limited to the pH-sensitive potassium ion (K+) channel genes kcnj10 (Kir4.1), kcnj16 (Kir5.1), kcnk1 (TWIK-1), kcnk3 (TASK-1) and kcnk9 (TASK-3). Intracellular contents isolated from single adult eGFP+ 5-HT neurons confirmed gene expression of Kir4.1, Kir5.1 and other K+ channels, but also showed heterogeneity in the expression of multiple genes. 5-HT neuron-enriched cell pools from selected post-natal ages showed increases in Kir4.1, Kir5.1, and TWIK-1, fitting with age-dependent increases in Kir4.1 and Kir5.1 protein expression in raphe tissue samples. Immunofluorescence imaging confirmed Kir5.1 protein was co-localized to brainstem neurons and glia including 5-HT neurons as expected. However, Kir4.1 protein expression was restricted to glia, suggesting that it may not contribute to 5-HT neuron pH sensitivity. Although there are caveats to this approach, the data suggest that pH-sensitive Kir5.1 channels may underlie cellular CO2/pH chemosensitivity in brainstem 5-HT neurons.


INTRODUCTION
Chronic physiological regulation of arterial pH is controlled by the kidney, but arterial oxygen, carbon dioxide (CO 2 ), and pH levels are acutely regulated by chemoreflex-mediated changes in breathing. Central (brain) respiratory chemoreceptors are resident CNS cells (neurons/glia) whose pH-or CO 2 -dependent activity can affect ventilation (Richerson et al., 2005). Some brainstem cells have been identified as CO 2 /pH-sensitive and involved in the control of breathing, and in some cases the location, neurochemical identity, and molecular mechanism(s) of pH sensitivity have been determined. For example, accumulating data has demonstrated that CO 2 /pH-sensitive glutamatergic (Phox2b + ) neurons in the retrotrapezoid nuclei (RTN) owe their intrinsic pH sensitivity to the combined expression of a specific Twik-related acid-sensitive potassium (K + ) ion channel (TASK-2), and a novel G protein-coupled receptor (GPR4) (Kumar et al., 2015). The activity of Phox2b + RTN neurons is also modulated by raphe-derived neuromodulators (serotonin and substance P), and via ATP released by local pH-sensitive glia (Mulkey et al., 2007a), which appear to rely on heteromeric, inwardly rectifying, pH-sensitive K + (Kir4.1/Kir5.1) channels for their chemosensitivity (Wenker et al., 2010). Other populations of CO 2 -sensitive neurons in the Locus Coeruleus (LC; catecholaminergic neurons) and solitary complex are also thought to rely on yet to be identified pHsensitive K + channels (Martino and Putnam, 2007;Li and Putnam, 2013).
Serotonin-producing (5-HT) raphe neurons have also been shown to be pH-and CO 2 -sensitive, directly responding to hypercapnic and/or metabolic acidosis with increases in action potential firing rates in a variety of conditions and experimental preparations including behaving animals (Veasey et al., 1995(Veasey et al., , 1997Wang et al., 2001;Bradley et al., 2002;Washburn et al., 2002;Iceman et al., 2013). The cellular response of 5-HT neurons to hypercapnic acidosis increases in magnitude on or after postnatal (P) day 12 in rats in vitro (Wang and Richerson, 1999;Wu et al., 2008). Specific subpopulations of 5-HT neurons in the rostral medullary raphe, which arise from the Erg2-expressing developmental cell lineage display an intrinsic pH sensitivity, and selectively innervate populations of brainstem neurons in respiratory-related nuclei (Brust et al., 2014). A modest pH sensitivity of young (<10 days of age) 5-HT neurons may depend upon their expression of pH-sensitive leak K + channels such as TASK-1 and TASK-3 channels, although these channels do not appear to contribute to the CO 2 chemoreflex in mature mice (Washburn et al., 2002;Mulkey et al., 2007b). Thus, the molecular determinants of the mature CO 2 /pH response in brainstem 5-HT neurons remain unclear Teran et al., 2014).
Failure of mechanisms involved in cardiorespiratory control, and specifically dysfunction in the brainstem 5-HT system and the CO 2 chemoreflex are believed to contribute to the pathophysiology of Sudden Infant Death Syndrome (SIDS; Panigrahy et al., 2000;Kinney et al., 2003Paterson et al., 2006). SIDS is a leading cause of post-neonatal infant mortality in USA, and is thought to result from exogenous stressors applied during a critical period of development in an infant with an underlying abnormality (Filiano and Kinney, 1994). Several 5-HT system abnormalities have been identified in SIDS cases (Paterson et al., 2006;Duncan et al., 2010), and thus defining the developmental and regional shifts in gene expression profiles within brainstem 5-HT neurons and the potential genetic determinants of cellular CO 2 /pH chemosensitivity could provide key insights into the underlying biological dysfunction in SIDS.
Our goal here was to use age-related changes in gene expression profiles in groups of brainstem 5-HT neurons to identify developmentally regulated genes that may encode putative pH-sensitive molecules, including ion channels. By combining fluorescence-activated cell sorting (FACS) of primary eGFP-expressing 5-HT neurons from an established transgenic rat (SS ePet−eGFP rats; Katter et al., 2013) and subsequent mRNA sequencing and differential gene expression analyses, we identified several genes that may encode proteins that contribute to cellular pH sensitivity in 5-HT neurons. In particular, the gene and protein expression of pH-sensitive Kir4.1 and Kir5.1 channel subunits increased with age in 5-HT neuron-enriched cell pools and raphe tissues, respectively. However, mRNA and/or protein for other pH sensitive K + channels (TASK-1 and TASK-3) decreased or did not change with increasing age in 5-HT cell pools/raphe tissues. Single cell qPCR in adult eGFP + 5-HT neurons confirmed heterogeneity of gene expression of Kir4.1, Kir5.1, TASK-1, TASK-3, and other genes. Immunofluorescence staining of adult brainstems showed Kir5.1 protein co-localized with brainstem neurons (including 5-HT neurons) and astrocytes, but Kir4.1 protein expression was restricted to brainstem astrocytes. Our data suggest a potential novel role for Kir5.1 within the medullary raphe. Additionally, our data suggest transcriptome comparisons can provide candidate genes that contribute to cellular pH sensitivity, but also suggest that more specific methods such as single-cell transcriptomics and subsequent RNA sequencing could provide an improved methodology for this purpose.

MATERIALS AND METHODS
Transgenic Dahl SS (SSMcwi) rats carrying the T2/ePet-eGFP transgene (Katter et al., 2013), which directs enhanced green fluorescent protein (eGFP) expression to Pet-1-expressing neurons [(SS-TgTn(T2ePet-eGFP)A1Mcwi rats; RGD strain ID: 8661234], were bred and maintained in-house and used for this study (n = 66). These rats have eGFP expression restricted to all 5-HT neurons, and will be referred to as SS ePet−eGFP hereafter. Additional wild type Dahl Salt-sensitive (SS/JrHsdMcwi; SS) rats, an inbred strain generated and maintained by the Medical College of Wisconsin were used as control animals for this study (n = 6). All rats were housed in the Biomedical Research Center and allowed access to chow (Teklad) and water ad libitum, and maintained on a 12:12 h light/dark cycle. All experimental protocols were approved by the Medical College of Wisconsin Institutional Animal Care and Use Committee prior to experimental use.

Cell Dissociation and Fluorescent Assisted Cell Sorting (FACS)
Young (P0-P3) or mature (P42-60) SS ePet−eGFP rats were anesthetized with isoflurane (20% in propylene glycol v/v), rapidly decapitated, whole brains removed and placed in icecold Hibernate A (Life Technologies, #A1247501) for surgical isolation of the medullary raphe similar to that previously described (Wang et al., 1998). Tissue "wedges" were sectioned into rostral and caudal halves, placed in fresh Hibernate A, minced into small pieces and incubated in ice cold 2 mL Hibernate A or 2 mL of Accutase (Millipore; 30 min), pelleted at 380 g (2 min), re-suspended in 3 mL of Hibernate A and mechanically dissociated with a large-diameter sterile polyethylene transfer pipet similar to methods previously described by others (Guez-Barber et al., 2012). Tissues were gently titrated three to five times, and 2 ml of supernatant was transferred to fresh Hibernate A, and the process repeated until dissociation was complete. Dissociated cells were pelleted at 480 g (3 min) and re-suspended in Hibernate A and filtered (40 µm) and 4 , 6-diamidino-2-phenylindole (DAPI) added as a dead/live indicator prior to cell sorting [FACS Aria IIIa (BD Biosciences)]. Collection gating was established based on neuronal health (DAPI), forward and side scatter (size and complexity, respectively), by negative control tissues obtained from naïve SS rats to determine the threshold for eGFP (488 nm wavelength), and control cell samples incubated with anti-NeuN (neuron-specific nuclear protein, Abcam #AB104225) to determine neurons from debris. All cells subjected to sorting into eGFP + or eGFP − pools were collected in Trizol (Life Technologies; 5000-8000 cells per pool), where tissues from three young rats were combined for dissociation per cell sort, and one adult rat was used for each cell sort.
All differentially expressed transcripts from a given comparison were further analyzed using Ingenuity Pathway Analysis software. We focused on differentially expressed transcripts between YRP vs. ARP, YCP vs. ACP, ACP vs. ACM, ARP vs. ARM, YCP vs. YCM, and YRP vs. YRM and performed "core analysis" for each which provided outputs including Top Canonical Pathways, Upstream Analysis, Diseases and Functions, Regulator Effects, and Networks. An in-depth analysis of all upregulated plasma membrane molecule transcripts reaching significance (P > 0.01) between each library comparison as well as a targeted comparison of all differentially expressed K + channels was performed.

qRT-PCR Validation
To validate the specificity of the cell sort and age-dependent changes in gene expression, we used quantitative real-time PCR on a separate set of FACS-isolated eGFP + and eGFP − cells to quantify expression of the 5-HT neuron-specific genes Tph2 and Slc6a4, and a housekeeping gene (18s) as an internal control for normalization.
To validate differentially expressed transcripts generated through mRNA-sequencing, we used quantitative PCR on a separate set of FACS-isolated eGFP + and eGFP − cells from a separate six animals each for four different age groups [P0-1, P7-8, P19-20, and adult (>P42)] of SS ePet−eGFP animals using primers to quantify multiple transcripts including Kcnj10 (Kir4.1), Kcnj16 (Kir5.1), Kcnk3 (Task-1), Kcnk9 (TASK-3), Kcna2 (Kv1.2), and 18S (internal control). Total RNA was extracted from FACS-isolated cells with Trizol (described above). Samples were then quantified and qualified (Nanodrop), converted to cDNA (RevertAid First Strand cDNA Synthesis Kit, #K1622, Thermo Scientific) and assayed using SYBR Green (iTaq Universal SYBR Green Supermix, Bio-Rad) technology as previously described (Yao et al., 2005) in triplicate. Gene expression was calculated using the delta-delta Ct method ( C T ) where 18s was the housekeeping control. C T of GFP + cells was the experimental group. C T of GFP − or P0 cells were the control groups. Data are expressed as log fold change. Primer and probe sequences are shown in Supplementary Table 1.

Isolation of Single Cell Contents and qRT-PCR
To verify 5-HT neuronal expression of candidate genes derived from FACS-RNA Sequencing, single intracellular 5-HT neuronal contents were isolated and gene expression for Kcnj10, Kcnj16, Kcnk1, Kcnk3, Kcnk5, Kcnk9, and Kcna2 was measured using a TaqMan chemistry based Single Cell-to-C T kit (#4458236, Ambion). SS ePet−eGFP rat pup (P18-23) brains were isolated and then dissected in ice cold oxygenated (95% O 2 /5% CO 2 ) slice solution (220 mM sucrose, 10 mM MgCl 2 , 2.5 mM KCl, 0.2 mM CaCl 2 , 10 mM dextrose, 1.3 mM NaPO 4 , and 26 mM NaHCO 3 at a pH of 7.4) and then coronally sectioned (300 µm) using a vibratome. Rostral raphe sections were incubated for 30 min in nine parts slice solution and one part Ringers solution (124 mM NaCl, 3 mM KCl, 2 mM MgCl 2 , 2 mM CaCl 2 , 10 mM dextrose, 1.3 mM NaH 2 PO 4 , and 26 mM NaHCO 3 ) before being placed into an electrophysiology recording chamber (no recordings were made). Single 5-HT neurons were identified under fluorescent live imaging. 0.5 µm diameter sterile micropipettes were backfilled with <1 µl of a modified intracellular solution according to Cadwell et al., 2016 [123 mM potassium gluconate, 12 mM KCl, 10 mM HEPES, 0.2 mM EGTA, 4 mM MgATP, 0.3 mM NaGTP, 10 mM sodium phosphocreatine, 20 µg/mL glycogen, and 1 U/µl recombinant RNase Inhibitor (#2313A, Takara) with pH 7.25] and using a micromanipulator (P97; Sutter Instruments), were placed onto the surface of single 5-HT fluorescent neurons such that there was a visible deflection or dimple in the membrane (seal resistance was not measured). Back pressure was applied and intracellular contents were collected using pulses negative pressure, which were visually verified using live fluorescent imaging (Cool Snap Nikon camera and Nikon elements D). Micropipettes were gently pulled out of the tissue and the micropipette tips (∼1 cm) were placed into an Alconox detergent solution (20 mg/ml) for 1 min to remove potentially attached tissue debris. The micropipette tip was then placed in Ringers solution for 1 min to remove the detergent. Isolated contents were then expelled into Lysis/DNase buffer provided by the Single Cell-to-CT kit (9 µl of lysis buffer and 1 µl of DNase), incubated for 5 min followed by a 2-min incubation after 1 µl of Stop Solution was added. Each sample was stored at −20 • C until all samples were collected and then processed according to the Single Cell-to-CT kit. We confirmed that washing the micropipette tip had no unwanted effects and actually improved the purity of our samples by comparing gene expression of 5-HT (Tph2, Ddc), neuronal (Map2, Nefl), and glial (Gfap, Aqp4, and Sox10) specific gene markers and Kcnj10 and Kcnj16 from a subset of samples washed with Ringers + Ringers (n = 5 cells) and with Alconox + Ringers (n = 6 cells). There was no difference in Ct-values between washing methods. We collected another seven samples that we used to measure the same 5-HT, neuronal, and glial markers in addition to K + channel genes (Knj10, Kcnj16, Kcnk1, Kcnk3, Kcnk5, Kcnk9, and Kcna2). To gauge the expression pattern of Kcnj10 and Kcnj16 across 5-HT neurons we pooled all cells together (n = 18) and determined if there was or was not expression of these genes. A sample of bulk tissue was obtained and RNA isolated using a kit (74034, Qiagen) in order to compare 5-HT, neuronal and glial gene expression patterns between bulk tissue and our single cell isolation approach. qRT-PCR was run on a Quant Studio 6 Flex (Life Technologies).
Raw 8 or 16 bit monochrome images of the red and green channels were also processed for each combination of proteins labeled to determine co-localization in ImageJ by calculating Manders' correlation coefficients. This was calculated by determining the fraction of pixels of protein A overlapping with protein B, and the reverse calculation to derive correlation coefficients where 1.0 is a perfect co-localization and 0 no colocalization.

Enriching for eGFP-Expressing Medullary Raphe 5-HT Neurons with FACS
In order to derive candidate molecules underlying cellular CO 2 /pH sensitivity, we measured age-related changes and regional differences in gene expression in brainstem 5-HT neurons. We took advantage of an established rat line that carries an eGFP transgene under the control of an enhancer element for a 5-HT neuron-specific transcription factor Pet-1 (Katter et al., 2013). Brainstems from SS ePet−eGFP rats were extracted, eGFP + (5-HT) and eGFP − (non-5-HT) neurons were dissociated then isolated into separate pools using FACS, similar to methods previously published (Guez-Barber et al., 2012). Conservative eGFP + gating was established by sorting unlabeled medullary raphe cells from naïve SS rats, allowing for the identification of the putative neuronal population based on the resulting forwardand side-scatter data (Figures 1A,B). The scatterplot revealed a distinct cell population that corresponds to isolated neurons (Guez-Barber et al., 2012), which we confirmed in our sorted cell dissociates with or without fixation and conjugated primary antibodies specific for neurons (NeuN/Alexa 546; Figure 1B). This gated sub-population (P1) ranged from 68% to as high as 99.9% NeuN positive cells (data not shown), and from this collected pool we estimated that eGFP + neurons represented 0.9-1.8% of these cells, supporting our assumption that we had established conservative gates for eGFP + cells (Figures 1C,C'). Total RNA extracted from the FAC-sorted eGFP + and eGFP − cell pools was of high quality and sufficient quantity (Nanodrop 2000), and RT-qPCR experiments validated the specificity of the sort. We found 5-HT neuron specific genes (Slc6a4 and Tph2), were significantly more expressed in eGFP + (5-HT + ) vs eGFP − (5-HT − ) cell pools of adult rats (Figure 1D), indicating significant 5-HT neuronal enrichment using this method.

RNA Sequencing Reveals Large Developmental Shifts in Gene Expression in Medullary 5-HT Neurons
After validating our FACS enrichment strategy for brainstem 5-HT neurons, cDNA libraries for RNA Sequencing were prepared from FAC-sorted cell pools (20,000-70,000 cells/pool) collected from young (Y) or adult (A) rats, rostral (R) or caudal (C) raphe regions that were either eGFP-positive (P) or eGPFnegative (N; YRP, YRN, YCP, YCN, ARP, ARN, ACP, ACN). We noted that biological replicates within each experimental group yielded similar numbers of transcripts detected, and across all cDNA libraries there were high transcript yields (2.25-4.00 Gb), read number (23.1-41.0 million/library), mean quality score (94.2 or greater), and mapping rate (89.5% or better; data not shown). Libraries from common experimental groups were combined and differentially expressed genes determined across group comparisons using an adjusted p-value of q < 0.05 as the threshold for significance. We noted a far greater number of differentially expressed genes when making comparisons across ages (6234.3 ± 200.5 genes), and relatively few differences comparing cell types or regions within an age group (183.8 ± 44.8 genes). This major difference could reflect that these cells continue to develop after birth into adulthood, whereas the relatively small differences across cell type and region might suggest significant heterogeneity in gene expression profiles that reduce the statistical power to identify differential expression among all genes. Regardless, changes in gene expression patterns are greatly affected by age compared to cell type and raphe region (caudal vs. rostral).
Given the preponderance of K + channels in the above comparisons and our assumption that the genetic determinants of cellular chemosensitivity will increase in expression with increasing age, we focused on all K + channels found in 5-HT neuron-enriched cell pools and determined which of these channel genes were expressed more highly in 5-HT versus non-5-HT neuron cell pools ( Figure 3B). Among the agerelated changes in expression of all K + channels detected in the rostral raphe 5-HT neurons, we found 25 genes were downregulated including TASK-1 and TASK-3 (red; Figure 3B), 13 were upregulated including TWIK-1, and Kir4.1 and Kir5.1 (green; Figure 3B), and 31 were unchanged with age. However, expressing these fold-changes in expression relative to their FIGURE 2 | All plasma membrane-associated genes that were upregulated with increasing age in 5-HT neuron-enriched cell pools. Plasma Membrane molecules identified in IPA were significantly increased (q < 0.01) with increasing age in 5-HT neuron-enriched cell pools by FACS. (A) Gene name, description, protein type and relative change in expression (greatest increase (red) to smallest increase (yellow) for log fold-change with increasing age) is shown for identified plasma membrane-associated molecules. (B) Representation of the type of plasma membrane molecules categorized by transporters, ion channel, G-protein coupled receptors (GPCR), enzymes, or other. abundance in adult 5-HT neurons ( Figure 3B) suggested that Kir4.1 and/or Kir5.1 may be ideal candidates to underlie pH sensitivity in 5-HT neurons as they were both increased with age and in much higher abundance relative to all other K + channels expressed.

Single Cell PCR Validation Confirms 5-HT Neuron Specific Expression of Kcnj10 and Kcnj16
The RNA-Seq data showed significantly greater abundance (q < 0.05) of 5-HT neuron-specific genes including Tph2, Ddc, and Slc6a4 (SERT) in eGFP + cell pools versus eGFP − pools from the rostral raphe of mature rats, confirming successful 5-HT enrichment. However, we did not measure significantly greater 5-HT neuron-specific genes in our other group comparisons. For unknown reasons, 5-HT neuron enrichment was not as robust in cells collected from the caudal raphe of young or mature rats, similar to other reports (Wylie et al., 2010). To our surprise, we also found multiple glial-specific genes in all 5-HT neuronenriched cell pools, including Gfap, Aqp4, and Sox10. Thus, the sequencing data indicated our dissociation and FACS methods had enriched for 5-HT neurons but also likely included fragments of neighboring glial cells, which brought into question the 5-HT neuron-specific expression of our candidate genes.
Developmental Increases in Kcnj10/Kir4.1 and Kcnj16/Kir5.1 in 5-HT Neuron-Enriched Cell Pools Chemosensitivity increases with increasing age in 5-HT neurons in vitro, and thus we hypothesized that candidate genes that FIGURE 3 | Assumption-based identification of candidate genes through unique and overlapping differential expression patterns highlight the potential involvement of potassium ion (K + ) channels. (A) Upregulated plasma membrane molecules [>0.5 log-fold increase) from age-specific (young vs adult) rostral 5-HT RNA sequencing comparison; (blue)] and region-specific [caudal versus rostral adult 5-HT neurons; (yellow)] with common genes in middle (tan). Note that the genes listed in red encode proteins with known pH sensitivity, and those with asterisks bind to and/or associate with pH-sensitive proteins. (B) Fold change in expression (young vs. old) of all K + channel genes plotted against overall abundance (fragments per kilobase of transcript per million reads; FPKM), where those in red are downregulated, green are upregulated, and those in blue are not differentially regulated with age.
underlie cellular chemosensitivity may also increase in expression with increasing age. To better define the temporal expression pattern of the identified K + channels, we performed additional FACS experiments in SS ePet−eGFP rats at 0, 7-8, 18-19, and 42-60 days of age. PCR amplification of selected genes in the 5-HT neuron-enriched cell pools was normalized to 18s and expressed as a fold-change ( C T ) from P0 values ( Figure 5A). There were significant increases in expression of kcnj10, kcnj16, and kcna2 in 5-HT neuron-enriched cell pools with increasing age (Figure 5A; P < 0.05; One-way ANOVA), where expression levels were increased by P19 or after. In contrast, mRNA expression levels of both kcnk3 and kcnk9 significantly decreased with increasing age in 5-HT neuron-enriched cell pools (Figure 5; P < 0.05; One-way ANOVA). These PCR data differ somewhat from the single cell PCR data from individual adult 5-HT neurons (Figure 4E), where there was robust gene expression of Kcnk3 and Kcnk9, which is more consistent with previous reports (Washburn et al., 2002).

Localization of Kir4.1 and Kir5.1 Channel Proteins in the Medullary Raphe
Kir4.1 and Kir5.1 channel subunits are largely thought to be co-expressed in glial cells throughout the CNS, although there is some evidence for Kir5.1 channel expression in neurons. Figure 6 shows images of immunofluorescence of fixed-frozen brainstem sections labeled with primary antibodies targeting Kir4.1 (Figures 6A-F) or Kir5.1 (Figures 6G-L) along with neuronal (NeuN) and astrocytic (GFAP) markers. Kir4.1 expression in the medullary raphe (and throughout the brainstem) was restricted to GFAP-expression astrocytes (Figures 6A-C), and was not found co-localized with the neuronal marker NeuN (Figures 6D-F). Kir5.1 protein was also found co-localized with the vast majority of GFAPexpression astrocytes in the medullary raphe (Figures 6G-I), but was also found in cells other than glia. Kir5.1 expression was co-localized to many neurons in the brainstem, and in particular appeared to be co-localized mostly with neurons in the rostral medullary raphe (Figures 6J-L). Consistent with these observations, there was no observable Kir4.1 protein that co-localized with 5-HT neurons (Figures 7A-C), but we found co-localization of Kir5.1 with 5-HT neurons in the medullary raphe (Figures 7D-F), where the 5-HT neuronspecific expression of Kir5.1 appeared to be perinuclear and to a lesser extend membrane bound (Figures 7G-I). These data demonstrate the presence of these Kir isoforms within the raphe, but suggest only Kir5.1 protein was expressed in 5-HT (and other) neurons while Kir4.1 protein was localized to neighboring glia.

DISCUSSION
The molecular mechanisms conferring intrinsic CO 2 and/or pH chemosensitivity of medullary raphe 5-HT neurons remain unknown. The ventilatory response to hypercapnia has generally been shown to increase throughout post-natal development in rodents and humans (Moss and Inman, 1989), and the intrinsic CO 2 /pH sensitivity of brainstem 5-HT neurons appear to also increase throughout post-natal development (Wang and Richerson, 1999;Wu et al., 2008). SIDS appears to result due to cardiorespiratory control failure (Thach, 2005), particularly during a critical window of development where SIDS rates are highest (Kinney and Thach, 2009). Understanding fundamental properties of 5-HT neurons, developmental aspects of their function, and their role in cardiorespiratory control mechanisms is critical to improving our understanding of SIDS given the identification of several 5-HT system abnormalities reported in SIDS cases (Paterson et al., 2006;Duncan et al., 2010).
Herein, we tested the hypothesis that gene transcripts encoding pH sensitive proteins, including known pH-sensitive ion channels are developmentally regulated in 5-HT neurons by applying fluorescent cell sorting and subsequent RNA sequencing to determine age-related and region-specific differential gene expression profiles. Further, we developed a robust technique to isolate the intracellular contents of individual adult 5-HT neurons to validate the mRNA expression of some of the candidate genes identified from RNA Sequencing. The data showed that adult brainstem 5-HT neurons express genes that encode pH sensitive molecules, including ion channels. Candidate genes were selected based on three major underlying assumptions, that: (1) genes encoding the CO 2and/or pH-sensitive proteins driving cellular chemosensitivity FIGURE 7 | Kir5.1 but not Kir4.1 protein is co-localized to brainstem 5-HT neurons. Epifluorescence images of fixed-frozen brainstem sections labeled with primary antibodies targeting Kir4.1 (green; A) and TPH (red; B) or Kir5.1 (green; D) and TPH (red; E). Overlay images of indicate that Kir5.1 (correlation coefficient = 0.605; F) but not Kir4.1 (correlation coefficient = 0.109; C) is co-localized to TPH + 5-HT neurons (yellow), where TPH and Kir5.1 co-localization appears primarily perinuclear and to a lesser extent along the plasma membrane (G-I). Confocal images (63×, 2.5 magnification) of brainstem tissues double-labeled with TPH (red; G) and Kir5.1 (green; H) also show intense co-localization (yellow; I) along the perinuclear and to a lesser extent plasma membrane regions. Scale bar = 20 µm (C,F) or 10 µm (I).
in 5-HT neurons increase in expression with increasing age (Wang and Richerson, 1999;Wang et al., 2001), (2) rostral raphe (magnus) regions contain significantly greater numbers of chemosensitive 5-HT neurons compared to caudal raphe (obscurus) regions (Brust et al., 2014), and (3) that the candidate genes would encode proteins associated with or embedded within the extracellular membrane. These assumptions gave rise to transcriptome comparisons from 5-HT neuron-enriched cell pools across age and raphe regions, allowing lists of >6000 developmentally regulated genes to be reduced to <100 candidate genes. Further restriction of the gene list to K + ion channels alone, we found that of the ∼70 K + channel genes expressed in 5-HT-enriched cell pools, there were a select few K + channel genes that increased with age and were in high abundance; Kcnj10 and Kcnj16. These genes were of particular interest as they encode the inwardly rectifying K + channel subunits Kir4.1 and Kir5.1, respectively, which can form heteromeric "leak" K + channel with pHsensitivity well within the physiologic range (pKa ∼6.8 -7.5;reviewed in Putnam et al., 2004). Gene expression of both Kcnj10 and Kcnj16 increased throughout postnatal age in 5-HT neuron-enriched cell pools with FACS, which agreed well with the developmental increases in protein expression of Kir4.1 and Kir5.1 measured in raphe tissues. Finally, qPCR from the intracellular extracts from individual adult eGFP + 5-HT neurons also showed expression of Kcnj10, Kcnj16 or both genes, further suggesting that these genes may encode pH-sensitive K + channels that could underlie cellular pH sensitivity of chemosensitive brainstem 5-HT neurons.
Although Kir4.1/5.1 channels and others containing Kir4.1 or Kir5.1 subunit-containing channels have previously been postulated (Putnam et al., 2004;Wu et al., 2004) or demonstrated to contribute to cellular pH sensitivity (Wenker et al., 2010;D'Adamo et al., 2011) in the CNS, their expression has been considered to be mainly if not exclusively expressed in glial cells (reviewed in Hibino et al., 2010). However, there are reports demonstrating mRNAs for both Kir4.1 and Kir5.1 in brainstem neurons (Wu et al., 2004) and specifically in FACsorted pools of 5-HT neurons (Wylie et al., 2010) or in individual mouse 5-HT neurons (Okaty et al., 2015), consistent with our RNA sequencing and single cell PCR data. However, the immunofluorescence staining in adult rat brainstems showed that Kir4.1 protein expression appeared exclusively glial, with complete co-localization with GFAP expression without any observable co-localization with NeuN expression. In contrast, we found Kir5.1 expression co-localized with GFAP throughout the brainstem, but also found Kir5.1 expression in many brainstem neurons. Kir5.1 (but not Kir4.1) immunoreactivity was also co-localized with raphe 5-HT neurons, which further supports the nomination of Kcnj16 (Kir5.1) subunit-containing channels as potential contributors to cellular pH-sensitivity in 5-HT neurons. Confocal microscopy indicated that the somatic subcellular localization for Kir5.1 was largely perinuclear and to a lesser extent distributed to the cell membrane, although this would need to be assessed with electron microscopy to be certain.
Additional functional evidence is needed for complete validation of the hypothesized role of Kir5.1 as pH sensors in 5-HT neurons, but there are data that suggest Kir5.1 channels contribute greatly to cellular pH responses in other chemosensitive neurons. Cellular responses to acidosis in catecholaminergic LC neurons rely on 4-AP-sensitive (A-type) currents implicating Kir channel involvement (Pineda and Aghajanian, 1997;Li and Putnam, 2013). Furthermore, D'Adamo et al. (2011) showed a reduced cellular responses to reductions in pHi upon removal of prepulses of NH 4 Cl in LC neurons in Kir5.1 knockout mice, and ventilatory data from Kir5.1 knockout mice show a blunted hypercapnic ventilatory response under physiologic conditions in mice challenged with 21% O 2 (although not with hyperoxic hypercapnia) (Trapp et al., 2011), consistent with a role of Kir5.1 channels in the ventilatory CO 2 chemoreflex (although this was not the ultimate conclusion of that study). Thus, Kir5.1 subunit-containing channels contribute to cellular pH responses in aminergic LC neurons and to the ventilatory CO 2 chemoreflex in vivo, but functional contributions to pH sensitivity in 5-HT neurons remains unknown.

Caveats and Considerations
There are several caveats and considerations with the current data, in particular the difficulties in determining the specificity of our cell sorting methods based on gene expression. We took advantage of established methods previously reported for sorting primary dissociated, antibody-tagged CNS cells (Guez-Barber et al., 2011;Guez-Barber et al., 2012) for our eGFP-based approach. Initial gating was based on cell viability (DAPI) and neuronal populations targeted by side-scatter selection confirmed by NeuN-tagged cell sorts. Further, conservative gates were set for eGFP inclusion in an effort to purify the collected 5-HT cell pools. However, the sequencing data suggest that for unknown reasons, the 5-HT neuron enrichment methods with FACS worked better for some conditions than others. In addition, we could readily find glial-specific genes in our 5-HT neuron-enriched cell pools, which could suggested that glial cells were likely included in our 5-HT neuron-enriched cell pools. This conclusion is, however, counterbalanced against the results of our single cell PCR data, which not only demonstrate a high degree of specificity for 5-HT neuronspecific gene expression (11 of 11 expressed Tph2, Ddc, and Nefl), they also show that many 5-HT neurons express mRNAs for proteins considered glial-specific, including Gfap and Aqp4 (9 of 11) and Sox10 (3 of 11). In addition, other groups have shown single cell transcriptome data from mature mouse 5-HT neurons, in which they found similar expression patterns of genes considered to be glial-specific (Okaty et al., 2015). Thus, the presence of "glial-specific" genes in the sequencing data do not necessarily speak to the efficacy of the cell sorting, but rather it speaks to a potential disconnect with measurements of gene expression and whether that translates to protein expression. Our data nominated both Kir4.1 and Kir5.1, but only one of these highly expressed genes gave rise to the corresponding protein in 5-HT neurons. There may be small amounts of Kir4.1 proteins being expressed, but it was not detectible with immunofluorescence.
Despite the inherent disconnect between gene and protein expression, the use of transcriptome comparisons to identify molecular determinants that distinguish between cellular properties or populations may still be useful as a discoverybased approach. For example, the molecular basis for infrared radiation detection in venomous pit vipers (TRPA1) was identified using transcriptome comparisons of trigeminal ganglia from rattlesnakes and other non-pit species. With regard to the molecular basis for chemosensitivity in brainstem 5-HT neurons, the approach may require more specificity than cell sorting of primary neurons can provide, in part due to known heterogeneity in this cell population and for the reasons stated above. Our single cell PCR data highlight this cell-to-cell variation in gene expression, which is also reflected in other single cell transcriptome data (Okaty et al., 2015). An improvement therefore would be to compare gene expression profiles among brainstem 5-HT neurons based on function (pH-sensitive vs. -insensitive for example) at the single cell or subpopulation level to enhance the identification of true positives. Finally, any approach using differential gene expression to identify the molecular basis of a function requires functional validation in addition to PCR and protein expression validation experiments. Thus, our gene candidate Kcnj16 (Kir5.1) requires additional functional validation in vitro and in vivo.
Another consideration is that functional Kir5.1 channels are thought to require heteromeric assembly with other Kir subunits, such as Kir4.1, Kir4.2, or potentially Kir2.1 (Xu et al., 2000). This is largely based on the data showing that heterologous expression of Kir5.1 channels alone does not yield a functional channel (Pessia et al., 1996). If this is the case, the additional co-expressed subunit in 5-HT neurons, which would bind to Kir5.1 subunits and form a functional heteromeric channel is not clear, as we found: (1) no Kir4.1 protein expression in 5-HT neurons, (2) little or no Kir4.2 gene expression in the RNA sequencing data, and (3) no evidence of Kir4.2 protein in the brainstem with immunohistochemistry (data not shown). Kir2.1 gene expression has been localized to neurons, but we and others (Okaty et al., 2015) found little or no evidence of Kir2.1 gene expression in 5-HT neurons. Although, we cannot identify a specific binding partner that would lead to a functional Kir5.1-containing channel, the addition of this subunit may be required for a Kir channel to exhibit enhanced pH sensitivity with a pKa in the physiologic range. But, these observations also lead us to question if indeed homomeric Kir5.1 channels, under certain circumstances, could functionally contribute to resting membrane potential and/or alterations in membrane excitability. Functional homomeric Kir5.1 channels have been demonstrated when they are co-expressed with the post-synaptic density 95 (PSD-95) protein, which directs Kir5.1 channel expression to the plasma membrane (Tanemoto et al., 2002). This membrane localization is reversibly prevented by activation of protein kinase A, suggesting regulation of the Kir5.1 channel localization may modulate its function. Other data indicate functional Kir5.1 channel homomers in the retina do not require colocalization with Kir4.1, Kir4.2, or PSD-95 but may require an additional factor for independent functionality (Ishii et al., 2003). Our confocal imaging of brainstem 5-HT neurons suggested that most Kir5.1 channel protein was perinuclear with only a fraction at or near the plasma membrane, similar to previous descriptions (Tanemoto et al., 2002). Further investigation is required to determine if this sub-cellular localization is the "default" state of the channel, and only under certain circumstances is it localized to the plasma membrane to alter cellular function.
Lastly, we have chosen to focus on the validation of a small number of identified ion channels that were nominated by the RNA sequencing data, although there were many more genes identified. Our focus on K + channels was based on the long-standing assumption that pH sensitive currents are most likely due to the closure of "leak" K + channels (Putnam et al., 2004). However, equally plausible is the concept of the activation of a cation current by hypercapnic acidosis leading to depolarization of respiratory chemoreceptors, which has previously been hypothesized for chemosensitive brainstem 5-HT neurons (Richerson, 2004). Moreover, we further limited our analyses to genes that encode or modulate channels with known pH-sensitivity, where there could be other novel candidates that could be functionally altered through various mechanisms (alternative splicing, post-translational modifications, etc.), or others prominently expressed in other organs that are repurposed in neurons. For example, one of the developmentally regulated genes in our sequencing data is Scnn1a, which encodes the alpha-ENaC a sodium channel that regulates electrolyte transport across membranes particularly in the kidney. However, alpha-ENaC activity is highly pHdependent, and was recently found to be co-localized to 5-HT neurons in the brainstem. Finally, it is quite possible that cellular pH sensing mechanisms require multiple pH-sensitive ion channels/proteins, including complementary expression of pH-insensitive ion channels. This concept has a great deal of merit, particularly in light of the requirements of combined expression of two pH-sensitive proteins TASK-2 and Gpr4 for direct pH sensing in chemosensitive RTN neurons (Kumar et al., 2015), and other data showing Nalcn, KCNQ, HCN, THIK-1 and others altered chemosensory responses in RTN chemoreceptors independent of a role in direct pH sensing (Lazarenko et al., 2010;Hawryluk et al., 2012;Hawkins et al., 2015;Shi et al., 2016). Thus, additional changes in gene expression of pH-insensitive ion channels could also contribute to pH sensing mechanisms, but may not be identified using the approach herein. Despite the stated limitations, this cell-specific gene expression-based approach could be ideal to identify potential panels of pH sensitive proteins (and those that have indirect, supportive roles) by comparing multiple cell populations, which would be difficult to do with the standard electrophysiology-based approach using pharmacological inhibitors. Moreover, a similar approach could be used to identify additional acidsensing mechanisms in other neural systems, including pain sensation/perception and fear/anxiety disorders. While the overall significance of these observations remains unclear (since the molecular underpinnings of pH sensitivity in brainstem 5-HT neurons remains unclear), further exploration of these and additional candidate genes and their functional validation is warranted.

AUTHOR CONTRIBUTIONS
MP and GM performed experiments, analyzed data, performed statistical analyses and contributed to manuscript preparation. PL analyzed sequencing data and performed statistical analyses. MH performed experiments and contributed to manuscript preparation.