Genetic variation regulates opioid-induced respiratory depression in mice

In the U.S., opioid prescription for treatment of pain nearly quadrupled from 1999 to 2014. The diversion and misuse of prescription opioids along with increased use of drugs like heroin and fentanyl, has led to an epidemic in addiction and overdose deaths. The most common cause of opioid overdose and death is opioid-induced respiratory depression (OIRD), a life-threatening depression in respiratory rate thought to be caused by stimulation of opioid receptors in the inspiratory-generating regions of the brain. Studies in mice have revealed that variation in opiate lethality is associated with strain differences, suggesting that sensitivity to OIRD is genetically determined. We first tested the hypothesis that genetic variation in inbred strains of mice influences the innate variability in opioid-induced responses in respiratory depression, recovery time and survival time. Using the founders of the advanced, high-diversity mouse population, the Diversity Outbred (DO), we found substantial sex and genetic effects on respiratory sensitivity and opiate lethality. We used DO mice treated with morphine to map quantitative trait loci for respiratory depression, recovery time and survival time. Trait mapping and integrative functional genomic analysis in GeneWeaver has allowed us to implicate Galnt11, an N-acetylgalactosaminyltransferase, as a gene that regulates OIRD.


Results
Development of a high-throughput method for measuring respiratory depression in response to opioids. The PiezoSleep system measures pressure on the mouse cage floor to determine sleep or wake behaviors over small time interval. In addition to responding to pressure from gross body movements, the sensors also respond to thorax motion. The current commercial system provides estimates of respiration rates during detected sleep states; however, for this study a custom algorithm was developed to estimate respiration rates on low activity segment during wake. Figure 1 illustrates the piezoelectric sensor capability by comparing simultaneous recordings from the plesthysmography and piezoelectric sensors. Very good signal agreement is observed between the plethysmography and piezoelectric signals during periods of low activity. Figure 1A compares signals when the mouse is likely resting (no body motion). The only motion from the mouse is from the thorax due to breathing, creating pressure on the cage floor that follows the changing air flow/pressure, as tracked by the plethysmography signal. Figure 1B compare signals when there is some form of low to moderate motion, such as slow locomotion or head movements. Amplitudes represent pressure but are normalized in the figure for comparing their dynamics. Amplitude changes for the piezoelectric signals depend on posture and contact with the cage floor, while the plethysmography signals directly measure the changes air pressure from breathing. The other body movements result in greater amplitude variation in the piezo signals; however, the breathing patterns from the plethysmography signals are still observable riding on top of these larger amplitude swings (Fig. 1B). Very high activity levels result in obscuring the breathing signal, eventually to a point that it cannot be estimated reliably.
The algorithm developed for this work extracts breath rates with short-time autocorrelation functions (4 s) to suppress aperiodic signal components from body motions. 27 This approach is particularly effective on signals like those shown in Fig. 1B. Estimates are only made when signal energy is below a threshold, where body motion does not obscure the thorax pressure changes. To assess the limitations of extracting breath rates this way, the algorithm was applied to 19 simultaneous recordings (about 19 total hours of data) on mice with the plethysmography and piezoelectric sensors. Results from comparing estimates from 4 s intervals over the whole data set, show a 3% underestimation for the piezoelectric method. This combines low activity (resting/sleeping) and moderate activity intervals. For just the low activity intervals (a little less than half the data) there was a slight over estimation of 0.2% (effectively zero, which is consistent with what is expected from the agreement shown Fig. 1). For the moderate activity intervals there was an underestimation of 6%. For this study it is likely that the baseline estimate for computing respiratory depression was similar in composition to the total data set, and a 3% negative bias is expected. In the high activity period after the injections, the bias value is expected to approach the 6% limit. Therefore, it is expected that a 3% increase in the respiratory depression estimate can result from the higher number of activity intervals during the opioid response.
For the application of the algorithm to this study, estimates were performed in the low to moderate activity regions and these values were averaged over much larger intervals to provide an average respiration rate every 12 min. This was used to obtain a baseline rate (averaged over a full 24 h pretreatment) and a stable estimate of percent respiratory depression post-treatment every 12 min. This interval was large enough so that even in high activity regions there were typically 30 or more intervals where estimates could be made and averaged together. The resulting derived measures enabled a quantitative evaluation of respiratory depression, survival time and recovery time (Fig. 2).
Respiratory depression, recovery time and survival time are largely uncorrelated and heritable traits. A bracketing approach was used to construct a dose-response curve with a minimal number of mice. Eight inbred strains of mice (CC and DO founders) of both sexes were tested with an initial probe dose of 436 mg/kg of morphine. We found a wide range of outcomes for respiratory depression (49% to 77%),  PiezoSleep output for a mouse that recovers and a mouse that fails to recover after opioid treatment. Baseline respiratory rate is first established for 24 h. Time 0, is the time at which morphine is administered. The blue line represents the derived respiratory rate trajectories relative to the baseline [defined as the average respiratory rate over the first 24 h (set at 100%)]. Respiratory depression is defined as the lowest percentage of baseline reached after morphine treatment (purple dotted arrow). For a mouse that recovers, the green vertical arrow indicates recovery time when the respiratory rate returns to baseline. For a mouse that does not recover, the red vertical arrow indicates survival time when breathing stops (i.e., piezo output becomes machine noise never returning to baseline).    Based upon this literature, we selected a probe dose of 436 mg/kg. To generate the LD 50 curve using the minimum number of mice, additional higher or lower doses of morphine were included based on the response to the probe dose until an LD 50 could be estimated accurately, whereby doses on each side of the 50% mark were tested. Using this approach, survival curves by strain and sex were established (Fig. 4A) and the LD 50 for each was determined (Fig. 4B). The morphine LD 50 ranged from 212.2 mg/kg in A/J females to 882.2 mg/kg in CAST/ EiJ females. There was not a consistent sex bias up or down across all strains; instead, some LD 50 values were higher in females than males (e.g., CAST/EiJ and NOD/ShiLtJ), but higher in males than females (e.g., PWK/ PhJ and WSB/EiJ).

