Unless otherwise stated, all experiments were performed on outbred CD-1 mice (Charles River, Germany). Mice in which the Ptger3 gene (EP3R) was selectively deleted had a C57BL/6J background (Fabre et al., 2001). All mice were reared by their mothers under standardized conditions with a 12:12-h light-dark cycle. Food and water were provided ad libitum. The studies were performed in accordance with European Community Guidelines and approved by the regional Ethics Committee. The animals were reared and kept at the Department of Comparative Medicine, Karolinska Institutet, Stockholm, Sweden.

Mice were sedated on postnatal day (P) 2–4 using 5% Vetflurane (Virbac Animal Health Ltd. Bury St. Edmunds, Suffolk, United Kingdom) dispersed by an anesthesia unit (AgnTho’s AB, Univentor 410). After decapitation, tissue was placed in ice-cold artificial cerebrospinal fluid (ACSF) containing (in mM): 25 D-glucose, 2.5 KCl, 1.25 NaH₂PO₄, 12 N-acetyl–L-cystein, 147 glycerol, 3 ethyl pyruvate, 20 HEPES (pH 7.4 at 4°C), 5 sodium ascorbate, 30 NaHCO₃, 10 MgSO₄, 0.5 CaCl₂. Final osmolarity was 300–320 mOsm after bubbling with 95% O₂/5% CO₂. The brainstem was removed and attached (cyanoacrylate glue) dorsal-side down to an angled agar block (tilted 20° from vertical). The agar block was then attached to the stage of a vibrating microtome (7000smz-2, Campden Instruments, Loughborough, United Kingdom). Serial transverse slices of 200–400 μm thickness were taken starting from the rostral end of the brainstem until the opening of the 4th ventricle was reached. Subsequent trim slices (100–200 μm thickness) were cut at low speeds (<50 μm/s advance) until landmarks associated with the preBötC could be seen (Ruangkittisakul et al., 2014). For in vitro recordings, rhythmic slices were cut 650 μm thick, with the preBötC exposed on the rostral surface of the slice. For single-cell RNA-Seq samples, slices were cut 300 μm thick, with the preBötC centered in the slice.

Transverse slices were transferred to a recording chamber (BSC1, Scientific Systems Design Inc., Mississauga, Canada) and perfused with recording ACSF (bubbled with 95% O₂/5% CO₂, 26°C, 300–320 mOsm, pH 7.3–7.4). Recording ACSF, containing (in mM): 25 D-glucose, 9 KCl, 1.25 NaH₂PO₄, 1 sodium ascorbate, 26 NaHCO₃, 115 NaCl, 1 MgCl₂, 1.5 CaCl₂. A suction electrode pipette with ca. 70–100 μm opening (GC120-15, borosilicate glass capillaries, 1.2 mm O.D. × 0.69 mm I.D., Clark Electromedical Instruments, Reading, United Kingdom) was gently attached to the hypoglossal (XII) rootlets, which were retained in the transverse slice. The extracellular recording signal was amplified with a differential AC Amplifier (Model 1700, A-M Systems, Sequim, United States). The signal was filtered (low-cut 300 Hz, high-cut 1,000 Hz, Gain 10,000×) and recorded with a digitizer (Digidata 1320A, Molecular Devices, San Jose, CA, United States) at 20,000 samples/s.

Slices were perfused with ACSF for 20–40 min prior to recording a 10-min baseline. Drug stocks were added to recording ringer and sequentially washed into the chamber using an automated pinch valve system (PC-16 controller, Bioscience Tools, San Diego, CA, United States). Slices were exposed to drugs for 20 min before moving on to the next drug. However, only the last 10 min of each exposure epoch were analyzed to allow for stabilization of the rhythm. The following drugs were used: EP1 antagonist GW848867X (Cayman Chemicals, Ann Arbor, MI, United States, 50 nM), EP2 Antagonist PF-04418948 (Cayman Chemical, 500 nM), EP3 antagonist DG041 (Tocris, Bristol, United Kingdom, 40 nM), EP4 antagonist MK2894 (APExBIO, Houston, TX, United States, 100 nM), Prostaglandin E2 (Sigma-Aldrich, St. Louis, MO, United States, 1 nM-1 μM), Prostaglandin F receptor antagonist AL-8810 (Cayman Chemical, 1 μM), Prostaglandin D2 receptor antagonist Setipiprant (Cayman Chemical, 500 nM), Non-selective EP2/3/4 agonist 11-deoxy PGE1 (Cayman Chemical, 10 nM), EP2R agonist Butaprost (Cayman Chemical, 5 μM), EP3 receptor agonist Sulprostone (Cayman Chemical, 1 μM).

