Transcriptome of neonatal preBötzinger complex neurones in Dbx1 reporter mice

We sequenced the transcriptome of brainstem interneurons in the specialized respiratory rhythmogenic site dubbed preBötzinger Complex (preBötC) from newborn mice. To distinguish molecular characteristics of the core oscillator we compared preBötC neurons derived from Dbx1-expressing progenitors that are respiratory rhythmogenic to neighbouring non-Dbx1-derived neurons, which support other respiratory and non-respiratory functions. Results in three categories are particularly salient. First, Dbx1 preBötC neurons express κ-opioid receptors in addition to μ-opioid receptors that heretofore have been associated with opiate respiratory depression, which may have clinical applications. Second, Dbx1 preBötC neurons express the hypoxia-inducible transcription factor Hif1a at levels three-times higher than non-Dbx1 neurons, which links core rhythmogenic microcircuits to O2-related chemosensation for the first time. Third, we detected a suite of transcription factors including Hoxa4 whose expression pattern may define the rostral preBötC border, Pbx3 that may influence ipsilateral connectivity, and Pax8 that may pertain to a ventrally-derived subset of Dbx1 preBötC neurons. These data establish the transcriptomic signature of the core respiratory oscillator at a perinatal stage of development.

Neural rhythms that drive inspiratory breathing movements in mammals originate from the brainstem preBötzinger complex (preBötC) 1,2 . Neurons derived from Dbx1-expressing progenitors comprise its rhythmogenic core [3][4][5][6][7][8][9] . Although we know the site and neuronal constituents at the point of origin of respiratory rhythm, the cellular and molecular mechanisms that generate and control respiration remain incompletely understood.
Electrophysiological recordings in preBötC neurons generally, and Dbx1-derived preBötC neurons in particular, have characterized intrinsic membrane properties, including ion channels, membrane pumps and transporters, as well as synaptic currents that influence the neural mechanisms of respiration 2,10,11 . However, testing their relative rhythm-and pattern-generating roles typically relies on promiscuous pharmacology and leads to inconclusive results. We argue that identifying specific subunits, isoforms, and genes that underlie putatively rhythmogenic conductances and integral membrane proteins would facilitate more conclusive experiments. Knowledge of the newborn mouse preBötC transcriptome -the expressed transcripts and their relative quantity -could be exploited to develop targeted physiological experiments, with the added benefit of uncovering novel genes that may influence preBötC development as well as regulate respiratory function.
Here we provide the first RNA-Seq gene expression profile for preBötC neurons in newborn mice. We analysed gene expression levels within the Dbx1 population as well as differential expression between Dbx1 and non-Dbx1-expressing populations, and we interpret their significance for defining the structure and function of the preBötC in the context of existing literature. These data are publicly available in an open access database (NCBI gene expression omnibus, https://www.ncbi.nlm.nih.gov/geo/) for custom analyses and applications that interrogate preBötC development as well as the cellular and molecular neural bases for breathing behaviour. Amino acid neurotransmitters. Excitatory preBötC neurons are respiratory rhythmogenic [19][20][21][22][23][24][25] . Inhibitory neurons, which populate the preBötC in roughly equal numbers, regulate respiratory rhythm and mediate sensorimotor integration [26][27][28][29] (q.v., refs 11, 30).
Given the well-established role of Dbx1 preBötC neurons in rhythm generation 3, 4, 6-9 , we expected the Dbx1 samples to express transcripts associated with excitatory transmitter phenotype. By contrast, we expected non-Dbx1 samples to exhibit markers of both excitatory and inhibitory transmitter phenotypes.

Short non-coding RNAs
Contrary to our expectations, gene expression for inhibitory amino acid-synthesizing enzymes and transporters was commensurate for Dbx1 and non-Dbx1 neurons ( Supplementary Fig. S1b  The RNA related to inhibitory synaptic transmission may remain untranslated in most Dbx1 preBötC neurons similar to excitatory CA1 hippocampal neurons where post-transcriptional regulation controls neurotransmitter phenotype 34 . Nevertheless, a non-negligible subset of Dbx1 preBötC neurons communicate via chloride-mediated synaptic inhibition. For example, Gray et al. (ref. 4) identified inhibitory Dbx1 preBötC neurons that expressed GAD1 or GlyT2. Dbx1 preBötC neurons with inhibitory transmitter phenotype would be well equipped to periodically suppress activity in respiratory nuclei such as the expiratory-related lateral parafacial respiratory group [35][36][37][38] or the post-inspiratory complex 39 , both of which discharge neural bursts out-of-sync with the preBötC. It is also possible that inhibitory Dbx1 preBötC neurons suppress orofacial behaviours during the inspiratory phase of the breathing cycle 40,41 or participate in inhibitory pulmonary stretch receptor feedback.
The endogenous ligands for neurokinin 1 receptors, neurokinin A and substance P, both result from tachykinin precursor 1 (Tac1), which is highly expressed in both sample populations (RPKM > 60, Supplementary  Fig. S2) but Tac1 did not meet the criterion for differential expression (L2FC = 1.36, p = 0.018, FDR = 0.48). This result is not surprising given that substance P is a cotransmitter in overlapping populations that modulate preBötC function 46,52,53 .
Transcription factors and other distinguishing markers of preBötC Dbx1 and non-Dbx1 neurons. Combinatorial codes of transcription factors, expressed in the developing spinal cord and hindbrain, govern the assembly of central pattern generating circuits for locomotion and breathing [65][66][67][68] . We cross-referenced our RNA-Seq results with a list of murine transcription factors (RIKEN Institute, Wako, Japan) as well as the curated TF-checkpoint project 69    HIF-1 (hypoxia-inducible factor 1, a protein complex) regulates gene transcription cascades in response to hypoxia [70][71][72] . Its expression increases within minutes under hypoxic conditions 73 . We constantly perfused oxygenated solution through the tissue as we isolated neurons from the superficial 40 µm of the tissue slices (e.g., Fig. 1b) wherein there is no diminution of the oxygen gradient 74,75 . Hif1a, differentially-expressed at higher levels in Dbx1 neurons, encodes the O 2 -sensing component of the complex (the HIF-1α subunit). In heterozygotic Hif1a knockouts, the carotid body, the primary sensor of arterial O 2 , no longer reacts to hypoxia, which is lethal 76 . To our knowledge HIF-1α has not been associated with respiratory rhythm generation directly, although its upregulation in brainstem may occur at high-altitudes during hypoxic conditions 77 . HIF-1α being highly-and differentially expressed in core rhythmogenic neurons suggests that Dbx1 preBötC neurons may be programmed to increase ventilation in response to hypoxia.
HIF-1α is widely associated with angiogenesis 78 . preBötC neurons are embedded within a network of arterioles 79 like chemosensitive retrotrapezoid neurons 80 , which suggests that preBötC neurons, like their retrotrapezoid counterparts, need access to blood for chemosensation.
Pbx3 is a Hox gene cofactor 81 . The expression pattern for Pbx3 is mosaic in transverse sections of the newborn mouse medulla at the level of the preBötC 82 as well as in parasagittal sections from the Allen Institute developing mouse brain atlas (AIDMBA) (Fig. 5a, AIDMBA: P4, Pbx3, slide 10). Dbx1 and Pbx3 knockout-mice both die perinatally due to central hypoventilation 3,4,82 . Dbx1 knockout mice never breathe and form no preBötC 3,4 . In contrast, the Pbx3 knockout mice breathe irregularly with periodic apnoea, which indicates that the preBötC is present but dysfunctional. Pbx3 was shown to interact with other Pbx genes and Hox genes to influence connectivity of motor neurones in ipsilateral motor pools 83 . We posit that Pbx3 may influence ipsilateral connectivity in the preBötC. Commissural connectivity, in contrast, is governed by Robo3 in Dbx1 preBötC neurons 3 , which was ostensibly undetected (Dbx1 samples, RPKM = 0.06 ± 0.06; non-Dbx1 samples, RPKM = 0.01 ± 0.02) most likely because commissural projections are mature by early post-natal development. Failure to properly interconnect locally could impair emergent network rhythmicity 6,24 , and impede normal respiration in Pbx3-deficient mice.
Homeobox (Hox) genes govern cell fate identity along the anterior-posterior axis 91 . We observed a cascade of expression among the Hox2-5 genes across the Hoxa, Hoxb, Hoxc, and Hoxd chromosomal clusters, which peaked with Hoxa4 (Fig. 4). Hoxa4 showed the second highest ∆RPKM among all transcription factors (Fig. 2d, lower plot, and Fig. 4), but it did not pass the threshold for differential expression (L2FC = 1.39, p = 0.0022, FDR = 0.16). In the anterior-posterior axis of the lower medulla, Hoxa4 expression stops at the caudal border of the compact division of the nucleus ambiguus, which coincides with the rostral limit of the preBötC (Fig. 5b, reproduced with permission from ref. 92), and thus Hoxa4 expression might influence the rostral preBötC boundary, which would apply to both Dbx1 and non-Dbx1 neurons. Hox genes that are also relevant to respiration, in the cervical spinal cord and not the brainstem, include Hoxa5 and Hoxc5 93 , which influence phrenic motoneuron development and survival. That these genes continue to be expressed postnatally to maintain the organization of phrenic motor columns may be generalizable to postnatal expression of Hox genes in the preBötC postnatally, but that remains to be investigated. Paired box (Pax) genes influence cell lineage specification 94 . Pax8 was the highest expressed transcription factor of its class (Fig. 4) in both Dbx1 and non-Dbx1 neurons (L2FC = 0.17, p = 0.79, FDR = 1.0), which is notable because Pax8 expression perdures in ventral hindbrain neurons in adult mice 95 . Data from the AIDMBA show Pax8 expression clustered in the vicinity of the preBötC at P4 (Fig. 5c; Pax8, slide 11) and Pax8 expression is on the ventral edge of Dbx1 expression at E11.5 ( Fig. 5d and e, AIDMBA: Dbx1, slide 9, and Pax8, slide 9). Pax7 was completely undetected (Fig. 4). Pax2 had the highest L2FC among Pax genes but did not meet the criterion for differential expression (L2FC = 1.39, p = 0.0067, FDR = 0.30). Even though Pax8 and Pax2 were not differentially expressed, these data are consistent with core preBötC rhythmogenic neurons making up the ventral subset of V0 interneurons (i.e., V0 V ), which express Dbx1 as well as Pax8 and Pax2, but not Pax7. In contrast, Dbx1-derived dorsal interneurons (V0 D ), which may reside in the reticular formation and not the preBötC, express Pax7 3, 96-98 .
Dbx1 was undetected in preBötC neurons (Fig. 4) because its expression is restricted to embryonic development (E9.5-E11.5) 3, 4, 16, 96 . Two of the Dbx1 samples (Pos2 and Pos3) exhibited Evx1 expression for a mean RPKM of 5.76. The third Dbx1 sample (Pos1) did not show Evx1 expression. Evx1 was not detected in any non-Dbx1 samples. Because four of the six samples did contain detectable Evx1, DESeq. 2 cannot compute statistics beyond raw reads. Nevertheless, we suspect that a significant subset of Dbx1-expressing progenitor cells differentiate into Evx1-expressing ventral (V0 V ) interneurons, which ultimately form the preBötC core. In the lumbar spinal cord, Evx1 expression is critical for V0 V interneuron identity 99 . Commissural V0 V interneurons in zebrafish spinal cord require Evx1 (and Evx2) to become glutamatergic 100 . We posit that Evx1 may play analogous roles in Dbx1-derived interneurons of the preBötC.
Neither Dbx1 nor non-Dbx1 preBötC neurons express Atoh1 (Fig. 4), which is an embryonic transcription factor associated with progenitors of the central chemoreceptive ventral parafacial (pF V ) or retrotrapezoid nucleus 43, 106-108 . We observed low expression of Tlx3 and Egr2 (a.k.a., Krox20 [ref. 109]). The combined RPKM of Tlx3 and Egr2 measured 0.78 ± 1.20. Those data are not surprising because Tlx3 and Egr2 are associated with the parafacial and serotonergic raphé neurons, but not preBötC neurons. Tshz3 and Mecp2, which are linked to respiratory-related dysfunction, and the latter specifically with Rett syndrome [110][111][112] , showed modest expression across both Dbx1 and non-Dbx1 neurons. Their combined RPKM measured 3.33 ± 1.62. A notable gene, which is not a transcription factor, but was significantly differentially expressed in Dbx1 samples was synaptotagmin-10 (Syt10, L2FC Syt10 = 2.34, p = 0.00034, FDR = 0.047), which is involved in synaptic exocytosis 113,114 , and could be involved in synaptic depression observed in preBötC rhythmogenic neurons 24, 115 . Neonatal preBötC transcriptome provides a baseline for comparative analyses during development. We present this transcriptome, including non-coding transcripts, from perinatal Dbx1 preBötC neurons, which are obligatory rhythm-and pattern-generating interneurons for breathing. We also present the transcriptome of non-Dbx1 preBötC neurons, which serve non-rhythmogenic respiratory and non-respiratory functions. We analysed the perinatal transcriptome because animals at this age are widely employed in respiratory neurobiology research, and provide testable predictions for experiments in adult rodents 2,11,116,117 . Further, brainstem tissue explants from neonatal fluorescent reporter mice are optimal for visual identification and isolation of single cells. These data can be exploited or meta-analysed to design new experiments and approaches that manipulate Dbx1 and non-Dbx1 neuron development and physiology and thus elucidate their roles in the neural generation and control of breathing.

Methods
Animals. All procedures were approved by the Institutional Animal Care and Use Committee at the College of William and Mary, which conforms to the guidelines of the US Public Health Service policy on humane care and use of laboratory animals (Office of Laboratory Animal Welfare, the National Institutes of Health, Bethesda, MD). To identify Dbx1-derived preBötC neurons, we crossed female tamoxifen-inducible Dbx1 Cre-driver mice (Dbx1 CreERT2 ; Hirata et al., 2009; stock no. 028131, The Jackson Laboratory, Bar Harbor, ME) with male reporter mice whose Rosa26 locus was modified by targeted insertion of a loxP-flanked STOP cassette followed by a gene for the fluorescent protein tdTomato (Rosa26 tdTomato , stock no. 007905; The Jackson Laboratory). Tamoxifen (22.5 mg/kg body mass) was administered to pregnant dams during gestation at embryonic day 9.5 (ref. 16). In their offspring, Dbx1 CreERT2 ; Rosa26 tdTomato mice, Cre-mediated recombination resulted in cytosolic tdTomato expression in cells derived from Dbx1-expressing progenitors.
Medullary slices. Neonatal Dbx1 CreERT2 ; Rosa26 tdTomato mice (postnatal day 2) of both sexes were anesthetized by hypothermia and their neuraxes were dissected in ice-cold artificial cerebrospinal fluid (aCSF) containing the following (in mM): 124 NaCl, 3 KCl, 1.5 CaCl 2 , 1 MgSO 4 , 25 NaHCO 3 , 0.5 NaH 2 PO 4 , and 30 dextrose equilibrated with 95% O 2 /5% CO 2 , pH 7.4. We prepared transverse brainstem sections (500 µm thick) whose rostral surface exposed the preBötC, which was verified by its position between the semi-compact division of the nucleus ambiguus and the dorsal boundary of the principal loop of the inferior olive 118 (Fig. 1a). Slices were perfused with ice-cold aCSF bubbled with 95% O 2 /5% CO 2 at 5 ml/min. Neuron isolation. Dbx1 preBötC neurons were identified using epifluorescence and removed under bright-field imaging (specifically, a version of differential interference contrast microscopy called 'Dodt' by Zeiss) on a fixed-stage microscope. Non-Dbx1 neurons were identified as neuronal somata without tdTomato fluorescence. We fabricated and heat-sterilized micropipettes from borosilicate capillary glass (1.50 mm outer diameter, 0.86 mm inner diameter). Micropipettes were devoid of solution prior to the isolating cells. After forming a seal with the plasma membrane, we applied negative pressure to draw single neurons into the tip of the micropipette (Fig. 1b). Tips containing the sampled neurons were immediately broken in sterile RNase/DNase-free tubes submersed in liquid nitrogen until a total of 15 neurons were collected. Each sample of 15 neurons was collected from a different animal (Fig. 1c, n = 6). We also performed mock cellular isolations as a control, bringing micropipettes into contact with the slice surface but without collecting cells. In mock cell isolations, we drew the same amount of fluid as real experiments, and performed cDNA synthesis in exactly the same manner as the real samples (n = 3). cDNA synthesis. Nuclease-free water and lysis buffer with RNase inhibitor, from the SMART-Seq v4 ultra low input RNA kit for Sequencing (634889, Clontech, Mountain View, CA), were added to each tube to raise the volume to 10.5 μL. The contents were incubated with sonication for 5 min to lyse the neurons and release cytoplasmic RNA, and then transferred to a 0.2 mL RNase-free PCR tube. First strand cDNA was synthesized by performing reverse transcription in a thermal cycler (42 °C for 90 min, 70 °C for 10 min). Then these cDNAs were primed with the 3′ SMART-Seq CDS Primer IIA, and we used the SMART-Seq v4 Oligonucleotide for template switching at the 5′ end of the transcript. The cDNA was then amplified by LD-PCR from the SMART sequences introduced by 3′ SMART-Seq CDS Primer IIA and the SMART-Seq v4 Oligonucleotide in a heated thermal cycler (95 °C for 1 min, 17 cycles of 98 °C for 10 s, 65 °C for 30 s, 68 °C for 3 min; 72 °C for 10 min). PCR-amplified cDNA was purified by immobilization on Agencourt AMPure XP beads (A63880, Beckman Coulter, Brea, CA), which were then washed with 80% ethanol and cDNA was eluted with elution buffer. Amplified cDNA was validated using the Agilent 2100 Bioanalyzer and Agilent's High Sensitivity Kit (5067-4626, Agilent Technologies, Santa Clara, CA), and its concentration was determined using Qubit dsDNA High-sensitivity Assay Kit (Molecular Probes). The full-length cDNA output was processed with the Nextera Library Preparation Kits (FC-131-1024, Illumina, San Diego, CA) to obtain cDNA libraries for RNA-Seq experiments. Dbx1 and non-Dbx1 samples contained an average of 1481 ± 352 pg/µl of amplified cDNA, whereas mock cell-isolation samples contained 93 pg/ µl cDNA (n = 3). Thus our transcriptome reflects cytoplasmic RNA from preBötC neurons. Figure 1d illustrates cDNA synthesis and quality control steps. Sequence analysis. The six cDNA libraries were submitted to the DNA Sequencing Center at Brigham Young University (Provo, UT) for sequencing on an Illumina -HiSeq. 2500. We received an average of 15,310,365 single-end reads per sample, with an average read length of 49.5 bp (range: 35-50 bp). We received these sequences and quality scores in the form of FASTQ files (first part of Fig. 1e) and aligned them to the mm10 murine genomic database from the University of California at Santa Cruz (http://hgdownload.soe.ucsc.edu/ downloads.html) using Tophat2 software 119 with the parameters specified in "TophatParameters.txt", running on the Galaxy cluster computer at Johns Hopkins University (https://usegalaxy.org/). Average mapping rate was 85%. We sorted the binary alignment/map (BAM) files using Samtools (http://www.htslib.org/download/), utilities that process short DNA sequence read alignments, and then we applied HTSeq-count 120 to map reads to genes with the following general command: htseq-count -f ba×m -s no -i gene_name BAMFILE.bam genes.gtf> BAMFILE.txt where genes.gtf is the gene annotation from the Ensembl (described below). This resulted in an average of 21,131,831 aligned reads across the six samples of which 7,815,666 uniquely mapped to genes. We implemented custom Python scripts to compute reads per kilobase of transcript per million mapped reads (RPKM) for each gene in each sample by normalizing for exon length according to the Ensembl mouse gene annotation database. These scripts read the Ensembl mouse gene annotation database (http://www.ensembl.org/Mus_musculus/Info/ Index), and then apply the gtf-to-genes−1.40 package (https://pypi.python.org/pypi/gtf_to_genes) to extract the total exon length for each gene. We analysed transcription factors by cross-referencing our transcriptome data to the Riken Institute (Wako, Saitama, Japan) list of transcription factors (http://genome.gsc.riken.jp/TFdb/tf_list. html) and the TF-checkpoint database (http://www.tfcheckpoint.org/data/TFCheckpoint_download_180515.txt) from the Norwegian University of Science and Technology (Trondheim, Norway). We evaluated differential gene expression between Dbx1 and non-Dbx1 samples via the DESeq. 2 121 algorithms performed on the HTSeq-count output (Fig. 1e).
We wrote more than 20 custom Python/R scripts to process the data displayed in Figs. 2-4 and Supplementary Figs S1-S9; and to process the DESeq. 2 results. These start from the BAM files generated from running the original FASTQ files on Tophat with the default parameters using the UCSC mm10 standard murine genome. They may be run sequentially from the beginning using the "runAll.sh" shell script that has been provided. The bioinformatics procedures are illustrated schematically in Fig. 1e. These scripts are freely available and open source (http://dbx1seq.sourceforge.net/). Differential expression. We evaluated differential gene expression between Dbx1 and non-Dbx1 neurons using DESeq. 2 (ref. 121), which computes the probability of obtaining the observed mean difference in gene expression if the Dbx1 and non-Dbx1 samples were drawn from the same underlying genetically homogenous population (the null hypothesis) as well as an adjusted probability, which reflects the false discovery rate (FDR) associated with multiple comparisons (i.e., ~50,000 genes). Significance level was set at FDR = 0.1 by convention. DESeq. 2 modifies the dataset to account for heteroscedasticity intrinsic to RNA-Seq data when calculating the raw p-value from a Wald test. DESeq. 2 calculates FDR 32 using a Benjamini-Hochberg correction 31 , which aggressively combats Type I errors, i.e., false positive discoveries 31, 32 at the expense of inflating false negatives (Type II errors) 33 . Acknowledging that p-value thresholds can be misleading if experiments are not analysed in context 122 , we interpret our data according to raw and adjusted probability (i.e., FDR) for genes that prior and corroborating literature suggests have respiratory neurobiological relevance. The associated log 2 fold change (L2FC) quantifies the degree of differential expression. L2FC > 0 for genes more highly expressed in the Dbx1 samples; L2FC < 0 for genes more highly expressed in the non-Dbx1 samples.
Data availability. The original data, which includes FASTQ files (raw nucleotide sequences and quality scores) and processed data files (sequenced reads that have been aligned and normalized to a murine reference genome) are publicly available in the NCBI Gene Expression Omnibus database, https://www.ncbi.nlm.nih.gov/ geo/, accession number GSE100356.