Genetic mapping in Do mice.
To find genetic loci that influence OIRD, quantitative trait loci (QTL) mapping was performed on the respiratory phenotypes using the high diversity, high precision, DO mouse population. A probe dose of 486 mg/kg, chosen based upon the average LD 50 of the eight founder strains and two sexes, was given to 300 DO mice, 150 of each sex. Of the 300 DO mice entered into the study, 193 (83 females, 110 males) recovered and 107 (67 females, 40 males) did not recover. The quantitative metrics for respiratory Respiratory response QtL. We next mapped QTL in the DO for the three respiratory traits using R/qtl2 30 which implements an eight-state additive haplotype model. A significant QTL was identified for respiratory depression, but no genome-wide significant QTLs were detected for recovery time or survival time using these sample sizes. One suggestive QTL was identified for overall morphine survival using a COXPH model as has previously been used for QTL mapping 32 (Fig. S1). For the respiratory depression trait, we identified a LOD 9.2 QTL with a 1.5 LOD drop interval on Chr 5:24.30-26.25 (Fig. 6A). The 95% LOD score threshold was 7.65 for p < 0.05. The QTL is called Rdro1 (respiratory depression, response to opioids 1). The Rdro1 QTL is driven by strong NOD vs. WSB/EiJ allele effects (Fig. 6B). A two-state SNP association analysis identifies potential causative SNPs and results in a reduction of the 1.5 LOD drop credible interval from the peak marker (24.98-26.16 MBp), purple SNPs in Fig. 6C. An ANOVA shows that the peak marker for the respiratory depression phenotype accounts for 42.8% of heritable variation in the DO population.
Identification of candidate genes. Genetic mapping studies are used to identify regions of interest containing variants that influence complex traits. To identify the relevant genes involved in complex trait regulatory mechanisms, there must be evidence of genetic polymorphisms segregating in the population that either influence protein structure or gene expression and evidence of a biological mechanism of action connecting them to the trait, such as expression in a trait-relevant tissue. We identified a 24.98-26.16 MBp credible interval containing 10,782 SNPs (including insertions and deletions) across the DO mice. More specifically using the SNP association mapping model, we found that 1,885 of these SNPs alleles were in the interval (Fig. 6B, Table S1). Of these 1,885 SNPs, eight were in coding regions of genes, 11 were in 5′UTRs, 16 were in 3′UTRs, 861 were intronic, 932 were intergenic, and 107 were in non-coding transcripts (some SNPs have multiple functions, Table S1). Non-coding variants typically influence phenotype through effects on regulation of gene expression via SNPs in regulatory regions, often located in the 5′ UTR, 3′ UTR, introns or within intergenic regulator features. These potential regulatory SNPs are all located in Speer4a, Actr3b, Xrcc2, Kmt2c, Galnt5, and Prkag2. While these SNPs remain candidates for regulation of the respiratory depression phenotype, we focused on coding SNPs because their impact is more readily predictable. The eight coding SNPs lie in four genes (Galtn11, Kmt2c, Speer4a, and Galnt5). None of the coding SNPs are the type with the most deleterious effects, such as a stop loss, stop gain www.nature.com/scientificreports/ or coding region insertion (frameshift). All four of these protein-coding genes are expressed in the mouse pre-Bötzinger complex, a group of brainstem interneurons that control regulation of respiratory rythmogenicity 33 , and are thus mechanistically plausible. Three of these four pre-Bötzinger complex-expressed protein-coding genes contain polymorphisms that cause non-synonymous (Cn) amino acid changes. Two of the three changes (serine to threonine at amino acid 702 in Kmt2c and asparagine to aspartic acid at amino acid 29 in Speer4a) occurred in residues that are not conserved across species. The remaining gene expressed in the pre-Bötzinger Galnt11 SNP analysis. The Galnt11 SNP, rs37913166 C/T, changes polar serine residue 545 to a hydrophobic leucine residue (Fig. 7). This amino acid is located within a functional domain that is conserved in vertebrates (Fig. 8A), suggesting that this base change regulates a key functionality. Indeed, this change places the hydrophobic residue, which are generally buried internally, onto the surface of the protein. The 3D protein structure analysis ( Fig. 8B) indicates that the non-synonymous serine-to-leucine SNP occurs in a conserved residue within the functional Ricin B lectin domain, a domain generally involved in sugar binding required for glycosylation.
integrative functional genomic analysis. GALTN11 is a member of a family of an N-acetylgalactosaminyltransferases, which function in O-linked glycosylation. The non-synonymous serine-to-leucine SNP changes a conserved residue within the Ricin B lectin domain of GALNT11 and would be predicted to impede O-linked glycosylation. GALTN11 has been shown to uniquely glycosylate at least 313 glycoproteins in HEK cells 34 . Using the Jaccard similarity tool in the GeneWeaver 35 system for integrative functional genomics, we determined that of the 313 human glycopeptides, 287 mouse orthologs are expressed in the mouse pre-Bötzinger complex 33 (Table S2). Four of those genes Hs6st2 36 , Fn1 37 , Lrp1 36 , and Sdc4 38 were identified by the Comparative Toxicogenomic Database 39 as morphine-associated genes, and one additional pre-Bötzinger-expressed glycoprotein, Cacna2d1, has been shown to have differential brain expression in mice following morphine treatment 40 . Therefore, variants in Galnt11 have a plausible mechanism of action through regulation of morphine-associated peptides in the pre-Bötzinger complex.

