Single cell transcriptomics reveals opioid usage evokes widespread suppression of antiviral gene program

Single cell transcriptomics reveals opioid usage evokes widespread suppression of antiviral gene program

PBMCs from opioid-dependent individuals


Frozen vials of PBMCs prepared from the fresh blood of opioid-dependent (mostly heroin dependent), and non-dependent neighborhood control individuals were collected in the Comorbidity and Trauma Study (CATS)37,38 and subsequently obtained from the biorepository of National Institute on Drug Abuse (NIDA, Rockville MD). We evaluated 14 age-matched subjects ranging from 24 to 45 years in age, with an equal number of male and female subjects (Supplementary Table 6). Cases were recruited from opioid replacement therapy (ORT) clinics in the greater Sydney, Australia region, while controls were recruited from areas in close proximity to ORT clinics (neighborhood controls). Cases and controls were required to be English speakers 18 years of age or older. Cases were participants in ORT for opioid dependence while controls were excluded for recreational opioid use 11 or more times lifetime. All subjects provided written informed consent37,38. All samples were stripped of personally identifying information and assigned sample ID numbers prior to receipt. No further ethical oversight was required from the Boston University IRB following de-identification of the procured samples.


scRNA-seq of LPS-treated patient PBMCs


Frozen PBMCs isolated from the blood of opioid-dependent and non-dependent neighborhood individuals were revived, and live cells were isolated via fluorescently activated cell sorting (FACS) using a Sony SH800 cell sorter and a live/dead cell stain (LIVE/DEAD Fixable Green Cell Stain Kit, for 488 nm Excitation, Thermo Fisher, L34969). The FACS gating strategy used for the isolation of live PBMC is illustrated in Supplementary Fig. 23. Dilutions were prepared from all 14 samples at a concentration of 1000 cells/μL as outlined in the Chromium Single Cell 3′ Reagent Kit v2 User Guide (10X Genomics, CG00052 Rev.B), and 7000 cells per sample were used to perform the droplet-based Chromium Single Cell 3′ scRNA-seq method (10X Genomics, Chromium Single Cell 3′ Library and Gel Bead Kit, Cat# PN-120237). In total, 200,000 cells from 6 of the 14 samples (three dependent and three non-dependent) were plated into a non-tissue culture-treated 96-well plate in a leukocyte-supporting complete RPMI medium (10% HI-FBS, 1% l-glutamate, 1% NEAA, 1% HEPES, 1% sodium pyruvate, 0.1% B-mercaptoethanol). Lipopolysaccharide (LPS) (Invivogen, LPS-EK Ultrapure, Cat# tlrl-pekpls) was then added to a final concentration of 100 ng/mL and the cells were incubated at 37 ˚C for 3 h. Cells were then collected, washed, and diluted to 1000 cells/μL before being used to perform the 10X Genomics Chromium Single Cell 3′ method as outlined in the Single Cell 3′ Reagent Kit v2 User Guide. Briefly, 20 μL of 1000 cells/μL PBMC suspension from each subject/condition were combined, and 33.8 μL of cell suspension (total cell number = 33,800) was mixed with 66.2 μL of RT reaction mix before being added to a chromium microfluidics chip already loaded with 40 μL of barcoded beads and 270 μL of partitioning oil. The chip was then placed within the chromium controller where single cells and barcoded beads were encapsulated together within oil droplets. Reverse transcription was then performed within the oil droplets to produce barcoded cDNA. Barcoded cDNA was isolated from the partitioning oil using Silane DynaBeads (Thermo Fisher Scientific, Dynabeads MyONE Silane, Cat# 37002D) before amplification by PCR. Cleanup/size selection was performed on amplified cDNA using SPRIselect beads (Beckman-Coulter, SPRIselect, Cat# B23317) and cDNA quality was assessed using an Agilent 2100 BioAnalyzer and the high-sensitivity DNA assay (Agilent, High-Sensitivity DNA Kit, Cat# 5067-4626). Sequencing libraries were generated from cleaned, amplified cDNA using the 10X Chromium Kit's including reagents for fragmentation, sequencing adaptor ligation, and sample index PCR. Between each of these steps, libraries were cleaned and size selected using SPRIselect beads. Final quality of cDNA libraries was once again assessed using the Agilent BioAnalyzer High-Sensitivity DNA assay, and quality-confirmed libraries were sequenced using Illumina’s NextSeq platform. All reagents are listed in Supplementary Table 2.


scRNA-seq of IFNβ-treated patient PBMCs


Frozen PBMCs isolated from the blood of three opioid-dependent and three non-dependent neighborhood individuals were revived, and live cells were isolated via FACS using a Sony SH800 cell sorter and a live/dead cell stain (LIVE/DEAD Fixable Green Cell Stain Kit, for 488 nm Excitation, Thermo Fisher—L34969). The FACS gating strategy used for the isolation of live PBMC is illustrated in Supplementary Fig. 23. Live cells were plated into a non-tissue culture-treated 96-well plate in a leukocyte-supporting complete RPMI medium (10% HI-FBS, 1% l-glutamate, 1% NEAA, 1% HEPES, 1% sodium pyruvate, 0.1% B-mercaptoethanol) at a density of 200,000 cells per well. Twenty-two microliters of a 100 U/mL solution of IFNβ (Recombinant Human IFN-beta Protein, R&D Systems, Cat# 8499-IF-010) was then added to each well (for a final IFNβ concentration of 10 U/mL) and cells were incubated for 3 h at 37 °C. After treatment, all cells were collected, and an equal number of cells per patient sample were collected, washed, and each sample was “hashed” using unique oligonucleotide-barcoded antibodies15 (Supplementary Table 4) to track the cells’ well/condition of origin. Briefly, cells were suspended in Cell Staining Buffer (BioLegend, Cat# 420201) and blocked using Human TruStain FcX reagent (BioLegend, Cat# 422301). Cells were then incubated with 1 μg of TotalSeq antibodies (BioLegend, Cat# 3964601, 394603, 394605, 394607, 394609, 394611, 394613, 394615, 394617, 394619, 394623, 394625), washed with PBS, and filtered through 40 μM cell strainers (Bel-Art, Flowmi Cell Strainer, Cat# H13680-0040). Samples were then normalized to 1000 cells/μL, mixed in equal measure (20 μL each), and used to perform the Chromium Single Cell 3′ scRNA-seq method. Briefly, 20 μL of 1000 cells/μL PBMC suspension from each subject/condition were combined, and 33.8 μL of cell suspension (total cell number = 33,800) was mixed with 66.2 μL of RT reaction mix before being added to a chromium microfluidics chip already loaded with 40 μL of barcoded beads and 270 μL of partitioning oil. The chip was then placed within the chromium controller where single cells and barcoded beads were encapsulated together within oil droplets. Reverse transcription was then performed within the oil droplets to produce barcoded cDNA. Barcoded cDNA was isolated from the partitioning oil using Silane DynaBeads (Thermo Fisher Scientific, Dynabeads MyONE Silane, Cat# 37002D) before amplification by PCR. Cleanup and size selection was performed on amplified cDNA using SPRIselect beads (Beckman-Coulter, SPRIselect, Cat# B23317) and cDNA quality was assessed using an Agilent 2100 BioAnalyzer and the high-sensitivity DNA assay (Agilent, High-Sensitivity DNA Kit, Cat# 5067-4626). Sequencing libraries were generated from cleaned, amplified cDNA using the 10X Chromium Kit’s including reagents for fragmentation, sequencing adaptor ligation, and sample index PCR. Between each of these steps, libraries were cleaned, and size selected using SPRIselect beads. Final quality of cDNA libraries was once again assessed using the Agilent BioAnalyzer High-Sensitivity DNA assay, and quality-confirmed libraries were sequenced using Illumina’s NextSeq platform. Additional primers were included in the cDNA amplification step to amplify the TotalSeq oligonucleotide tags. During the post-amplification cleanup, supernatant containing amplified TotalSeq tags was collected and processed parallel to the standard 10X library fraction. All reagents are listed in Supplementary Table 2.


PBMCs from healthy individuals used in in vitro assays


Fifteen milliliters of fresh whole blood from healthy donors (Research Blood Components, Boston MA) were diluted 1:1 with warm PBS + 2% FBS, mixed, and gently layered atop 30 mL Ficoll-Paque density gradient medium (GE Healthcare, Ficoll-Paque PLUS, Cat# 17-1440) in 50 mL conical tubes. This process was repeated seven times in processing 100 mL of blood. Tubes were centrifuged for 20 min at 1200 × g to separate leukocytes from red blood cells and plasma. The leukocyte-containing buffy coat was carefully transferred into new tubes, washed with warm PBS + 2% FBS, counted, resuspended in DMSO, and aliquoted. Isolated cells were then stored in liquid nitrogen until later experimental use. Due to the anonymous nature of the procured whole-blood samples, no ethical oversight was required from the Boston University IRB for these samples. All reagents are listed in Supplementary Table 2.


Controlled substances


Solid morphine sulfate (Sigma Aldritch, Cat# M8777-25G) was obtained with approval and oversight from the controlled substances sub-office of the Boston University Department of Environmental Health and Safety. Aliquots of a 10 mM stock solution were prepared and stored for further use in experimentation. All reagents are listed in Supplementary Table 2


Morphine titration in PBMCs from healthy individuals


Normal PBMCs were revived in leukocyte-supporting complete RPMI medium and plated onto non-tissue culture-treated 96-well plates at a density of 2.0e5 cells/well (two wells per condition, 4.0e5 cells total). Cells were treated either with a mock treatment or a titration of morphine sulfate in RPMI complete medium (0, 10, 100 μM) for 24 h. At the end of the morphine incubation either medium or LPS (final concentration 100 ng/mL) was added to the wells and cells were incubated for further 3 h, at the end of which the cells were collected, washed, and processed for total RNA using the ZymoPure QuickRNA MiniPrep kit (Zymo Research, Cat# R1055). RNA samples were then used to perform RT-qPCR. All reagents are listed in Supplementary Table 2


RT-qPCR analysis


Total RNA was isolated from cells using the ZymoPure QuickRNA MiniPrep kit. cDNA was synthesized using ~50 ng of total RNA per sample (Thermo Fisher, SuperScript IV First-Strand Synthesis System, Cat# 18091200). Two microliters of cDNA per reaction was then used to perform qPCR (Fisher Scientific, PowerUp SYBR Greβen Master Mix, Cat# A25742) with primers against transcripts of the Interferon target gene ISG15 (IDT; Supplementary Table 3) We used primers against ActB (IDT, primer sequences in supplement) as a housekeeping gene control. All reagents are listed in Supplementary Table 2.


scRNA-seq of in vitro morphine treatment with healthy PBMCs


Cells were treated either with a mock treatment or 100 μM of morphine sulfate in RPMI complete medium for 24 h. At the end of the morphine incubation either medium or LPS (final concentration 100 ng/mL) was added to the wells and cells were incubated for additional 3 h, at the end of which the cells were collected and washed. After treatment, all cells were collected, and an equal number of cells per patient sample was subjected to cell hashing for scRNA-seq using 1 μg of TotalSeq antibodies (Supplementary Table 5). Hashtagged cells were then washed, diluted to 1000 cells/μL, and pooled before being used to perform the 10X Genomics Chromium Single Cell 3′ method. Briefly, 20 μL of 1000 cells/μL PBMC suspension from each subject/condition were combined, and 33.8 μL of cell suspension (total cell number = 33,800) was mixed with 66.2 μL of RT reaction mix before being added to a chromium microfluidics chip already loaded with 40 μL of barcoded beads and 270 μL of partitioning oil. The chip was then placed within the chromium controller where single cells and barcoded beads were encapsulated together within oil droplets. Reverse transcription was then performed within the oil droplets to produce barcoded cDNA. Barcoded cDNA was isolated from the partitioning oil using Silane DynaBeads (Thermo Fisher Scientific, Dynabeads MyONE Silane, Cat# 37002D) before amplification by PCR. Cleanup/size selection was performed on amplified cDNA using SPRIselect beads (Beckman-Coulter, SPRIselect, Cat# B23317) and cDNA quality was assessed using an Agilent 2100 BioAnalyzer and the High-Sensitivity DNA assay (Agilent, High-Sensitivity DNA Kit, Cat# 5067-4626). Sequencing libraries were generated from cleaned, amplified cDNA using the 10X Chromium Kit’s including reagents for fragmentation, sequencing adaptor ligation, and sample index PCR. Between each of these steps, libraries were cleaned and size selected using SPRIselect beads. Final quality of cDNA libraries was once again assessed using the Agilent BioAnalyzer High-Sensitivity DNA assay, and quality-confirmed libraries were sequenced using Illumina’s NextSeq platform. All reagents are listed in Supplementary Table 2.


Single cell analysis


LPS treatment: RNA-seq processing and downstream analysis: We used CellRanger version 2.1.0 (10X Genomics) to pool and process the raw RNA sequencing data. First, using CellRanger mkfastqc pipeline, each sample sequencing library was demultiplexed based on the sample index read to generate FASTQ files for the paired-end reads. STAR aligner39 was used to align reads to the human reference genome (GRCh38) through the CellRanger count pipeline. After alignment, all sample libraries were equalized to the same sequencing depth (each sample cell is downsampled to have the same confidently mapped reads per cell) and aggregated together subsequently to generate a gene-cell barcode matrix using CellRanger aggr pipeline.


After data aggregation, we performed all filtering, normalization, and scaling of data using Seurat suite version 2.3 (refs. 12,13). Cells with less than 300 and greater than 2000 detected genes were filtered out, as well as cells with greater than 10,000 UMIs and greater than 10% mitchondrial counts were filtered out. Genes that were detected in less than 10 cells were removed. Gene counts for each cell were normalized by total expression, multiplied by a scale factor of 10,000 and transformed to log scale.


PCA based on the highly variable genes detected (dispersion of 2) was performed for dimension reduction and the top 20 principal components (PCs) were selected. We clustered cells based on graph-based methods (KNN and Louvain community detection method) implemented in Seurat. The clusters and other known annotations were visualized using t-stochastic neighbor embedding (t-SNE)40.


Cell hashing processing and analysis: For hashtag oligo (HTO) quantification, we first ran Cite-seq-Count15,41 on the HTO fastq files to process the HTO reads with the parameters specific to 10X Genomics single-cell 3′ v2 data as stated in https://github.com/Hoohm/CITE-seq-Count. In addition, we used CellRanger v.3.0.2 (10X Genomics) to process the raw sequencing RNA reads and Seurat suite for downstream analyses. To identify the cells sample-of-origin, we demultiplexed the HTOs and removed doublets and ambiguous cells using the Seurat pipeline for demultiplexing as mentioned in https://satijalab.org/seurat/hashing_vignette.html.


IFNβ treatment: After demultiplexing the HTOs, we performed all downstream analyses as described above. Cells with less than 200 and greater than 2500 detected genes were filtered out, as well as cells with greater than 5% mitochondrial counts were filtered out.


In vitro morphine treatment: After demultiplexing the HTOs, we performed all downstream analyses as described above. Cells with less than 500 and greater than 3000 detected genes were filtered out, as well as cells with greater than 5% mitochondrial counts were filtered out.


RNA sequencing DE analysis: To identify peripheral immune subpopulations, we performed differential expression analysis using Wilcoxon rank-sum test between clusters to identify top expressing genes for each cluster for cell type identification implemented in Seurat. Cell-type-specific gene signatures were determined from the overlap of more highly expressed and canonical gene markers.


We performed differential expression analysis for each cell type between control cells and opioid-dependent cells using Model-based Analysis of Single Cell Transcriptomics (MAST)42. Utilizing this method, we fit a hurdle model modeling condition and the centered cellular detection rate (cngeneson), and then performed a likelihood ratio test dropping the condition term to identify genes upregulated and downregulated in opioid-dependent samples compared to controls. Differentially expressed genes were evaluated according to their log fold change (greater than log2(1.5)) and adjusted p values (0.05). All figures generated using ggplot2 R package43.


Enrichment analysis: We performed gene enrichment analysis of the list of differential genes between opioid-dependent individuals and non-dependent controls for each cell type using Metascape44 online tool (http://metascape.org/). The enrichment analysis was run using default settings, and was assessed and visualized through a heatmap of significance (−log(p value)). All heatmaps were generated using ComplexHeatmap R package and color scale generated using dependent R package circilize45.


Reporting summary


Further information on research design is available in the Nature Research Reporting Summary linked to this article.