Opioid agonist and antagonist use and the gut microbiota: associations among people in addiction treatment

Murine models suggest that opioids alter the gut microbiota, which may impact opioid tolerance and psychopathology. We examined how gut microbiota characteristics related to use of opioid agonists and antagonists among people receiving outpatient addiction treatment. Patients (n = 46) collected stool samples and were grouped by use of opioid agonists (heroin, prescription opioids), antagonists (naltrexone), agonist–antagonist combinations (buprenorphine–naloxone), or neither agonists nor antagonists within the month before enrollment. We sequenced the V4 region of the 16S rRNA gene using Illumina MiSeq to examine how alpha diversity, enterotypes, and relative abundance of bacterial genera varied by opioid agonist and antagonist exposures. Compared to 31 participants who used neither agonists nor antagonists, 5 participants who used opioid agonists (without antagonists) had lower microbiota diversity, Bacteroides enterotypes, and lower relative abundance of Roseburia, a butyrate producing genus, and Bilophila, a bile acid metabolizing genus. There were no differences in gut microbiota features between those using agonist + antagonists (n = 4), antagonists only (n = 6), and neither agonists nor antagonists. Similar to murine morphine exposure models, opioid agonist use was associated with lower microbiota diversity. Lower abundance of Roseburia and Bilophila may relate to the gut inflammation/permeability and dysregulated bile acid metabolism observed in opioid-exposed mice.


Supplementary Methods Methods Overview
We consented 46 patients aged 18-60 years who were receiving treatment in an outpatient addiction treatment facility to participate in a gut microbiota and substance use study ( Supplementary Fig. S1). We surveyed participants about the substances they used in the 30 days prior and reviewed participant medical records to characterize their opioid agonist and antagonist use. Participants submitted a stool sample and we sequenced the V4 hypervariable region of the 16S ribosomal ribonucleic acid (rRNA) gene to characterize the gut microbiota. Each participant's opioid agonist and antagonist use patterns from survey and medical record data were summarized as agonist only (Ag), combined agonist-antagonist (AgAt), antagonist only (At), or neither agonist nor antagonist (N). We compared gut microbiota diversity, enterotypes, and genera relative abundance between participants in these four opioid agonist and antagonist groups. The study was approved by the Institutional Review Board at the University of Michigan (HUM00113964).

Substance Use Assessment for Microbiota Study Eligibility
We used an abbreviated version of the Alcohol, Smoking, and Substance Involvement Screening Test (ASSIST) to assess past 30-day substance use. ASSIST is a 12-item questionnaire that asks about use of street opioids, opioid pain medications (with or without a prescription), methadone or buprenorphine-naloxone (with or without a prescription), tobacco, alcohol, cannabis, cocaine (including crack), amphetamines, inhalants, sedatives or sleeping pills, and hallucinogens 1,2 . Participants who only endorsed tobacco were not eligible for the microbiota study without concurrent reported use of another substance. Participants reporting past 30-day use (nonmedical or as prescribed) of medications for opioid use disorder (buprenorphinenaloxone or methadone) were eligible, regardless of past 30-day use of other substances. We used a modified version of the Current Opioid Misuse Measure (COMM) to capture self-reported misuse of prescribed opioids 3 . The COMM is an eight-item scale that assesses self-reported opioid misuse not captured by ASSIST 3 . Participants with a score >0 for questions 4 and 6-8 were eligible based on self-reported prescription opioid misuse regardless of their past 30-day substance use from the ASSIST 3 .

Opioid Agonist and Antagonist Use
We used a combination of the ASSIST and medical record review to identify any opioid agonist use. We defined opioid agonist use as a binary indicator for either of the following: 1) self-reported opioid use (heroin, methadone, buprenorphine, or prescription opioids used as prescribed or not as prescribed) in the 30 days before study enrollment on the modified ASSIST, or 2) buprenorphine-naloxone use during the day of sample collection in the medical record. We supplemented the ASSIST data on buprenorphine use with data from the medical record because survey questions about medications for opioid use disorder were added to the modified ASSIST during January 2017 and were not available prior. Of note, one participant who was prescribed opioids to manage pain was included in the group of participants who used opioid agonists.
We also assessed opioid antagonist use. For the purpose of this study, we defined opioid antagonist use as either of the following according to the medical record: 1) buprenorphinenaloxone use during the day of sample collection, or 2) naltrexone use during the day of sample collection. We included any formulation of either opioid antagonist. Finally, we created a categorical variable describing the overlap in opioid agonist and antagonist use (i.e., agonist only, combined agonist-antagonist, antagonist only, and neither agonist nor antagonist).