Discussion
In this study, we found heritable strain differences in the quantitative metrics of respiratory response to morphine, including respiratory depression, recovery time and survival time, using an advanced, high-throughput, behavioral phenotyping protocol. We further identified genomic loci involved in morphine-induced respiratory depression using an unbiased genetic approach. Mapping these traits in the DO mice and evaluation of sequence www.nature.com/scientificreports/ variants and protein structure, followed by integrative functional genomic analysis in GeneWeaver, has allowed us to implicate Galnt11 as a candidate gene for respiratory depression in response to morphine. We identified specific inbred strains of mice that were more sensitive to morphine than other inbred strains of mice. The effect of sex was bidirectional in the inbred CC/DO founder strain population where some males and females differed from each other in either direction depending on strain. The traits of respiratory depression, recovery time and survival time were all shown to have a high degree of heritability. In determining our probe dose for the outbred population, we observed that the LD 50 for morphine differed by four-fold between these eight parental strains harboring 45 million SNPs, or an equivalent genetic variation as found in the human population.
The traits of respiratory depression, recovery time and survival time were largely independent traits, as seen by their lack of correlation. Strains such as AJ, which had the lowest LD 50 for both males and females, did not demonstrate the highest degree of respiratory depression, suggesting that factors other than respiratory depression www.nature.com/scientificreports/ may play a role in opioid overdose. The lack of correlation between percent respiratory depression and survival time suggest other mechanisms of death in conjunction with OIRD. OIRD results in hypoxia and contributes to cardiac edema 41 a feature found in 96% of fentanyl and 94% of non-fentanyl opioid deaths 42 . Dolinak suggests that other baseline characteristics, such as obstructive sleep apnea, obesity, heart disease and lung disease, make an individual more susceptible to opioid toxicity 43 . Many of the mouse strains involved in this study display these traits naturally, for example NZO/HILt develop severe polygenic obesity 44 , C57BL/6 mice show spontaneous apnea and post-sigh apnea 45 , and C57BL/6J vs A/J show differences in traits related to cardiac structure and function 46 . The presence of the alleles segregating in the DO population are encouraging for finding additional QTLs related to recovery time and survival time in larger cohorts and possibly using alternatative opioids. The mechanisms underlying sensitivity to morphine and fentanyl are known to differ in many respects and this same work should be performed for the more potent fentanyl. Morphine is a µ-opioid receptor agonist with significant activities at δ-opioid receptor and κ-opioid receptor subtypes 3,5,47,48 , producing potent analgesia. In contrast, fentanyl is a potent (80-100x) high efficacy µ-opioid recepter agonist with much less activity at other opioid receptors 3,5,47,48 . Fentanyl has similarities to morphine with respect to the recruitment of intracellular signaling mechanisms but there are also key differences. For example, fentanyl causes internalization of µ-opioid receptors whereas morphine does not and the development of tolerance to morphine is JNK dependent whereas to fentanyl it is JNK independent. OR events initiated by fentanyl cause greater activation of β-arrestin complexes than do those for morphine 5,20,[49][50][51] .
To date, there have been no human GWAS or linkage studies for OIRD. Several GWAS studies have evaluated opioid use disorder (OUD) and/or opioid dependence, and have identified a variety of loci including SNPs linked to OPRM1, KCNG2, CNIH3, LOC647946, LOC101927293, CREB1, PIK3C3, and RGMA [52][53][54][55][56][57] Additional human linkage studies have identified sex-specific traits 58 , epigenetic biomarkers 59 or copy number variations 60 for risk of OUD. Other studies have sought to separate opioid use from opioid dependence, and have thus far identified SNPs associated with SDCCAG8, SLC30A9, and BEND4 61 . Only one study has looked at human opioid overdose risk, specifically by scoring overdose status and determining the number of times that medical treatment was needed in European American populations 54 . In this study, SNPs near MCOLN1, PNPLA6 and DDX18 were identified as overdose risk alleles. Human genes have thus been mapped to opioid use, opioid dependence and opioid overdose susceptibility but human studies are not able to assess opioid-induced respiratory depression, specifically the LD 50 of an opioid.
Animal studies have allowed us the opportunity to assess the LD 50 of a drug in a variety of genetic backgrounds and then map those sources of variation. These types of controlled exposure experiments cannot be conducted in humans for which exquisite control of environment is not feasible and prior exposure history is unknown. Our genetic approach of QTL mapping in the DO mouse population has allowed us to identify a genomic region containing no genes previously known to function in opioid pharmacodynamics or pharmacokinetic processes, or implicated in OUD. The genetically diverse structure of this population allows for the identification of narrow genomic intervals often with very few candidate genes. This approach of using advanced mouse populations together with integrative functional genomics has been useful for the prioritization of candidate genes in a variety of different disciplines [62][63][64] .
The identification of Galnt11 as functioning within the morphine respiratory response reveals a potential new target for therapeutic development. GALNT11 is an N-acetylgalactosaminyltransferase that initiates O-linked glycosylation whereby an N-acetyl-d-galactosamine residue is transferred to a serine or threonine residue on the target protein. The lectin domain of GALNT11 is the portion that functions to recognize partially glycosylated substrates and direct the glycosylation at nearby sites. This type of post-translational modification controls many phamacokinetic and pharmacodynamic processes as well as the regulation of delta opioid receptor (OPRD1) membrane insertions as O-linked glycosylation is required for proper export of OPRD1 from the ER 65 . O-linked glycosylation is also required for opioid binding peptides, increasing their ability to cross the blood brain barrier 66 . The integrative functional analysis in GeneWeaver identified Hs6st2 36 , Fn1 37 , Lrp1 36 , and Sdc4 38 as glycosylation targets of Galnt11. All are expressed in the pre-Bötzinger complex and are known to be responsive to morphine further supporting the concept that Galnt11 is involved in morphine-related responses. Both Hs6st2 and Lrp1 were identified as increasing and decreasing, respectively, in responsive to both morphine and stress in C57BL/6J mice 36 . Fn1 was upregulated six hours post morphine but downregulated four days later 37 , and Sdc4, a known-µ-µopioid receptor-dependent gene was also upregulated by morphine 38 . The drug-regulated expression of these known glycosylation targets of GALNT11 in relevant tissues further supports the functional relevance of Galnt11 to OIRD.
Our findings demonstrate the initial mapping of a locus involved in OIRD in mice, for which the likely candidates do not act via the opioid receptor, thereby providing a potential new target for remedial measures. Although it is through mouse genetic variation that we identified this gene, it should be noted that this gene or its glycosylation targets need not vary in humans to be a viable target mechanism for therapeutic discovery and development. Characterization of the role of Galnt11 and its variants along with other viable candidates will resolve the mechanism further, and continued mapping studies in larger populations will enable detection of additional loci for various aspects of the opioid-induced respiratory response. These findings suggest that phenotypic and genetic variation in the laboratory mouse provides a useful discovery tool for identification of previously unknown biological mechanisms of OIRD.