Raw data from XII rootlet recordings was analyzed using custom software written in Python. In short, XII recordings were resampled at 1/10 rate and artifacts were excluded using a Hampel filter with a 0.5 s window and a threshold of 5 standard deviations. Bursts were defined as peaks of minimum prominence of 0.2 and distance of 15 s. Results were manually checked for false positives. XII burst area was implemented as a leaky integrator with a decay constant of 0.98. Burst characteristics were calculated for each burst individually and subsequently averaged for each sample. If percentages are given, values were normalized to the mean of the baseline period for the given variable. Two-sided Smirnov-Grubb’s test was performed to exclude outliers with p < 0.05. Multimodality was tested using the Hartigan’s Dip test (Hartigan and Hartigan, 1985).

Statistics for rootlet recordings are the result of General Estimating Equations (GEE) performed with a custom R script using the “geepack” and “emmeans” packages. GEE was chosen as an alternative to an ANOVA, since the assumption of homogeneity of variances was not satisfied. This can be explained by the increase in variability after addition of PGE2, leading to different variances between the experimental groups. In short, a Gaussian GEE model was fitted to the data set and subsequently Bonferroni adjusted contrasts were calculated, if permitted by the results of the GEE model.

Outbred CD-1 littermates (age P4) were dissected on the same day, using the methods described in Slice Preparation. All procedures were carried out under aseptic conditions. We cut slices of 300 μm thickness, with the preBötC approximately centered in the slice according to anatomical landmarks (Ruangkittisakul et al., 2014). Using a sterilized metal tool (NIH Style Neuro Punches, Agnthos, Lidingö, Sweden), circular punchouts (ø 690 μm) were cut from the ventral-lateral border on both sides of each slice, centered over the approximate location of the preBötC in the transverse plane. After extraction, the punchouts were immediately placed into dissociation ACSF containing (in mM): 25 D-glucose, 2.5 KCl, 1.25 NaH₂PO₄, 12 N-acetyl-L-cysteine, 3 ethyl pyruvate, 20 HEPES, 5 sodium ascorbate, 30 NaHCO₃, 10 MgSO₄, 0.5 CaCl₂, and 0.1 U/μL RNase inhibitor (RNasin Plus, Promega, Madison, WI, United States). In addition, we added 5% w/v D-trehalose dihydrate (VWR Chemicals, Radnor, PA, United States) to improve cell viability (Saxena et al., 2012). The final osmolarity of the dissociation ACSF was 300–320 mOsm and pH was 7.3–7.4 after bubbling with 95% O₂/5% CO₂. We enzymatically dissociated the cells by adding 1 mg/mL pronase (Roche, Basel, Switzerland) to the dissociation buffer. Dissociation was carried out in silanized glass conical tubes at room temperature (21–22°C) for 25 min with periodic agitation of the solution. During this time, the head space of the conical tube was continuously gassed with 95% O₂/5% CO₂. We then titrated the tissue with silanized glass Pasteur pipettes, which had been fire polished to an inner diameter of 200, 400, or 800 μm. Further protease action was then inhibited by addition of 1% fetal bovine serum (FBS) to the dissociation buffer. Additionally, we added a fixable fluorescent viability indicator to stain dead cells (Live-or-Dye NucFix Red, Biotium, Fremont, CA, United States). Cells were incubated on ice in the FBS/dye-containing solution for 5 min and then centrifuged at 300 × g for 5 min. The cell pellet was resuspended in ice-cold dissociation ACSF containing 2% FBS that had been pre-equilibrated with 95% O₂/5% CO₂. Finally, we cryopreserved the cells in methanol (Alles et al., 2017). In brief, methanol chilled to −20°C was added dropwise to the resuspended cells (final concentration 80% methanol). Methanol-fixed samples were allowed to incubate at 4°C for at least 15 min before being transferred to a −80°C freezer for storage. Samples were stored for less than 1 week before subsequent library preparation.