Depression, Anxiety, and Cravings to Use
We summarized scores from the Patient Health Questionnaire-9 (PHQ-9), a reliable and

Microbiota Study Population Characteristics
We summarized average age, self-identified gender (female, male, and other), race (black, white, other, or multiple races), ethnicity (Hispanic vs. non-Hispanic), and self-reported days in SUD treatment at the study site. We used the ASSIST to describe alcohol use during the 30 days before completing the substance use survey 1 . Finally, we summarized self-reported antibiotic use during the week of sample collection.

Dietary Fiber Intake Estimation
Dietary fiber intake has previously been associated with the gut microbiota, transit time, and stool water content 9-12 . We therefore calculated predicted past-month dietary fiber intake (grams per day) from food frequency questions completed at enrollment in the microbiota study using validated predictive models developed from the National Health and Nutrition Examination Survey (NHANES) 13 . Upon enrollment in the microbiota study, participants completed a validated food frequency questionnaire about general past month intake of 25 food and drink items 14 . Due to a skipped question on the survey, participants did not report their cheese intake. We therefore assigned participants the survey weighted mean values of cheese intake using age and race-specific values from the 2009-2010 cycle of the NHANES, which included the food frequency questionnaire used in our study. We did not use sex-specific values as cheese intake did not differ by sex in NHANES. We calculated predicted past month dietary fiber intake (predicted grams per day) from reported dietary item intake and validated predictive models developed from the National Health and Nutrition Examination Survey 13 .

Stool Sample Collection
We based our stool collection protocols on those from shaking, and stored the sample at room temperature for up to two days before returning the sample to a research assistant. Previous research supported that RNAlater TM preserved the composition of the stool bacterial community at room temperature for up to three days after collection 16,17 . Stool samples were then frozen in 1 mL aliquots at -80C.

