Postoperative nausea and vomiting (PONV) is a key driver of unplanned admission and patient satisfaction after surgery. Because traditional risk factors do not completely explain variability in risk, this study hypothesized that genetics may contribute to the overall risk for this complication. The objective of this research is to perform a genome-wide association study of PONV, derive a polygenic risk score for PONV, assess associations between the risk score and PONV in a validation cohort, and compare any genetic contributions to known clinical risks for PONV.
Surgeries with integrated genetic and perioperative data performed under general anesthesia at Michigan Medicine (Ann Arbor, Michigan) and Vanderbilt University Medical Center (Nashville, Tennessee) were studied. PONV was defined as nausea or emesis occurring and documented in the postanesthesia care unit. In the discovery phase, genome-wide association studies were performed on each genetic cohort, and the results were meta-analyzed. Next, the polygenic phase assessed whether a polygenic score, derived from genome-wide association study in a derivation cohort from Vanderbilt University Medical Center, improved prediction within a validation cohort from Michigan Medicine, as quantified by discrimination (c-statistic) and net reclassification index.
Of 64,523 total patients, 5,703 developed PONV (8.8%). The study identified 46 genetic variants exceeding the threshold of P < 1 × 10−5, occurring with minor allele frequency greater than 1%, and demonstrating concordant effects in both cohorts. Standardized polygenic score was associated with PONV in a basic model, controlling for age and sex (adjusted odds ratio, 1.027 per SD increase in overall genetic risk; 95% CI, 1.001 to 1.053; P = 0.044), a model based on known clinical risks (adjusted odds ratio, 1.029; 95% CI, 1.003 to 1.055; P = 0.030), and a full clinical regression, controlling for 21 demographic, surgical, and anesthetic factors, (adjusted odds ratio, 1.029; 95% CI, 1.002 to 1.056; P = 0.033). The addition of polygenic score improved overall discrimination in models based on known clinical risk factors (c-statistic, 0.616 compared to 0.613; P = 0.028) and improved net reclassification of 4.6% of cases.
Standardized polygenic risk was associated with PONV in all three of the study’s models, but the genetic influence was smaller than exerted by clinical risk factors. Specifically, a patient with a polygenic risk score greater than 1 SD above the mean has 2 to 3% greater odds of developing PONV when compared to the baseline population, which is at least an order of magnitude smaller than the increase associated with having prior PONV or motion sickness (55%), having a history of migraines (17%), or being female (83%) and is not clinically significant. Furthermore, the use of a polygenic risk score does not meaningfully improve discrimination compared to clinical risk factors and is not clinically useful.
Postoperative nausea and vomiting is a common complication after surgery, which decreases patient satisfaction and can lead to unplanned admission
Known risk factors only partially explain interindividual variation in risk, suggesting that genetics may play a role
The investigators identified genetic variants associated with postoperative nausea and vomiting (PONV) by performing a genome-wide association study
Based on this study, they created a polygenic risk score for PONV in a derivation cohort; in a validation cohort, when added to traditional risk factors, the polygenic risk score only modestly augmented the prediction of PONV
The use of this polygenic risk score did not result in a clinically meaningful improvement in PONV prediction when added to traditional risk factors
Postoperative nausea and vomiting (PONV) is associated with as many as 30% of surgical procedures and, despite improvements in patient screening and prophylaxis, remain the most frequently encountered side effects of general anesthesia.1 PONV affects the surgical experience through discharge delays, increased resource utilization, and decreased patient satisfaction.2 As the number of surgical patients with genetic information continues to increase,3 an improved understanding of the genetic mechanisms of PONV4 could enable genotype to be incorporated into risk stratification and treatment.
Genome-wide association studies for perioperative phenotypes have mostly identified genetic loci with small effects leading to difficulty integrating this information with patient and operative factors, which frequently exert stronger influence on overall risk.5,6 Polygenic risk scores aggregate associations (derived from genome-wide association studies) across millions of genetic variants to provide a baseline estimate of an individual’s genetic burden for a specific disease or trait.7 Although the impact of any individual single genetic variant on perioperative complications is often lower than traditional clinical risk factors and may not have clinical impact by itself,5,8 the additive effect of the genetic burden, as quantified by polygenic risk, can have a meaningful impact.9,10 Furthermore, polygenic risk score can be easily incorporated as a single variable (with normal distribution) into existing risk models, instead of using thousands of variables with weak effect on the disease risk (jeopardizing the validity of the statistical model).7 Attempts to elucidate the genetics of PONV risk have largely used the candidate gene approach, studying genes and variants already hypothesized to play a role, such as polymorphisms in serotonin receptors,11–14 muscarinic acetylcholine receptors,1 dopamine receptor,15,16 and cytochrome P450.17 The largest genome-wide association studies involving PONV have been limited to a few hundred patients18,19 and lacked statistical power to detect associations at a genome-wide significance threshold. Other studies have been limited by discrepancies in how PONV outcomes are defined20,21 and negative or conflicting results.14,22 Large biorepositories greatly expand the scale of genetic studies possible.23,24 However, most of these biorepositories lack integration with perioperative data, which is necessary to study perioperative phenotypes such as PONV.
We hypothesize that genetics may contribute to the overall risk for PONV. The objective of this research is to (1) perform a genome-wide association study of PONV, (2) derive a polygenic risk score for PONV, (3) assess associations between the risk score and PONV, and (4) compare any genetic contributions to known clinical risks for PONV. Accomplishing this objective first required a discovery phase, during which genome-wide association studies were performed on each genetic cohort and the results were meta-analyzed. Next, in the polygenic phase, we assessed whether a polygenic score derived from genome-wide association study in a derivation cohort from Vanderbilt University Medical Center (Nashville, Tennessee) improved prediction within a validation cohort from Michigan Medicine (Ann Arbor, Michigan), as quantified by discrimination (c-statistic) and net reclassification index. Finally in the polygenic phase, we characterized the significance and contribution of polygenic score to overall PONV risk relative to both traditional25 and advanced (i.e., nongenetic risk factors, extending beyond widely utilized, published models) clinical risk factors within the validation cohort.5 This study was enabled by standardization of perioperative phenotypes and integration of data from two large academic biorepositories.
Materials and Methods
All study procedures were approved by the Institutional Review Boards at Michigan Medicine (approval No. HUM00127909) and Vanderbilt University Medical Center (approval No. 220622). The STrengthening the REporting of Genetic Association Studies (STREGA)26 and the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD) reporting guidelines were followed.27 Informed patient consent was obtained at the time of enrollment. Before the data were accessed, the study methods including data collection, outcomes, and statistical analyses were written and filed at institutional peer-review committee meetings (Anesthesia Clinical Research Committee28 at Michigan Medicine, January 12, 2022; and Vanderbilt Institute for Clinical and Translational Research29 at Vanderbilt University Medical Center, July 8, 2022).
Patient Populations
The discovery cohort was created from patients with available perioperative and genetic data from enrollment in one of two institutional biobanks (39,261 from Michigan Medicine and 25,262 from Vanderbilt University Medical Center). Our study utilized perioperative data (local repository) from two participating Multicenter Perioperative Outcomes Group sites (Michigan Medicine and Vanderbilt University Medical Center).30
Patients aged 18 yr and older who underwent a surgical procedure lasting greater than 45 min and requiring general anesthesia between 2014 and 2020 were included. Additionally, selected cases met a previously validated criteria for data completeness and quality (appendix 1) and passed site-specific genotype quality control (appendix 2). Patients transported directly to the intensive care unit, those receiving labor epidurals, and those undergoing deceased donor organ procurement or liver transplantation were excluded. A flow diagram depicting cohort generation at each institution and how many patients were excluded at each step can be found in supplemental figure S1 (https://links.lww.com/ALN/D687).
Choice of preoperative/intraoperative antiemetics was left to provider discretion at each institution. Michigan Medicine and Vanderbilt University have institutional practice recommendations for prophylaxis of PONV, which are based upon the Fourth Consensus Guideline for the Management of Postoperative Nausea and Vomiting.36 In brief, the guidelines recommend that patients with one or two risk factors (female sex, younger age, nonsmoker, surgery type, history of PONV or motion sickness, and opioid analgesia) receive prophylaxis with two agents, and those with three or more risk factors receive three or four agents. In practice, at Michigan Medicine, standard-risk patients without contraindication typically receive dexamethasone (4 mg) IV at induction and diphenhydramine (12.5 to 25 mg) IV at induction. Higher-risk patients receive aprepitant (40 mg) orally preoperatively. Ondansetron (4 mg) IV is typically reserved for administration postoperatively. At Vanderbilt University Medical Center, ondansetron is administered as first-line therapy (barring contraindications). Higher-risk patients receive preoperative scopolamine patch. Dexamethasone and haloperidol (0.5 to 1 mg) are commonly used second- and third-line agents.
Data Collection
Clinical Variables
Variables were derived from the local Multicenter Perioperative Outcomes Group data set. The variables were divided into (1) patient factors (age, sex, body mass index, comorbidities, and others); (2) known clinical risk factors (female sex, history of motion sickness or PONV, nonsmoking status, and postoperative opioid administration)25 ; (3) intraoperative details (emergency surgery, surgery classification, case duration, and others); and (4) anesthetic factors (neuraxial anesthesia, regional block, nitrous oxide, intraoperative opioids, and preoperative/intraoperative antiemetics; appendix 3). A full list of variables can be found in supplemental table S1 (https://links.lww.com/ALN/D687). Note that sex (extracted from the electronic health record and not genotypic sex) was used in the place of “female gender” from the original study by Apfel et al.25
Genotype
Genetic data were obtained from Michigan Medicine’s Michigan Genomics Initiative and Vanderbilt’s DNA biobank (BioVU). Study participants provided written informed consent for obtaining genetic data, and protocols were reviewed and approved by the institutional review boards.5 Due to the retrospective design of this study and the use of deidentified data, the need for additional informed patient consent was waived by both boards. Full details and characterization of genotype for both cohorts can be found in appendix 4.
Primary Outcome: PONV in the Postanesthesia Care Unit
The primary outcome was nausea or emesis occurring and documented in the postanesthesia care unit (PACU). Because we could not say with certainty whether antiemetics administered in the PACU were given as “rescue” treatment for PONV or “prophylaxis” as part of routine postoperative care, patients receiving antiemetics (as detailed in appendix 3) in the PACU (but no documented nausea or vomiting) were excluded from subsequent analysis (i.e., not included in either the “PONV” or “No PONV” groups).
Secondary Outcome: PONV or Antiemetics in the PACU
We repeated the full analysis (both genome-wide association study and polygenic scoring) using a composite of (1) nausea or emesis occurring and documented in the PACU and (2) antiemetics (appendix 3) administered in the PACU. The secondary outcome was selected to match the phenotype developed by the quality improvement arm of the Multicenter Perioperative Outcomes Group, known as the Anesthesiology Performance Improvement and Reporting Exchange.38,39 The intersection and union of each outcome (PONV in the PACU and antiemetics in the PACU) and how each contributes to the primary and secondary outcomes can be seen in figure 1B.
Discovery Phase
The aim of the discovery phase is to identify genetic variants (and associated genes) that are associated with PONV. A combined cohort, consisting of both Michigan Medicine and Vanderbilt University Medical Center patients (meta-analyzed) provides higher statistical power for novel genetic discovery (compared to the genome-wide association study repeated in the more limited, derivation cohort of the polygenic phase). The methods employed in the discovery phase include (1) genome-wide association study and (2) gene prioritization.
Outcome in the Discovery Phase
A genome-wide association study was conducted at the patient level (a patient’s genetic risk is static; fig. 1); however, a patient may have multiple qualifying surgical procedures within our inclusion window and may develop PONV after some procedures but not others. For the discovery phase (i.e., genome-wide association study and gene prioritization), patients had PONV if they met the criteria for primary outcome (or secondary outcome) on all procedures, and patients were considered “No PONV” if they met criteria on none of their surgical procedures. Patients with qualifying outcome on some, but not all, procedures were considered equivocal and omitted from analysis in the discovery phase of the study.
Genome-wide Association Study
Genome-wide association studies were performed on each genetic cohort independently using a logistic mixed model with saddle point estimation implemented to account for sample relatedness and case-control imbalance (appendix 5). Quantile–quantile plots were generated to assess inflation. A genome-wide significance threshold (P < 5 × 10−8) was selected a priori for identification of novel loci,40 and a lower threshold (P < 1 × 10−5)41,42 was selected as a suggestive threshold (and for replication of previously reported loci). Details on the harmonization of summary statistics and meta-analysis can be found in appendix 6.
Power Analysis
We used a convenience sample of qualifying cases. No a priori power analysis was performed. There were not enough patients of non-European ancestry (enrolled in Michigan Genomics Initiative or BioVU) to perform a more comprehensive analysis in cohorts of non-European ancestry; therefore, our study population was restricted to patients of European descent.
Gene Prioritization
Genome-wide association studies frequently identify variants located in intronic and intergenic regions, which hinders their functional evaluation (i.e., nonexpressed variants are unlikely to exert causal or mechanistic effect on PONV). Instead, these associations potentially reflect linkage disequilibrium (nonrandom association of multiple DNA markers that results from proximity within chromosome and resulting coinheritance) with other variants that were not detected. Polygenic priority scoring addresses this possibility by combining genome-wide association study results (variant-based) with a linkage disequilibrium reference panel and other gene–gene correlations matrices to rank (or “prioritize”) genes that may play a role in the pathogenesis of PONV.46 This leverages the full genetic signal and incorporates data about genes from a variety of sources to prioritize causal genes with increased precision. This method assessed whether PONV-associated variants overlap with genes or other regulatory features more often than would be expected based upon control sets prioritizing candidate genes. Polygenic priority scoring collapses summary statistics (obtained from genome-wide association studies) of variants from the same genetic region to estimate the gene-level contribution, using a linkage disequilibrium reference panel from the UK Biobank European-ancestry individuals47 and further generated gene–gene correlation matrix to perform enrichment analysis. Genes considered as prioritized candidate genes were defined by having a polygenic priority score in the top 10% across the whole genome in our study.
Polygenic Phase (Derivation and Validation)
We hypothesized that genetics may contribute to the overall risk for PONV. Testing this hypothesis required two aims in the polygenic phase: (1) calculating a polygenic score for PONV (“derivation”) and then (2) testing whether this polygenic score is associated with PONV and improves risk prediction (“validation”).
Outcome Differences in the Polygenic Phase
For the polygenic analysis, the patient’s static polygenic score for PONV was included with other procedure-specific variables (for example, surgery duration and nitrous oxide administration). If the patient underwent multiple surgical procedures within a 30-day time window, only the first procedure was considered. If a patient underwent more than five surgical procedures, additional procedures were excluded to not bias the overall analysis.
Derivation of Polygenic Score
The aim of the derivation phase is to perform a genome-wide association study in an independent data set, the results of which can subsequently be used to generate a polygenic score. The BioVU population was selected for derivation. Summary statistics and observed prevalence from genome-wide association studies in the BioVU population were combined with an external linkage disequilibrium reference panel (UK Biobank, European population)47 to generate a polygenic score for both PONV phenotypes in the validation Michigan Medicine population (fig. 1).48 Polygenic scores were calculated via the method of Bayesian regression and continuous shrinkage priors (the PRS-CS method) using default Strawderman–Berger prior (A = 1, B = 0.5).49 Polygenic scores were calculated using the effective sample size for BioVU (calculated as 4/[1/Ncases + 1/Ncontrols]) and with age, genotype inferred sex, genotyping array, and the first four principal components as covariates. The polygenic score for PONV for each patient in the Michigan Medicine population was then converted into a standardized score with mean 0 and SD of 1 to facilitate interpretation.
The derivation and validation cohorts were selected based upon the availability of additional clinical features available in the Michigan Medicine cohort that were incorporated into the full clinical model. Prediction accuracy of polygenic scoring depends mostly on the heritability of the analyzed traits, the number of single-nucleotide polymorphisms, and the size of the discovery sample.6 The full set of genetic associations with minor allele frequency greater than or equal to 0.01 from the derivation cohort were utilized in polygenic score generation.
Validation of Polygenic Score
Multivariable Regressions
Perioperative characteristics were summarized using means and standard deviations for normally distributed continuous covariates, medians and interquartile range for nonnormally distributed continuous variables, and counts and percentages for categorical covariates. Statistical analysis was performed in R version 3.5.1 (R Foundation for Statistical Computing, Austria). With a missing rate of less than 5% across study variables, complete case analyses were used and considered unbiased. The major known clinical risk factors (female sex, nonsmoker, opioid administered, and prior history of PONV or motion sickness) were collected along with the PONV outcome as part of a quality metric; therefore, there were no missing data for these covariates within the cases that PONV outcomes were present (supplemental fig. S1A, https://links.lww.com/ALN/D687). Missingness in other key variables included American Society of Anesthesiologists (Schaumburg, Illinois) Physical Status (N = 147), emergency surgery (N = 2), and surgical service (N = 0). We tested for independent association between polygenic score and PONV in a basic multivariable regression (controlling for the patient characteristics age and sex [model 1]). Next, we assessed association between polygenic score and PONV while controlling for known clinical risk factors (as described by Apfel et al.):25 (1) female sex, (2) prior history of either PONV or motion sickness, (3) nonsmoking status, and (4) opioid administration (model 2). Finally, we assessed association in a full multivariable regression (controlling for demographic details, comorbidities, intraoperative, and anesthetic factors [model 3]). We also tested each of our models using a secondary outcome. No additional sensitivity analyses were performed. The results are presented as the adjusted odds ratio, a measure of association between an exposure and an outcome, where 1 means no association.
Area under the Receiver Operating Characteristic Curves (c-statistic)
Performance characteristics for each of the three models (with and without polygenic score for PONV) were assessed using receiver operating characteristic curves. The c-statistics resulting from the receiver operating characteristic curves were compared using the Hanley and McNeil method, a statistical method for comparing two c-statistics calculated from different models with the same data.50, 51
Net Reclassification Index
We used a continuous net reclassification index to quantify change in classification by adding the polygenic score for PONV to the model. The continuous net reclassification index uses the continuous predicted values from each regression, evaluating events with higher predicted values as having improved and nonevents with lower predicted values as having improved.52,53 This allows for a category-free approach to net reclassification index, which removes some of the decisions of cutoff values and categories that can result in misleading values.
Results
Of the 64,523 patients included in our study, 5,703 (8.8%) developed PONV in the PACU (primary outcome). Of 69,706 patients considered in the secondary analysis, 10,886 (15.6%) developed PONV or were administered an antiemetic while in the PACU (secondary outcome; fig. 1A).
Discovery Phase
Cohort (Michigan Medicine and Vanderbilt University Medical Center)
The average age of the Michigan Medicine cohort was 54 yr (SD, 16 yr), and 53.3% were female. The average age of the Vanderbilt cohort was 56 yr (SD, 16 yr), and 51.5% were female (supplemental table S2, https://links.lww.com/ALN/D687).
Genome-wide Association Study
We tested the association between 42,124,633 and 37,973,293 genetic variants and PONV for genome-wide association studies 1 and 2, respectively. Of these variants, 16,564,986 and 8,004,272 appeared in both Michigan Medicine and Vanderbilt University Medical Center for genome-wide association studies 1 and 2, respectively. The difference in the number of genetic variants tested in study 1 versus study 2 can be explained by the overlapping but not identical cohorts (illustrated in fig. 1) used for each study. Because the participants in each study differed, the number of variants passing minor allele count and minor allele frequency thresholds (described in appendix 4) also was different. The Manhattan plot showing the results of the genome-wide association study meta-analysis and the quantile–quantile plots for genome-wide association studies 1 (primary outcome) and 2 (secondary outcome) are shown in figure 2. Quantile–quantile plots for single-nucleotide polymorphisms can be found in supplemental figure S2 (https://links.lww.com/ALN/D687).
Genome-wide association study 1 (primary outcome: nausea or emesis in the PACU) identified 46 variants exceeding the threshold of P < 1 × 10−5, occurring with minor allele frequency greater than 1%, and demonstrating concordant effects in both cohorts. Genome-wide association study 2 (secondary outcome: nausea/emesis or antiemetics in the PACU), identified one low-frequency (minor allele frequency = 0.002) variant (rs760998109) on chromosome 2 associated with PONV at the level of genome-wide significance (P = 4.4 × 10−8). However, this variant was only present in the Michigan Medicine data set and could not be assessed in our meta-analysis. Genome-wide association study 2 also had 33 variants exceeding the threshold of P < 1 × 10−5, occurring with minor allele frequency greater than 1%, and demonstrating concordant effects in both cohorts. The full list of selected associations (P < 1.0 × 10−5; minor allele frequency greater than 1%, and concordant effect in both independent biobanks) can be found in table 1.
Summary statistics are available for download in the genome-wide association study catalog (National Human Genome Research Institute-European Bioinformatics Institute Catalog of Human Genome-wide Association Studies). We also specifically report effect and significance for all variants previously reported to be associated with PONV (supplemental table S3, https://links.lww.com/ALN/D687).
Gene Prioritization
To highlight biology likely involved in this phenotype, we searched for candidate functional genes using polygenic priority scoring.54 Polygenic priority scoring is a method for prioritizing genes by integrating descriptive statistics from genome-wide association studies with gene expression, biologic pathways, and protein–protein interaction data. Genome-wide association studies commonly identify genetic variants in noncoding regions, which limits direct inference of biologic insights and mechanisms.55 ,56 Gene prioritization is a computational approach that attempts to identify causal genes by (1) accounting for noncoding variants in linkage disequilibrium with the causal gene and (2) mapping gene regulatory elements acting on the causal gene.54 Polygenic priority scoring was selected as the gene prioritization method for use in our study because it outperformed other commonly used gene prioritization algorithms in a meta-analysis of 46 traits.54 The genes prioritized are listed in supplemental table S4 (https://links.lww.com/ALN/D687) and include MAPT, PLA2G4A, THRB, NRP1, IL10, GJA1, EDNRB, and TGFB1 (which have been previously reported)56 ; SKP1, ACTN2, SKP1, ACTN2, COL19A1, and PLXNA4 (which have homologous genes reported)56 ; and HEY2, MITF, FCRL2, SPI1, BBS4, and PCGF2 (which are novel candidate genes).
Polygenic Phase
Univariate Analysis: Derivation
The derivation cohort is the subset of the discovery cohort from BioVU (fig. 1). Additional details on the derivation population are reported in supplemental table S2 (https://links.lww.com/ALN/D687).
Univariate Analysis: Validation
Younger patients, females, and nonsmokers were more likely to develop PONV (table 2). Additionally, heavier patients (body mass index mean, 30.9 ± 8.0 kg/m2vs. 29.8 ± 7.0 kg/m2; P < 0.001) were more likely to develop PONV based on the results of the univariate analysis. Patients with comorbidities including anxiety and migraines had a higher incidence of PONV. PONV occurred more frequently after longer procedures (213 ± 123 min vs. 183 ± 122 min; P < 0.001), procedures utilizing nitrous oxide (76.9% vs. 74.6%; P = 0.001), and procedures with higher administration of opioids (0.25 ± 0.16 oral morphine equivalents/[kg × min] vs. 0.23 ± 0.16; P < 0.001). Note that the calculated difference in normalized oral morphine equivalents (0.02 mg/[kg × min]) would equate to 588 μg of IV fentanyl, assuming a 70-kg person undergoing a 2-h procedure.
On a procedural level, 6,786 (11.0%) of 61,503 qualifying procedures at Michigan Medicine met the primary PONV outcome, and 14,425 (20.9%) of 69,142 qualifying procedures met the secondary, composite diagnostic criteria (fig. 1B). The 7,639 patients meeting criteria for the secondary outcome but not meeting criteria for the primary outcome (i.e., those receiving antiemetics in the PACU, without documentation of PONV) were excluded from analyses using the primary outcome. These 7,639 patients account for the difference in total number of qualifying procedures between the analyses using the primary and secondary outcomes (61,503 and 69,142, respectively). The intersection and union of the primary and secondary outcomes can be seen in figure 1B. Otolaryngology and oral and maxillofacial surgery (15.5%), urology (14.7%), and orthopedic surgery (13.8%) were the most common classes of procedure (table 2). Additional details on the validation population can be found in table 2.
Polygenic Score for PONV
Associations across 1,048,398 genetic variants (genome-wide association study 1) and 1,061,699 genetic variants (genome-wide association study 2) from the derivation cohort were used to calculate a polygenic score for PONV for every patient in the validation cohort. The association between polygenic score for PONV and PONV in the validation cohort was tested using multivariable regression models.5 The mean polygenic score was 0.025 ± 1.014 for patients developing PONV and −0.002 ± 0.999 for those not developing PONV (P = 0.020; table 2). Patients with a polygenic score for PONV greater than 1 (1 SD above the mean) had a 0.6% greater risk of developing PONV on univariate analysis compared to the baseline population (defined as polygenic score for PONV between −1 and 1), with absolute risks of 11.5% and 10.9%, respectively. Furthermore, the patients with the highest quintile of risk have 0.5% increased risk of developing PONV compared to those in the lowest quintile of risk (absolute risk, 11.4% and 10.9%, respectively). The polygenic score for PONV was significantly associated with the primary outcome without adjustment for any other covariates (odds ratio, 1.028; 95% CI, 1.002 to 1.054; P = 0.034; c-statistic, 0.509; 95% CI, 0.501 to 0.516). The polygenic score for PONV was also significantly associated with the secondary outcome without adjustment for other covariates (odds ratio, 1.021; 95% CI, 1.003 to 1.040; P = 0.025; c-statistic, 0.506; 95% CI, 0.501 to 0.5112).
Incorporation of Polygenic Score into Multivariable Regression Models
The polygenic score for PONV was incorporated into three regression models: (1) a basic model controlling for age and sex, (2) a model incorporating known clinical risk factors,25 and (3) a model built from the full clinical data (patient-level characteristics, known clinical risk factors, intraoperative surgical factors, and intraoperative anesthetic factors). These regressions were repeated for both the primary (nausea or emesis in the PACU) and secondary (nausea or emesis or antiemetics in the PACU) outcomes (table 3). The standardized polygenic score was significantly associated with the primary PONV metric in the basic (adjusted odds ratio, 1.027; 95% CI, 1.001 to 1.053; P = 0.044), the known clinical risk factors model (adjusted odds ratio, 1.029, 95% CI, 1.003 to 1.055; P = 0.030), and the full clinical model (adjusted odds ratio, 1.029, 95% CI, 1.002 to 1.056; P = 0.033). Furthermore, the standardized polygenic score was significantly associated with the secondary outcome in the basic model (adjusted odds ratio, 1.020; 95% CI, 1.001 to 1.039; P = 0.034) and the known clinical risk factors model (adjusted odds ratio, 1.022; 95% CI, 1.003 to 1.041; P = 0.024) but not the full clinical model (adjusted odds ratio, 1.019; 95% CI, 1.000 to 1.039; P = 0.051). The results of each model can be found in table 3.
The four variables validated in the study by Apfel et al.25 were all highly significant in our regression. In the full clinical model, a variety of patient factors (including body mass index, history of anxiety, and history of migraines), surgical (including surgical classification and case duration), and anesthetic factors (including inhalational anesthetic, neuraxial procedure, and intraoperative opioid administration) were shown to be significantly associated with PONV.
Characterizing Model Performance by c-Statistic
Incorporation of polygenic score for PONV improved model discrimination (assessed by c-statistic) of the known clinical risk factors model but not the basic or the full model. The c-statistic of the full clinical model was 0.699 (95% CI, 0.693 to 0.705) compared to the basic model, which had a discrimination of 0.624 (95% CI, 0.617 to 0.631), or known clinical risk factors model, which had a discrimination of 0.616 (95% CI, 0.609 to 0.623; table 3). Inclusion of polygenic score for PONV with the known clinical risk factors model increased discrimination from the comparable model without polygenic score for PONV for both the primary and secondary outcomes (c-statistic. 0.616 vs. 0.613; P = 0.028; and c-statistic. 0.610 vs. 0.608; P = 0.006, respectively; table 4). Inclusion of polygenic score for PONV did not improve discrimination from the comparable model without the polygenic scores for PONV for the basic or full clinical models. The full clinical model had improved discrimination compared to the basic or known clinical risk factors models. A comparison of model discrimination between analogous models (with and without inclusion of polygenic score for PONV) is reported in table 4.
Characterizing Model Performance by Net Reclassification Index
Net reclassification indices suggest that the inclusion of polygenic score for PONV improved classification in between 4 and 5% of procedures for the primary outcome and between 1 and 2% of procedures for the secondary outcome. Net reclassification index values quantifying the change in classification between each of the three models are reported in supplemental table S5 (https://links.lww.com/ALN/D687).
Discussion
We developed a rigorous PONV phenotype by combining a validated observational database with a standardized and validated quality-improvement metric. We then integrated this perioperative phenotype with genotypes from two institutional biorepositories to enable genome-wide association in a 64,523-patient discovery cohort. This genome-wide association study identified 46 variants exceeding the P < 1 × 10−5 threshold, occurring with minor allele frequency greater than 1% and demonstrating concordant effects in both cohorts (table 1). We repeated the genome-wide association study in a 25,262-patient derivation cohort and used the associations across 1,048,398 genetic variants (minor allele frequency greater than or equal to 0.01) to calculate a polygenic score for PONV for 61,503 surgical procedures in the validation cohort.
When controlling for a diverse set of covariates (basic, known clinical risk factors, and full) and considering primary and secondary diagnostic criteria for PONV, the effect size remained significant and exhibited very consistent effect size (adjusted odds ratio, 1.019 to 1.029). This implies that a patient with a polygenic risk score more than 1 SD above the mean, has 2 to 3% greater odds of developing PONV when compared to the baseline population (polygenic score for PONV less than 1 and greater than −1). This increase is much smaller than the increase associated with having prior PONV/motion sickness (55%), having a history of migraines (17%), or being female (83%; table 3) and is not likely of clinical significance.
The addition of polygenic scores for PONV improved the discrimination of the known clinical risk factors model but did not improve the discrimination of the basic or full clinical models. This could be because the area under the receiver operating characteristic curve is an insensitive marker of incremental model improvement when the baseline model performs well,52,53 the effect size of polygenic score for PONV for predicting PONV is small, or clinical markers in the model already capture the risk that is also captured by genetics (clinical markers are heritable and are the ones underlying the genetic association with PONV). Although addition of a polygenic score for PONV led to statistically significant improvement in discrimination (c-statistic, 0.616 vs. 0.613; P = 0.028) in the model consisting of known clinical risk factors, this improvement is not clinically meaningful. Furthermore, comparing c-statistics across models is inherently a low-power procedure. The full clinical model also had improved discrimination compared to the basic or known clinical risk factors models, suggesting that a more diverse set of clinical variables leads to improved prediction.
Genome-wide Significance
On primary analysis, we did not identify any variants reaching genome-wide significance (P < 5 × 10−8). Analysis of secondary outcome identified one low-frequency (minor allele frequency = 0.002) variant (rs760998109) on chromosome 2 (P = 4.4 × 10−8); however, this variant was only present in the Michigan Medicine data set and could not be assessed in our meta-analysis. Although numerous loci reached the “suggestive” threshold (P < 1 × 10−5),40 these should not be viewed as significant novel associations, because adoption of less stringent thresholds leads to increased rate of false discovery.40
Concordance with Expected Results
Notably, no previously described associations were confirmed at the replication threshold (P < 1 × 10−5). A mutation in exon 1 of the μ-opioid receptor OPRM1 (rs1799971), leading to substitution of the amino acid asparagine for aspartate at position 40, was the only previously described variant reaching significance at P < 0.05 threshold, although it was not significant at the P < 1 × 10−5 threshold selected a priori for replication of previously described associations.57,58 Our study has a much greater statistical power than these previous publications, which employed candidate gene approach (64,523 in our study compared to 454,1 306,59 171,11 and 18913 ).
The five genes receiving the highest priority score using the polygenic priority scoring algorithm (MAPT, PLA2G4A, THRB, NRP1, and IL10) were all previously shown to be associated with nausea in the Comparative Toxicogenomics Database data set.60 We did not find association between any of the dopamine receptor genes (most notably DRD2)15 or serotonin receptor genes (most notably HTR3A and HTR3B),11,13 which have been previously reported and have a strong mechanistic link. This may result from key differences in preoperative or intraoperative antiemetic administration between our study methodology (left to provider discretion with 98.9% receiving and patients routinely receiving multiple intraoperative antiemetics) and previous studies (which excluded antiemetic prophylaxis).11,13,15,19 If providers use dopamine antagonists or serotonin antagonists and are successful at preventing PONV, this might give the appearance that the associated receptor genes are not involved in PONV. Therefore, routine use of intraoperative antiemetics (defined as administration before PACU) in our study population may select for genes that are not antagonized by antiemetics.
Our multivariable regressions confirmed the strong association between female sex, nonsmoking status, prior history of motion sickness or PONV, and opioid administration and the PONV phenotype. Other variables like, age, duration of anesthesia, history of migraine, and intraoperative opioid administration found to be significant in our full clinical model agreed with a large meta-analysis of 22 studies and 95,154 patients.61 Association between history of anxiety62 has been also previously demonstrated but was not confirmed on our primary analysis. Neuraxial anesthesia was associated with a lower incidence of PONV, even when adjusting for opioid administered (adjusted odds ratio, 0.836; 95% CI, 0.738 and 0.947; P < 0.001).63
Clinical Implications
Polygenic score is associated with PONV and improves discrimination in the known clinical risk factors model. These findings are of statistical but unlikely clinical significance. Furthermore, polygenic score for PONV no longer improves discrimination when incorporated into a full clinical model. Finally, the influence of polygenic score for PONV is negligible compared with established clinical risk factors. Incorporating a polygenic score for PONV led to improved reclassification of 4 to 5% of patients, which although statistically significant may be of limited clinical significance.
Study Limitations
Our study is constrained by the major limitations inherent to observational studies. As such, we cannot posit causation between any of our covariates and PONV. Although efforts have been made by the Anesthesiology Performance Improvement and Reporting Exchange team to ensure reliability of clinical factors and standardization of PONV reporting, overall accuracy is limited by documentation compliance. Furthermore, we do not specifically account for changes in practice pattern over the 6-yr inclusion period.
Although the c-statistic for the known clinical risk factors model with and without polygenic score may have a statistically significant difference (0.616, 0.609 to 0.623 vs. 0.613, 0.607 to 0.620; P = 0.028), the practical and clinical significance is limited. Another inherent limitation of our data set is the lack of genetic diversity (92% European ancestry). Although the data within Michigan Genomics Initiative and BioVU are reflective of the patient population served by Michigan Medicine and Vanderbilt University Medical Center and lack of ancestral diversity is commonly encountered problem across many biobanks,64 integration of our data with global biobanks composed of individuals from more diverse ancestries will enable even more equitable, impactful research. Excluding patients with ancestry other than white European due to the relatively small number of these patients in the Vanderbilt and Michigan cohorts prevented us from studying the potential association of genetic predictors in these groups of surgical patients. For example, a recent study demonstrated differential rates of antiemetic prophylaxis based upon patient race,65 raising important questions whether social or biologic factors might contribute to this disparity—that can be addressed through increasingly diverse cohorts.66 Furthermore, analysis of patients at two tertiary academic centers in the United States is likely to produce different results compared to international or community centers.
Another limitation is that the severity of PONV is not stratified by the phenotype definition. Furthermore, routine use of intraoperative antiemetics in our study population may select for genes that are not antagonized by antiemetics. Alternatively, the routine use of intraoperative prophylactic antiemetics may mask the genetic influence on PONV; however, testing this hypothesis would be challenging due to ethical ramifications of withholding intraoperative prophylactic antiemetics. Future studies may segregate the worst cases from more mild cases, with the hope of obtaining a stronger genetic linkage. Models created from the more restrictive primary outcome had higher discrimination compared to models from a secondary (composite) outcome; these models correctly reclassified a higher number of cases (net reclassification index for primary outcome was greater than net reclassification index for composite outcome), and there was greater concordance between the Michigan Medicine and Vanderbilt University Medical Center data in the primary outcome compared to secondary outcome. Antiemetic administration is likely documented with higher fidelity compared to subjective findings of PONV. As such, our analysis of the broader, secondary outcome is limited by decreased specificity, whereas analysis of the more restrictive, primary outcome sacrifices sensitivity.67 The optimal balance between these two criteria remains unresolved; however, our findings suggest that future studies using Multicenter Perioperative Outcomes Group data might focus only on cases with documented nausea or vomiting in the PACU, instead of also including cases requiring antiemetics in the PACU.
A final limitation surrounds the way (1) nitrous oxide administration (yes/no), (2) intraoperative opioid (oral morphine equivalents · kg−1 · min−1), and (3) inhalational anesthetic (percentage of total case) are defined within the full clinical model. Administration of 10 min of nitrous oxide at the end of the case compared to 6 h of 70% nitrous oxide (alternatively: 100 μg fentanyl at induction vs. 10 mg morphine at the end of the case) would be expected to have very different influence on a patient’s overall PONV risk; however, current models treat them equivalently. Additionally, all volatile anesthetics (isoflurane, sevoflurane, desflurane, halothane, and enflurane) are treated as equivalent, despite a potentially differential impact on PONV risk. Despite this limitation, covariates for these three anesthetic factors were selected to match standardized definitions (and data extraction algorithms) within the Multicenter Perioperative Outcomes Group known as “phenotypes.” This ensures (1) validation of data quality, (2) generalizability across Multicenter Perioperative Outcomes Group sites and studies, and (3) ease of data extraction.
Conclusions
Standardized polygenic risk was associated with PONV in all three of our models, but the genetic influence was smaller than exerted by clinical risk factors. Furthermore, the use of a polygenic risk score does not meaningfully improve discrimination compared to clinical risk factors. We also noted significantly improved discrimination in our full model (compared to the known clinical risk factors model), suggesting that PONV risk stratification may benefit from inclusion of more diverse, nongenetic data, when extended to current anesthetic practice.
Appendix 1. Multicenter Perioperative Outcomes Group Perioperative Research Standard
The Multicenter Perioperative Outcomes Group perioperative research standard is defined by the following criteria:
Has valid anesthesia start time
Has valid anesthesia end time
Has a valid institution identification
Case duration longer than 15 min if anesthesia technique general = true
Case duration longer than 5 min if anesthesia technique general = false
Has an actual or predicted anesthesia Current Procedural Terminology code
Has age data
Has sex data
Has valid American Society of Anesthesiologists class
Has baseline blood pressure mean
Has at least one intraoperative medicine administered
Has at least one International Classification of Diseases-9/10 discharge diagnosis
Has at least one creatinine or hematocrit within 365 days before/after surgery
Appendix 2. Characterization and Genetic Quality Control
As part of the independent data flow for both individual sites, sample and genotype quality control included assessment of call rates, sex check, cryptic relatedness, single-nucleotide polymorphism missingness and the Hardy–Weinberg equilibrium. The standardized, predetermined quality control metrics of each respective biobank were retained, leading to minor differences in quality control metrics between BioVU and Michigan Genomics Initiative. BioVU excludes single-nucleotide polymorphisms with an overall call rate of less than 95%, high sample missingness (greater than 5%), or significant deviation from Hardy–Weinberg equilibrium (P < 10−6).31 The Michigan Genomics Initiative excludes single-nucleotide polymorphisms with an overall call rate of less than 99%, high sample missingness (greater than 1%), or significant deviation from Hardy–Weinberg equilibrium (P < 10−4) within each array.32,33 Principal components were generated using the Eigenstrat method in the EigenSoft version 7.2.1 program to control for population stratification. Before meta-analysis and polygenic risk calculation, we harmonized the Michigan Genomics Initiative and Vanderbilt University Medical Center summary statistics as a quality control and standardization step (i.e., matching reference or alternate alleles, removing single-nucleotide polymorphisms with invalid alleles and nonpalindromic strand flips). To do this, we used EasyQC34 to match BioVU and Michigan Genomics Initiative summary statistics to the Genome Aggregation Database (version 3.1.2) with FILTER=PASS for the non-neuro and non-Finnish European population. This removed nonpalindromic strand flips, invalid alleles (for example, chr1:855480:CAT:GAT) and sites with absolute allele frequency greater than 0.2 compared to the Genome Aggregation Database. Additionally, the variants used for the polygenic risk score were selected from imputed data with R2 greater than 0.3, hard call threshold ±0.1, and minor allele frequency greater than 1%.35 Only chromosomes 1 to 22 are included in the analysis.
Appendix 3. Medication List of Antiemetics
The following antiemetic medication list was defined by the Anesthesiology Performance Improvement and Reporting Exchange (ASPIRE) quality improvement metric.
Class: 5-Hydroxytryptamine receptor antagonists
Ondansetron
Dolasetron
Granisetron
Palonosetron
Anticholinergics
Scopolamine patch
Scopolamine
Butylscopolamine
Antihistamines
Dimenhydrinate
Diphenhydramine
Meclizine
Butyrophenones
Droperidol
Haloperidol
Neurokinin-1 receptor agonists
Aprepitant
Fosaprepitant
Phenothiazines
Promethazine
Prochlorperazine
Steroid
Dexamethasone
Methylprednisolone
Prokinetic
Metoclopramide
Other
Propofol (infusion only)
Note: Amisulpride was not included in the standardized ASPIRE metric utilized for this project, although the medication may be used at Michigan or Vanderbilt.
Appendix 4. Genotyping Protocol and Characterization of Genetic Data
Michigan Cohort: Michigan Genomics Initiative Freeze 5
Subjects were genotyped using Illumina Infinium CoreExome-24 bead or global screening arrays (Illumina Inc., USA) and imputed using the Trans-Omics for Precision Medicine panel. Variants with R2 greater than or equal to 0.3, minimum minor allele count greater than or equal to 4, and minor allele frequency of greater than or equal to 0.01% were selected, resulting in more than 50 million imputed variants after quality control and filtering. The European population was defined using the majority Admixture global ancestry fraction (Q value).
Vanderbilt Cohort
Subjects were genotyped using DNA extracted from discarded plasma samples with the Illumina Infinium expanded multiethnic genotyping array (MEGAEX, Illumina Inc.). Imputation of the Vanderbilt MEGAEX genotyped data set was performed using the Minimac4 and 1000 Genomes Phase 3 version 5 with the Michigan Imputation Server. Variants with R2 greater than or equal to 0.3 and minor allele frequency greater than or equal to 1% were selected, resulting in more than 21.1 million overall imputed variants available and 11.7 million imputed variants available for the acute kidney injury genome-wide association study with white European ancestry after quality control and filtering.
Ancestry Selection
We used supervised learning mode in Admixture to calculate ancestry coefficients for our data set.37 The 1000 Genomes Phase 3 data, which contains 2,504 reference individuals (of 5 ancestries), were used as training samples in the model. The five ancestry fractions (K_EUR, K_AFR, K_AMR, K_EAS, and K_SAS) are produced, and we used a threshold of K_EUR of greater than 0.9 to select the data set of European ancestry.
Appendix 5. Site-specific Protocols for Genome-wide Association Study
All patient and genetic data were retained at the local site (either Michigan Medicine for Michigan Genomics Initiative or Vanderbilt University Medical Center for BioVU). Results from the two genome-wide association studies were then meta-analyzed. Standard workflow for each site was maintained leading to minor differences (software version, imputation methodology, and quality control) between the two sites. The specific genome-wide association study methodology at each site is detailed below:
For the Michigan Genomics Initiative population, we implemented the Scalable and Accurate Implementation of Generalized mixed model (SAIGE) version 0.45.1 and controlled for age, genotype inferred sex, genotyping array, and the first four principal components. Genome-wide association studies were performed on genetic data from a European cohort in Michigan Genomics Initiative Freeze 5. The European population was defined using the majority Admixture global ancestry fraction (Q value.) Before fitting the null model in SAIGE step 1, we used PLINK to linkage disequilibrium (LD) prune directly assayed genotypes by squared correlation greater than 0.5 within a walking window of 500 variants with a step length of 5 variants. LD pruning was done for the cohort in each PONV phenotype independently. For association tests, we used Trans-Omics for Precision Medicine Imputed genotypes and excluded variants with low imputation quality (Rsq less than 0.3) or very rare minor alleles (minor allele frequency less than 0.0001 or minor allele count less than 4).
For the BioVU population, we implemented the SAIGE version 1.0.0 and controlled for age, sex, and the first four principal components. Genome-wide association studies were performed on genetic data from a European cohort in BioVU. The European population was defined using Admixture and with a global ancestry fraction (Q value) greater than or equal to 0.9. Before fitting the null model in SAIGE step 1, we used PLINK to LD prune directly assayed genotypes by squared correlation greater than 0.5 within a walking window of 500 variants with a step length of 5 variants. LD pruning was done for the cohort in each PONV phenotype independently. For association tests, we used imputed genotypes, and imputation was performed using Minimac4 and 1000 Genomes Phase 3 version 5 reference panel on the Michigan Imputation Server and excluded variants with low imputation quality (Rsq less than 0.7) or low allele frequency (minor allele frequency less than 0.01).
Appendix 6. Harmonization of Summary Statistics and Meta-analysis
We harmonized genome-wide association study summary statistics from BioVU and Michigan Genomics Initiative data to the Genome Aggregation Database43 version 3.1.2 as a reference using EasyQC44 version 23.8. Briefly, we harmonized the allele order and strand convention of validation summary statistics to reference and excluded sites that were not in reference, had invalid alleles, had mismatched alleles, or had an allele frequency in validation that differed from the reference population “non_neuro_nfe” by more than 0.2. We then used METAL in inverse variance weighted mode to meta-analyze variants represented in both the Michigan Genomics Initiative and BioVU cohorts.45
Acknowledgments
The authors acknowledge the help of Kathryn Buehler, M.S., R.N., C.P.P.S. (Clinical Program Manager), and Meridith Bailey, M.S.N., R.N. (Quality Initiative Coordinator) at the Multicenter Perioperative Outcomes Group/Anesthesiology Performance Improvement and Reporting Exchange (Ann Arbor, Michigan) for assistance curating the PONV phenotype. The authors acknowledge Baorong Shi, M.S. (Department of Anesthesiology, University of Michigan Health System, Ann Arbor, Michigan), for contributions in data acquisition and electronic search query programming for this project. The authors also acknowledge Justin Ortwine, B.S. (Department of Anesthesiology, University of Michigan Health System), and Erin O. Kaleba, B.A., M.P.H. (Data Office for Clinical and Translational Research, University of Michigan Medical School, Ann Arbor, Michigan), for help with acquisition of Michigan Genomics Initiative data. The authors thank Anisa Driscoll and Anita Pandit from the Precision Health Initiative (Michigan Medicine, Ann Arbor, Michigan) for technical expertise with genetic data acquisition and imputation. The authors acknowledge the Michigan Genomics Initiative participants, Precision Health at the University of Michigan, the University of Michigan Medical School Central Biorepository, and the University of Michigan Advanced Genomics Core for providing data and specimen storage, management, processing, and distribution services and the Center for Statistical Genetics in the Department of Biostatistics at the School of Public Health for genotype data curation, imputation, and management in support of the research reported in this publication. The authors thank the clinicians, staff, and study participants from the Michigan Genomics Initiative.
Research Support
Supported by the Department of Anesthesiology, University of Michigan Medical School (Ann Arbor, Michigan) and the Department of Anesthesiology, Vanderbilt University Medical Center (Nashville, Tennessee). Dr. Douville received support from National Institute of Diabetes and Digestive and Kidney Diseases (Bethesda, Maryland) grant No. 1 K08 DK131346-01 and a Foundation for Anesthesia Education and Research (Schaumburg, Illinois) Mentored Research Training Grant. Ms. Bastarache received funding from National Library of Medicine (Bethesda, Maryland) grant No. R01-LM010685-11. Dr. Mathis received grant funding support from the National Institutes of Health, National Heart, Lung, and Blood Institute and National Institute of Diabetes and Digestive and Kidney Diseases.
Competing Interests
Dr. C. B. Douville is a consultant to Exact Sciences (Madison, Wisconsin) and the founder of Belay Diagnostics (Chicago, Illinois) and is compensated with income and equity. Exact Sciences is a molecular diagnostics company specializing in the detection of early-stage cancers. Belay Diagnostics is a platform for the detection of brain and spinal cord cancers using cerebrospinal fluid. These companies, as well as other companies, have licensed previously described technologies. Licenses to these technologies are or will be associated with equity or royalty payments to the inventors, as well as Johns Hopkins University (Baltimore, Maryland). The terms of all of these arrangements are being managed by Johns Hopkins in accordance with its conflict-of-interest policies. Dr. C. B. Douville’s involvement with both companies is unrelated to the scope of the research presented in this article. Dr. Willer is currently employed by Regeneron Pharmaceuticals Inc. (Tarrytown, New York) and holds stock and options. Dr. Wu is currently employed by Regeneron Pharmaceutical Inc. Dr. Mathis received funding paid to the University of Michigan from Chiesi USA (Cary, North Carolina) for work unrelated to the scope of the research presented in this article. The other authors declare no competing interests.
Reproducible Science
Full protocol available at: ndouvill@med.umich.edu. Raw data available at: ndouvill@med.umich.edu.
Supplemental Digital Content
Supplemental Information, https://links.lww.com/ALN/D687
Figure S1: Flow diagram depicting derivation of study cohort
Figure S2. Quantile–quantile plots
Table S1. List of available clinical variables
Table S2. Descriptive statistics on BioVU cohort
Table S3. Comparison to previous genome-wide association study results
Table S4. Genes prioritized by polygenic priority scoring
Table S5. Net reclassification index for each model