Single-cell transcriptomic libraries were generated with the STRT-seq-2i method (Hochgerner et al., 2017). First, methanol-fixed samples were pooled into a single 15 mL conical tube and centrifuged at 2,000 × g for 10 min. Cells were resuspended in a solution of ice-cold PBS without Ca²⁺or Mg²⁺, including 20 μg/mL Hoechst 33342 (Merck), 2% bovine serum albumin (Merck, Kenilworth, NJ, United States) and 0.5 U/μL RNase inhibitor (RNasin Plus, Promega). Cells were washed a second time, counted, and then diluted to a concentration of 20 cells/μL and less than 0.1% BSA. The cells were dispensed via limiting dilution (iCell8 MSND, Wafergen, Fremont, CA, United States) into a 9600 microwell array (Wafergen). The chip was imaged on a microscope equipped with epifluorescence and screened via software detection (CellSelect, Wafergen) for empty wells, doublets (via Hoechst fluorescence), and dead cells (via red viability stain). Reagents for library preparation were only dispensed to wells called as viable singlets. All remaining library preparation, quality control, and tagmentation steps were performed as described in Hochgerner et al. (2017) Final libraries were sequenced on an Illumina HiSeq 2500, 50 bp single-end reads. Raw sequencing data was pre-processed as previously described (Hochgerner et al., 2017). In brief, UMIs containing any bases of Phred score <17 were rejected. Reads with less than 3 Gs following the 6 bp UMI were filtered out. After trimming poly-A tails, only reads with at least 25 remaining bases were kept. Filtered reads were mapped and demultiplexed with zUMIs (Parekh et al., 2018). Detected barcodes were filtered according to the barcode combinations dispensed during library tagmentation. Total reads per cell were adaptively downsampled to within three times the median absolute deviation.

The raw count matrix was further processed and analyzed using “Scater” and “Seurat” libraries in R (McCarthy et al., 2017; Hao et al., 2021). Cells expressing less than 1,000 genes, and genes expressed in fewer than 10 cells were excluded from analysis. Of the remaining cells, we eliminated outliers with log-library sizes and percent mitochondrial reads exceeding three times the median absolute deviation. Counts were normalized and scaled using “SCTransform” (Hafemeister and Satija, 2019). We removed doublets using “DoubletFinder” (McGinnis et al., 2019). Finally, the data was dimensionally reduced via PCA, and the first 20 PCs were used to cluster cells via the graph-based Leiden algorithm (Traag et al., 2019). The data was also dimensionally reduced to a 2D visualization via uniform manifold approximation and projection (UMAP) using the first 20 PCs (McInnes et al., 2020). We tested for differential expression with a non-parametric Wilcoxon rank sum test. Resulting p-values were adjusted with a Bonferroni correction, using the total number of genes in the count matrix. We tested DE only on genes which were expressed by at least 10% of the cells in a cluster and with log2 fold-change between groups greater than 0.1. We compared each cluster to all other cells in the data (“FindAllMarkers” function in “Seurat”). Non-neuronal clusters were annotated by differential expression for astrocytes, microglia, endothelial cells, oligodendrocytes, and oligodendrocyte precursor cells. We also compared differential expression in each of the remaining clusters with respect to all cells in the 5 non-neuronal clusters. Differential expression values are reported as the log2 fold-change in average expression between the comparison groups, followed by the percentage of cells expressing transcripts in the first group versus percentage of cells in the comparison group expressing gene transcripts.

We measured the dose response of PGE2 on inspiratory-related motor activity via extracellular recordings taken from hypoglossal nerve rootlets in rhythmically active transverse slice preparations containing the preBötC (Figure 1A). Lower concentrations of PGE2 (1–10 nM) increased the mean duration of the inspiratory burst period relative to baseline values (+21% above baseline in 1 nM PGE2, +25% above baseline in 10 nM PGE2), which was significantly different from the mean burst period measured in untreated control slices (Figures 1B–C). Additionally, the duration of the burst period was more irregular in 10 nM PGE2 compared to control slices (Figure 1G). In 1 μM PGE2, the mean duration of the burst period decreased compared to baseline values (−24% below baseline in 1 μM PGE2), which was also significantly different from untreated control slices (Figure 1C).

In all the recordings from PGE2-treated and control slices, burst amplitude and burst width gradually decreased over the duration of the recording (Figures 1D,E). However, this rundown was consistent and not statistically different between control and treatment groups. Similarly, there was no difference in the XII burst area (Figure 1F). These results suggest that PGE2 primarily influences the timing and regularity of inspiratory-related burst generation, but not the shape (i.e., pattern) of inspiratory motor output. However, changes in the timing of inspiratory bursts measured on XII rootlets could also result if synchronized rhythmic activity in the preBötC fails to propagate to downstream motor nuclei.