DNA Extraction and Illumina MiSeq Sequencing
After all participants were enrolled, we thawed a 1 mL aliquot of each stool sample, centrifuged at 10,000xg, and resuspended in 1 mL 1X phosphate buffered saline (Gibco, Thermo Two samples with fewer than 1,000 reads were removed from further analysis; however, we were able to retain these two participants for analysis using their duplicate sequenced sample. We verified that all mock communities resembled their known compositions (data not shown)

Genus Relative Abundance and Diversity Metrics
We calculated the relative abundance of genera in each sample (i.e., the number of sequencing reads from each genus divided by the total number of sequencing reads per sample).
For alpha diversity analyses, we normalized sequencing depth across samples using rarefaction 24,25 . Briefly, the oligotype table was sampled without replacement to 90% of the minimum sequencing depth across all 129 samples sequenced as part of the study (i.e., sampled to a sequencing depth of 6,338 reads). Rarefaction was repeated 100 times, and alpha diversity metrics were averaged across the 100 replicates. Because the application of rarefaction in microbiome research is the subject of some debate 26 , we also analyzed alpha diversity metrics on the non-rarefied oligotype table as a sensitivity analysis. Alpha diversity metrics included Shannon diversity and the Chao1 Index, which summarized within-sample oligotype diversity (number and evenness of oligotypes) and richness (number of oligotypes), respectively. Further, we visualized between-sample (beta diversity) differences using principal component analysis of the Aitchison Distance metric, which visualizes Euclidean distance using a centered log ratio transformation of the oligotype table, an approach consistent with the compositional nature of the data (described further below) 27 .
Microbiota sequencing data are compositional; the total number of sequences from each sample is bound by an arbitrary depth of sequencing coverage that does not reflect the true abundance of bacteria in the participant's gut 28 . Thus, taxa comparisons between samples must be analyzed on the relative (not absolute) scale 28 . Failure to analyze these data as compositional results in spurious correlations between taxa and biases associations of taxa abundance with other taxa and covariates 28 . We therefore applied three analytic approaches to compare taxa distributions by covariates of interest that accounted for the compositional nature of the data: the Aitchison Distance described above, and Dirichlet multinomial mixture modeling and ALDEx2, described further below [29][30][31] .

De Novo and Reference-Based Enterotyping
Enterotyping distills highly dimensional microbiota taxa data into clusters (i.e., groups of samples exhibiting similar taxa distributions) by leveraging the information from the covariance structure of taxa 9,30,32 . The high dimensionality of microbiome data (number of taxa) poses a challenge during analysis, as testing for associations of single taxa with metadata requires many tests and results can be difficult to interpret without acknowledging co-occurring changes in other taxa. Enterotyping addresses these challenges through data reduction and summary by clustering samples with similar bacterial profiles into groups based on the distribution of taxa in samples 9 .
We applied two clustering approaches that classified each stool sample's genus-level read counts into enterotypes 9,30 . First, we used Dirichlet multinomial mixture (DMM) models, an extension of latent profile analysis adapted for microbiota data by Holmes et al. 30 . Like traditional latent profile analysis, this technique recovers unobserved (i.e., latent) subgroups based on joint distributions of bacterial genera. DMM-based enterotypes were assigned de novo (i.e., based on the data and participants included in our study) using data from all 129 samples collected during the study to standardize DMM assignments across all samples. We compared model fit using the Laplace approximation of negative log models for DMM models with one to five enterotypes and chose the number of enterotypes that optimized model fit (i.e., minimized the Laplace approximation) 30 . We then assigned each sample to its most likely enterotype based on posterior probabilities of enterotype assignment. All 46 samples used in the main analysis had posterior probabilities ³90% and were therefore all assigned to an enterotype (minimum posterior probability: 95.3%). We examined the average relative abundance of all genera and the top 20 taxa that were most influential in distinguishing between enterotypes to summarize taxa distributions representative of each enterotype. We uploaded genus relative abundance data from our study to http://enterotypes.org and obtained enterotype assignments and a binary variable indicating whether each sample had similar genera to the observed patterns in reference samples from HMP and MetaHIT. Two of 46 samples in the main analysis were not comparable to reference samples based on this indicator and were therefore assigned as "missing" for their assignment.

Genus Differential Abundance by Opioid Agonist and Antagonist Use
We compared abundance of specific genera by opioid use status and other covariates of interest using ALDEx2, an analysis of variance-like tool for compositional data 31 . For each genus, ALDEx2 inferred absolute abundance given the observed abundance matrix over 1,000 Monte-Carlo simulations from a Dirichlet distribution. To account for the compositional structure of the data, genera abundances for each sample and simulation were transformed to centered log ratios (i.e., log of the ratio of taxai abundance for sample j divided by the geometric mean for the abundance of all taxa [i=1:K ] in sample j). ALDEx2 then calculated each genus' median centered log ratio by group (e.g., opioid agonist use vs. no use) for each simulation.
Within-group variability in each taxa's log ratios reflected sampling variation whereas betweengroup differences represented the biological variation of interest (e.g., by opioid agonistantagonist use status). These were used to calculate an effect size for each genus as the median difference in centered log ratios between groups across all simulations divided by the median of the largest detected difference within groups for each condition (e.g., the median of two items: 1) the largest difference in centered log ratios for taxai among opioid agonist use group across 1,000 simulations and 2) the largest difference in centered log ratios for taxai among people not using opioid agonists across 1,000 simulations). Statistical significance for each effect size was summarized as a p-value corrected for the false discovery rate (FDR) using the Benjamini-Hochberg procedure. We summarized genera centered log ratios and relative abundance for all samples as heat maps (Fig. 3) for taxa with FDR corrected p-values<0.05 identified from the ALDEx2 Wilcoxon rank sum test.

Supplementary Fig. S1 Inclusion Criteria and Opioid Use among 46 Participants Enrolled from an Outpatient Addiction Treatment Facility into a Study of the Gut Microbiota and Opioid Use, 2016-2017
We included 46 participants who were eligible and provided informed consent for both study stages and provided stool samples in the present analysis. Among the 9 using opioids, 5 used only agonists (prescription opioids or heroin) and 4 used an agonist-antagonist combination during the time their sample was provided. Among the 37 not using opioids, 6 used an opioid antagonist (naltrexone) and 31 were not exposed to neither opioid agonists nor antagonists during the time they provided their stool sample.
• 5 lost to follow-up or declined to participate after consent  We fit Dirichlet multinomial mixture models with one to five enterotypes. Model fit was optimized by a three enterotype model, which minimized the Laplace approximation. We identified two Bacteroides enterotypes, one with increased Faecalibacterium and a second with increased Clostridium cluster XIVa. The third enterotype was dominated by Prevotella. As such, Prevotella and Bacteroides were the first and second-most influential genera in assigning enterotypes, respecitively. We compared alpha diversity using the Shannon diversity and Chao1 metrics by past 30-day alcohol use (top) and dietary fiber (bottom). Shannon diversity was marginally lower among participants who did not use alcohol in the past 30 days (Wilcoxon rank sum p=0.052). Gut microbiota richness was positively associated with fiber intake and the Pearson correlation (r=0.35) between richness and fiber intake reached statistical significance (p=0.02). Other comparisons by alcohol use and fiber intake, including Spearman correlation coefficients for fiber and alpha diversity metrics, did not reach statistical significance. 12  p=0.04 (Fig. 1) p=0.05 p=0.08