Methods
Mice. Male (n = 6) and female (n = 6) mice from eight inbred strains including the founders that were used www.nature.com/scientificreports/ J:DO, JAX stock number 009376) from generation 28 of outcrossing were used. All mice were acquired from The Jackson Laboratory (JAX) and were housed in duplex polycarbonate cages and maintained in a climatecontrolled room under a standard 12:12 light-dark cycle (lights on at 0,600 h). Bedding was changed weekly and mice had free access to acidified water throughout the study. Mice were provided free access to food (NIH31 5K52 chow, LabDiet/PMI Nutrition, St. Louis, MO). A Nestlet and Shepherd Shack were provided in each cage for enrichment. Mice were housed in same sex groups of three to five mice per cage. All procedures and protocols were approved by JAX Animal Care and Use Committee and were conducted in compliance with the NIH Guidelines for the Care and Use of Laboratory Animals.
comparison of piezoelectric sensors and plethysmography for measuring breath rate in mice. A mouse plethysmography chamber was built consisting of a 160 ml plexiglass chamber with ports for air supply and pressure measurement, and end openings outfitted with rubber gaskets to create an airtight seal after closing. A piezo film sensor was sealed into the bottom of the chamber for data collection with PiezoSleep software, with piezo signal sampling set at 120 Hz. Air supply to the chamber occurs through a bias flow regulator supplying air at a constant flow of 300 ml/min, drawn through the chamber via a vacuum. A differential pressure transducer (Biopac) measures the pressure difference between the plethysmography chamber and reference chamber sampled at 200 Hz.
Because the data were recorded on two different systems, an alignment program "BreathCompare", was written to automatically compensate for the difference between subjects internal clocks and to compare breath rates. Signal alignment was done through a simple correlation, which typically indicates time differences of several seconds. The program graphically displays an overlay of the piezo and plethysmography signals for easy visual confirmation. A graphical breath rate overlay allows navigation to intervals of signal disagreement to inspect the signals at these areas.
Morphine. Morphine sulfate pentahydrate (NIDA Drug Supply) was prepared at varying concentrations in sterile saline to deliver doses (200-1,200 mg/kg s.c.) in a manner not to exceed 0.2 ml/10 g body weight. Starting with a dose of 436 mg/kg for all strains and depending upon the result of that dose, the next dose was either increased (536 mg/kg) or decreased (336 mg/kg) for the next cohort, such that doses flanking both sides of the 50% survival point (LD 50 ) were tested. This was repeated with increasing and decreasing doses of 100 mg/kg as necessary depending upon the results of the previous dose. If smaller increments were needed, based upon the results of flanking doses, the dosage was altered by 50 mg/kg or 25 mg/kg. Not all strains received all doses but each strain received at least three doses such that two flanked (one above, one below) the LD 50 . Using this approach and testing 3-4 doses per strain, the LD 50 for the eight strains and both sexes was determined and ranged from 212 to 882 mg/kg. piezoelectric sleep monitoring to determine respiratory depression, recovery time and survival time. Mice were placed individually into the 7 × 7 inch piezoelectric grid and chamber system (Signal Solutions) for 24 h to equilibrate to the apparatus and collect baseline activity data 27,67 . The mice were housed with approximately ½ cup of pine bedding and ALPHA-dri (Shepard). Individual testing is necessary due to the known enhanced lethality of cage mates during morphine exposure, which has been shown to affect survival 68 . The mice had access to food and water ad libitum while in the chamber. The room was maintained on a 12:12-h light:dark cycle. To control for known circadian effects 69 mice were placed in the chambers between 9 am − 12 pm on Day 1 and were injected with morphine 24 h later. They remained in their chambers undisturbed until 24 h after injection. Whenever possible, complete balanced cohorts of the eight strains and both sexes were run during each of nine replicates of the experiment. The data acquisition computer, food and water were checked daily; otherwise, the mice remained undisturbed. Breath rates were estimated from 4-s intervals in which animal activity dropped low (i.e. during sleep and brief rest periods and pauses during wake), and averaged over 24-min overlapping intervals to provide an average respiratory rate every 12 min. The respiratory rate baseline consisted of the average respiratory rate over the first 24 h, which included both sleep and wake periods. Respiratory rate was then measured in the same way after injection of opioid. These measures were then used to determine thresholds for obtaining the recovery time (respiratory rate returns to baseline, see Fig. 1) or survival time (animal stops moving and breathing, never returning to threshold, see Fig. 1). The 300 DO mice were tested in random cohorts of 36 using the PiezoSleep system with a single 486 mg/kg dose of morphine. This dose was determined as the average LD 50 dose across the eight strains and two sexes, 16 samples. calculating respiratory depression, recovery time and survival time. To test for difference in the respiratory depression, recovery time and survival time across the strains and sexes a linear model was fit, the full model was: where phenotype was respiratory depression, recovery time or survival time and where ε is random error. The β-parameters were estimated by ordinary least squares, and the type III sum of squares was considered ε in the ANOVA model. In all cases, the full model was fit and reduced by dropping non-significant interactions followed by main effects. calculating broad sense heritability (H 2 ). As a measure of broad sense heritability in the founder strains, the ICC was determined using ICCest from the ICC 2.3.4 70  www.nature.com/scientificreports/ Morphine LD 50 data analysis. The LD 50 was calculated using the drc 3.0-1 library 71 in R using the total tested and the observed dead at each dose as a binary or binomial response. A logistic regression model was fit, and a goodness of fit test (based upon Bates and Watts 72 ) performed. In addition, a regression model assuming equal LD 50 across strains was compared by chi square to an LD 50 assumed different across strains. An adjusted and unadjusted 95% credible interval was also calculated. To graph the data a non-linear 2-parameter model was fit in JMP 14.2.0 (RRID:SCR_014242) with the formula: where a = growth rate and b = inflection point.
Genotyping. Tail samples were collected at the conclusion of the experiment and all mice were genotyped using the Giga-MUGA genotyping array (NeoGene). Data were deposited in the DO Database (DODB) RRID:SCR_018180). Genotypes were imputed to a 69 K grid to allow for equal representation across the genome.
QtL mapping. Genotype probabilities were calculated according to the founder genotypes and then converted to allele probabilities. We then interpolated allele probabilities into a grid of 69,000 evenly-spaced genetic intervals 73 . We performed genome scans using R/qtl2 RRID:SCR_018181 31 . Sex and date of test were included as additive covariates. The model includes the random effect of kinship among the DO animals computed using the LOCO method 74 . The significance thresholds were determined for each trait by permutation mapping 75 . The confidence interval around the peak makers was determined using Bayesian support intervals. To determine the percent of variation accounted for by the QTL the mice were classified at the variant into one of eight states based on genotype probabilities of mice at that locus. Following this, a one-way ANOVA was fit to ascertain the strain variation relative to total variation towards estimating heritable variation at that locus. SNP association mapping was also performed using R/qtl2 to test the association of individual SNPs alleles in the region of the locus with the phenotype. Briefly, SNP data were obtained from the SANGER and MGI databases 76,77 for the interval. Using the genotype probabilities and the founder SNP genotypes to infer the SNP genotypes of the DO mice. At each SNP location the eight allele state probabilities are collapsed to two state SNP probabilities and the Cox proportional hazards regression was performed by coxph function in the survival (3.1-12) R package. Based on the output of the log (base e) likelihood for the null model and for the alternative model (with covariates and genotype probabilities), we took the difference of both log likelihoods and then divided by ln (10) to convert the results into the LOD scale. Full QTL mapping scripts are available https ://theja ckson labor atory .githu b.io/DO_Opioi d/index .html.
candidate gene analysis. In order to assess the plausibility of genes in the QTL interval we identified all SNPs within the additive SNP model segregating between the high and low alleles of the DO founder strains. Next we identified those that were within protein coding region that were most deleterious. Differential coding sequence non-synonymous amino acid substitution SNPs (Cn) that differed between the high and low allele groups were identified. GeneWeaver's database (RRID:SCR_003009) was searched to identify the overlap among tissue-specific expression profiles from Allen Brain Atlas as well as datasets derived from Entrez GEO profiles (RRID:SCR_004584) for pre-Bötzinger neurons 33 (GS 273275), diaphragm (GS273269) 78 and lung 79 with the QTL positional candidates. The gene sets were overlapped using the Jaccard similarity and GeneSet graph tools 80 .
In order to determine if the Cn SNPS were in areas of evolutionary conservation we aligned the sequence of several species. Representative sequences for each species Drosophila (Q8MVS5), Xenopus (Q6DJR8), Danio (Q08CC3), Rattus (Q6P6V1), Mus (Q921L8), Homo (Q8NCW6) were acquired from Uniprot (RRID:SCR_002380) and aligned. The Clustal Omega program was used with default parameters 81 . The transition matrix is Gonnet, gap opening penalty of six bits, gap extension of one bit. Clustal-Omega uses the HHalign algorithm and its default settings as its core alignment engine 82 . To determine where the three-dimensional effects of the amino acid change would be we obtained the 3D crystal structure (1XHB) 83 from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RRID:SCR_012820) and visualized it with Jmol (RRID:SCR_003796) 84 . integrative functional genomics. In order to assess the functional sufficiency of Galnt11 as a candidate gene the literature was searched to identify genome-wide studies characterizing glycosylation targets of GALTN11, one study was identified 33 and these genes were added to the GeneWeaver Database (GS356053). Using the Jaccard similarity tool, we overlapped the glycosylation targets with the genes expressed in the mouse pre-Bötzinger complex. We next overlapped these gene with genes identified by the Comparative Toxicogenomic Database (RRID:SCR_006530) as morphine-associated genes.

Data availability
The datasets generated during and/or analysed during the current study are available in the Mouse Phenome Database ProjectID:Bubier3; and in the DO Database (DODB).

code availability
The R markdown document for the QTL mapping is available at https ://theja ckson labor atory .githu b.io/DO_ Opioi d/Plot_DO_morph ine.html.