To determine whether PGE2 modulates rhythm- or pattern-generating cell populations in the preBötC, we further analyzed the distribution of individual burst periods in each slice recording (Figure 2). Burstlet theory posits that inspiratory rhythm and pattern generation are separable processes that arise from the activity of specialized neuronal subpopulations within the preBötC (Kam et al., 2013; Cui et al., 2016). If PGE2 predominately modulates the excitability of pattern-generating neurons, we would expect to see the burst period distribution become multimodal. The additional peaks in a multimodal distribution reflect dropout events where burstlets fail to elicit network-wide bursts that propagate to the XII nucleus. However, if PGE2 predominately modulates the rhythm-generating population, we would expect to see a shift or broadening of the burst period distribution that is not multimodal. We therefore tested whether the burst period distribution in our recordings became multimodal after exposure to PGE2. The distribution at baseline and after exposure to PGE2 of all bursts, grouped by experimental condition, was unimodal (Hartigan’s Dip test: p > 0.05, Figure 2). Aggregating measurements across multiple slices into a single distribution might mask effects on the individual slice level, so we further checked for multimodality in each slice and condition independently. The burst period distributions in 20 out of 21 individual slices were unimodal. The modes of the multimodal slice recording were centered around 0.93 ± 0.08 and 1.43 ± 0.16 at 100 nM PGE2 (Hartigan’s Dip test: p < 0.05) and 0.62 ± 0.16 and 1.01 ± 0.06 at 1 μM PGE2 (Hartigan’s Dip test: p < 0.05). These results suggest that PGE2 modulates inspiratory rhythm generation (i.e., timing of burstlets) rather than pattern generation (burst generation and downstream propagation).

The biphasic dose response in inspiratory rhythm frequency might be explained by differences in receptor binding affinity among Prostaglandin E2 receptors (Table 1). We attempted to distinguish the relative contributions of the EP receptor subtypes by utilizing selective antagonists while exposing rhythmic slices to PGE2 (10 nM or 1 μM, Figure 3 and Supplementary Table 3). Following a baseline recording, we preincubated each slice with EP1-4, FP, or DP2 receptor antagonists for 10 min before bath application of PGE2 (Figure 3A). For all the antagonists that we tested, the duration of the inspiratory burst period did not change in the presence of the antagonist by itself (Supplementary Table 4). Inhibiting EP2R attenuated the change in burst period normally induced by in both 10 nM PGE2 (−132% effect size compared to uninhibited controls) and 1 μM PGE2 (−100% effect size, Figure 3B). Inhibition of EP3R also attenuated the change in burst period at both 10 nM PGE2 (−72% effect size) and 1 μM (−64% effect size, Figure 3B). In contrast, inhibiting EP1 or EP4 receptors did not modify the biphasic response to PGE2. Neither XII burst area or variability was affected by incubation with PGE2 and antagonists (Figures 3C,D). Higher PGE2 concentrations (e.g., 1 μM) might additionally stimulate prostanoid receptors like DP2 and FP, which have dissociation constants in the sub-micromolar range (Table 1). In slices pre-incubated with a selective DP2 receptor antagonist (Setipiprant), the inspiratory burst period did not change after bath application of 10 nM PGE2 but it did decrease as expected in response to 1 μM. These results suggest that binding of EP2, EP3, and DP2 receptors is necessary to lengthen the burst period in 10 nM PGE2, whereas only EP2 and EP3 receptors are necessary to shorten the burst period in 1 μM PGE2.

Finally, we attempted to measure the response to PGE2 in mice that do not express EP3 receptors. We recorded inspiratory-related motor output in acute slices taken from transgenic Ptger3^–/– knock-out mice. These Ptger3^–/– mice were bred on a C57BL/6 background strain, which was different from the CD-1 wildtype mice used for our other experiments. The inspiratory burst period in both C57BL/6 wildtype controls and Ptger3^–/– knockouts did not increase in 10 nM PGE2 as it did with CD-1 mice (Figures 3E,F). In 1 μM PGE2 the inspiratory burst period decreased in both C57BL/6 wildtype (−25% below baseline) and Ptger3^–/– knockouts (−33% below baseline Figure 3F). This suggest that the stimulatory effect of PGE2 at 1 μM is conserved between CD-1 and C57BL/6 mice, and that EP3 receptors may not be developmentally necessary for 1 μM PGE2 to modulate inspiratory rhythm. However, the depressive effect of 10 nM PGE2 was not consistent between these two mouse strains, which could mean that there are strain-related differences in the sensitivity of the preBötC to low concentrations of PGE2.

Antagonizing a single receptor with specific pharmacology is useful in identifying receptor subtypes that are necessary for the biphasic PGE2 response. However, the responses we measured during receptor blockade could arise from mixed activity among the multiple prostanoid receptors that were left uninhibited. To better discriminate the actions of EP2 and EP3 receptors, we used specific agonists to stimulate EP2 and EP3 receptors independently (Figure 4A). Application of the EP2R agonist Butaprost shortened the burst period (−25% below baseline, Figure 4B). The EP3 receptor agonist, Sulprostone, lengthened the burst period (+24% above baseline, Figure 4B). A non-selective agonist that binds EP1R, EP2R, and EP3R (11-deoxy PGE1) did not change the burst period (Figure 4B). EP2R activation also decreased burst amplitude (−42% below baseline, Figure 3C), and EP3R activation decreased burst width (−21%, Figure 4D). XII burst area and the CV of the burst period also decreased after application of EP2R agonist Butaprost (area −61% below baseline, Figures 4E,F. These results suggest that biased activation of EP2 receptors is sufficient to shorten the inspiratory burst period, while biased activation of EP3 receptors is sufficient to lengthen the inspiratory burst period.

PreBötC slice recordings of hypoglossal motor output cannot readily distinguish which cell populations contained in the slice mediate the acute effects of PGE2 on inspiratory burst period. We therefore measured the relative abundances of the gene transcripts corresponding to prostanoid receptors (Ptger1-4, Ptgdr2, and Ptgfr) in single-cell RNA-Seq libraries of cells (n = 2119) that were isolated from punchouts of brainstem slices containing the preBötC (n = 10 pooled samples from littermates, age postnatal day 4, single 9600-microwell PCR plate). On average, we detected 11,247 unique molecules and 4,898 genes per cell. We detected transcripts for Ptger4 in 418 cells, Ptger3 in 292 cells, Ptger2 in 73 cells, Ptger1 in 26 cells, and Ptgfr in 48 cells (Figure 5A). We did not find any cells expressing transcripts for DP2 receptors (Ptgdr2). Clustering of the data partitioned cells into 13 groups (Figure 5B). Using cell-type-specific markers, we identified individual clusters that differentially expressed gene transcripts specific to astrocytes, microglia, oligodendrocytes, oligodendrocyte precursor cells, and endothelial cells (McKenzie et al., 2018; Supplementary Figure 1). Each of the remaining 8 clusters differentially expressed Stmn2, a neuron-specific marker, when compared to cells from non-neuronal clusters (Supplementary Table 1). We found two clusters that differentially expressed EP-receptor transcripts. Neuronal cluster 1 differentially expressed Ptger3 and Ptger4 (Ptger3: p = 0.018; Ptger4: p = 5.15e-6; Figure 5C). Cluster 1 was also the only cluster to differentially express GABAergic (Gad2, p = 5.3e-123; Slc32a1, p = 2.46e-128; Figure 5D and Supplementary Table 2) and glycinergic markers (Slc6a5, p = 1.99e-95; Figure 4D and Supplementary Table 2). Neuronal cluster 5 differentially expressed Ptger2 (p = 0.007; Figure 4C), and was also the only cluster with overlapping differential expression for Slc17a6 (VGLUT2), Tacr1 (NK1R), Sstr2 (SST receptor 2), and Cdh9 (Cadherin-9) (Figure 5D and Supplementary Table 2). These data suggest that, within the anatomical vicinity of the preBötC, cells expressing EP3 and EP4 receptors are generally more abundant than cells expressing EP1 and EP2 receptors. EP3 and EP4 receptors may be most enriched among inhibitory neurons, but our data cannot distinguish what proportion of EP3R⁺/EP4R⁺ inhibitory neurons are constituents of the preBötC rhythm-generating network. Finally, despite being detected in only 3.5% of all cells, transcripts for EP2 receptors were enriched among excitatory neurons that also differentially express markers known to overlap with the core of the preBötC (Gray et al., 2010; Gray, 2013; Yackle et al., 2017).