Prenatal methadone exposure selectively alters protein expression in primary motor cortex: Implications for synaptic function

As problematic opioid use has reached epidemic levels over the past 2 decades, the annual prevalence of opioid use disorder (OUD) in pregnant women has also increased 333%. Yet, how opioids affect the developing brain of offspring from mothers experiencing OUD remains understudied and not fully understood. Animal models of prenatal opioid exposure have discovered many deficits in the offspring of prenatal opioid exposed mothers, such as delays in the development of sensorimotor function and long-term locomotive hyperactivity. In attempt to further understand these deficits and link them with protein changes driven by prenatal opioid exposure, we used a mouse model of prenatal methadone exposure (PME) and preformed an unbiased multi-omic analysis across many sensoriomotor brain regions known to interact with opioid exposure. The effects of PME exposure on the primary motor cortex (M1), primary somatosensory cortex (S1), the dorsomedial striatum (DMS), and dorsolateral striatum (DLS) were assessed using quantitative proteomics and phosphoproteomics. PME drove many changes in protein and phosphopeptide abundance across all brain regions sampled. Gene and gene ontology enrichments were used to assess how protein and phosphopeptide changes in each brain region were altered. Our findings showed that M1 was uniquely affected by PME in comparison to other brain regions. PME uniquely drove changes in M1 glutamatergic synapses and synaptic function. Immunohistochemical analysis also identified anatomical differences in M1 for upregulating the density of glutamatergic and downregulating the density of GABAergic synapses due to PME. Lastly, comparisons between M1 and non-M1 multi-omics revealed conserved brain wide changes in phosphopeptides associated with synaptic activity and assembly, but only specific protein changes in synapse activity and assembly were represented in M1. Together, our studies show that lasting changes in synaptic function driven by PME are largely represented by protein and anatomical changes in M1, which may serve as a starting point for future experimental and translational interventions that aim to reverse the adverse effects of PME on offspring.

Introduction

As problematic opioid use has reached epidemic levels over the past 2 decades, the annual prevalence of opioid use disorder (OUD) in pregnant women has also increased 333% (Haight et al., 2018). A product of this increased opioid use has also lead to a seven-fold increase in neonates experiencing opioid withdrawal syndrome (NOWS) following prenatal opioid exposure (POE) (Epstein et al., 2013; Patrick et al., 2015; Ko et al., 2016; Ko et al., 2017; Duffy et al., 2018; Winkelman et al., 2018; Honein et al., 2019). POE negatively impacts neurodevelopment, producing dysfunction of motor control and persistent behavioral hyperactivity, with potential long-term negative effects on quality of life (Sundelin Wahlsten and Sarman 2013; Lee et al., 2019; Yeoh et al., 2019). Recent data also show that POE neonates who display NOWS symptoms at birth have better outcomes than POE neonates that did not display NOWS, likely a function of the treatments used to alleviate NOWS potentially having lasting protective effects (Leyenaar et al., 2021). Despite the proven potential for treatments to improve the lives of neonates with POE, there has been a lack of new treatments as current first-line treatments for both pregnant women with OUD and neonates displaying NOWS remain opioid maintenance therapies, such as methadone and buprenorphine (Committee Opinion, 2017; Ghazanfarpour et al., 2019; Tang et al., 2021; Zankl et al., 2021). As a consequence, it is difficult to separate the pharmacology of POE and NOWS treatments as well as other environmental factors, both independently and in combination with each other across the human lifespan. This complexity raises questions regarding the specific impacts of opioids and NOWS treatments alone. Ultimately, it is unknown exactly how POE impacts the brain and behavioral development of these neonates as they mature (Brogly et al., 2018; Rose-Jacobs et al., 2019; Labella et al., 2021).

Previous rodent studies have revealed some important findings regarding the specific impacts of POE on neurobiological function and/or behavioral development through the lifespan (Grecco and Atwood 2020; Boggess and Risher 2022). However, many of these prior studies have inadequately modelled human clinical cases, making their findings difficult to interpret and hindering translation to human therapeutics. Recently, our lab has developed a translational mouse model that more closely resembles the typical pattern of opioid use in a pregnant woman who is first dependent on oxycodone, then begins methadone maintenance treatment, and subsequently becomes pregnant while maintained on methadone (Grecco et al., 2021b). This prenatal methadone exposure (PME) model produces many outcomes reminiscent of human neonatal opioid withdrawal. In this model, we found that PME reduces physical growth in offspring persisting into adolescence and delays the development of sensorimotor function during the pre-weaning period. Furthermore, PME produced locomotor hyperactivity that persisted into adolescence as well.

The potential neural mechanisms that mediate the sensorimotor deficits and hyperactivity observed in our model are diverse. We previously characterized PME-induced physiological and biochemical changes in the dorsal striatum, a brain region that integrates input from sensorimotor cortices (Grecco et al., 2022b). We have also described the effects of PME on the biochemical, anatomical and electrophysiological properties of the primary somatosensory cortex (S1) (Grecco et al., 2022a). The effects we measured in these brain regions likely contribute to the sensorimotor behavioral changes we observed as a result of PME. However, the primary motor cortex (M1) likely also contributes, given its central role in the preparation, execution, and adaptation of motor movements (Peters et al., 2017; Svoboda and Li 2018). Indeed, we showed that PME alters the physiological parameters of layer five pyramidal neurons within M1 with similar effects in both males and females (Grecco et al., 2021b). Specifically, PME enhanced excitatory glutamatergic input to these neurons from outer M1 layers and enhanced voltage sag, an electrophysiological property related to the function of hyperpolarization-activated cyclic nucleotide-gated (HCN) channels (Sheets et al., 2011). However, in our previous study of M1 we failed to explore the biochemical changes produced by PME in M1 that could underlie the physiological changes it induced. Thus, we sought to profile the protein expression and phosphorylated protein states in this brain region and compare these changes to those in S1 and the lateral (DLS) and medial (DMS) divisions of the dorsal striatum. By tracking the protein expression and protein phosphorylation states altered by PME across distinct brain regions, we can learn if PME has differential effects across different brain regions, and further elucidate how POE uniquely alters infants in hopes of refining our approaches to treating NOWS and mitigate POE-induced long-lasting negative outcomes.

Materials and methods

Animals

Animal care and research were conducted in accordance with guidelines established by the National Institutes of Health and protocols were approved by the Indiana University School of Medicine Institutional Animal Care and Use Committee. Eight-week-old female C57BL/6 J mice were acquired from Jackson Laboratories (Bar Harbor, Maine), single housed, and randomly assigned to either saline (10 mL/kg) or oxycodone treatments. Oxycodone dependence was induced by a dose-ramping procedure with a dose of 10 mg/kg administered on pregestational day (PG) 14, 20 mg/kg on PG13, and then maintained on 30 mg/kg for PG12-6. All saline or oxycodone doses were administered subcutaneously twice daily at least 7 h apart. On PG5, oxycodone-treated mice were transitioned to 10 mg/kg methadone while saline-treated animals continued to receive saline injections (s.c. b. i.d.). All oxycodone and methadone solutions were prepared in saline. Oxycodone and methadone were obtained from the National Institute on Drug Abuse Drug Supply Program. Five days following initiation of methadone treatment, an 8-week-old C57BL/6 J male mouse (also acquired from Jackson Laboratories) was placed into the cage of each female for 4 days. Mucous plugs were assessed each morning to approximate gestational day G) 0.

Cages were examined for the presence of pups at the time of each morning and afternoon injection, and the day of birth was designated postnatal day P) 0. Only litters between three and eight pups were used in subsequent studies of offspring. To maintain a sufficient level of methadone exposure to the newborn pups during P0-P7, the dams were maintained at their highest gestational dose after giving birth and up to P7. After P7, the dose of methadone administered to dams was adjusted to their body weight. All treatments to dams were discontinued at weaning (∼P28). For more detailed methods on methadone exposure and the associated behavioral changes, please see our previous publication that defines the PME model (Grecco et al., 2021b).

Early adolescent offspring (P21–P36) were used for the proteomics and immunohistochemical studies. Sample preparation, mass spectrometry analysis, bioinformatics, and data evaluation for quantitative proteomics and phosphoproteomics experiments were performed in collaboration with the Indiana University Proteomics Core similar to our previous studies (Grecco et al., 2021a; Grecco et al., 2022a). Our previous studies of motor behavior and M1 physiology did not reveal sex differences. In addition, statistical power calculations indicated that the number of samples available to us for our proteomics assessments across all the brain regions sampled in this manuscript would not be sufficient to detect sex by prenatal exposure interactions for every single brain region sampled. Therefore, in this study we limited our comparisons to prenatal treatment alone to be conservative in our findings and consistent across all data types presented, but also included data from both males and females in each dataset. For immunohistochemistry data, we have illustrated data that came from males and females throughout the paper so that the reader might discern whether future, more powered studies may reveal sex differences. All data in this paper is also provided in Supplementary Materials should any reader want to re-analyze these findings for their own hypotheses.

M1 proteomics and phosphoproteomics

Protein preparation

For proteomics and phosphoproteomics, we followed similar methods on which we have previously published, with the exception that we used the more recently updated mouse reference proteome database (Grecco et al., 2021a; Grecco et al., 2022a). Animals were rapidly decapitated without anesthesia by blinded researcher and tissue was dissected. Slices were cut in a 0.5 mm coronal mouse brain matrix and whole M1, S1, DMS, and DLS were dissected from each slice. Tissue was immediately snap frozen in isopentane on dry ice and stored at −80°C. Flash frozen brain lysates were homogenized using a BeadBug™ 6 (Benchmark scientific Cat No: D1036, 3 mm zirconium beads Cat No: D1032-30, 10 rounds of 30 × 30 s, 4°C) in 1 mL of 8 M urea (CHEBI: 16199) in 100 mM Tris, pH 8.5 (CHEBI: 9754). Samples were then sonicated on a Bioruptor^® sonication system (Diagenode Inc. United States, North America Cat No: B01020001) with 30 s/30 s on/off cycles for 15 min in a water bath at 4°C. After subsequent centrifugation at 14,000 rcf for 20 min, protein concentrations were determined by Bradford protein assay (BioRad Cat No: 5000006). 100 µg equivalent of protein from each sample were then treated with 5 mM tris (2-carboxyethyl)phosphine hydrochloride (Sigma-Aldrich Cat No: C4706) to reduce disulfide bonds and the resulting free cysteine thiols were alkylated with 10 mM chloroacetamide (Sigma Aldrich Cat No: C0267). Samples were diluted with 50 mM Tris pH 8.5 (Sigma-Aldrich Cat No: 10812846001) to a final urea concentration of 2 M for overnight Trypsin/Lys-C digestion at 35°C (1:100 protease:substrate ratio, Mass Spectrometry grade, Promega Corporation, Cat No: V5072 (Levasseur et al., 2019; Li et al., 2020).

Peptide purification and labeling

We performed peptide purification and labeling according to previously published methods (Grecco et al., 2021a; Grecco et al., 2022a). To summarize our previous methods, digestion was halted by addition of 0.3% final v/v trifluoroacetic acid (TFA), and peptides were desalted on Waters Sep-Pak^® Vac cartridges (Waters™ Cat No: WAT054955) with a wash of 1 mL 0.1% TFA followed by elution in 0.6 mL of 70% acetonitrile 0.1% formic acid (FA). Peptides were dried by speed vacuum and resuspended 50 mM triethylammonium bicarbonate. Peptide concentrations were checked by Pierce Quantitative colorimetric assay (Cat No: 23275). The same amount of peptide from each sample was then labeled for 2 h at room temperature, with 0.5 mg of Tandem Mass Tag Pro (TMTpro) reagent (16-plex kit, manufactures instructions Thermo Fisher Scientific, TMTpro™ Isobaric Label Reagent Set; Cat No: 44520, lot no. VI310352) (Li et al., 2020). Labeling reactions were quenched with 0.3% hydroxylamine (final v/v) at room temperature for 15 min. Labeled peptides were then mixed and dried by speed vacuum. The TMT-labeled peptide mix was desalted to remove excess label using a 100 mg Waters SepPak cartridge, eluted in 70% acetonitrile, 0.1% formic acid and lyophilized to dryness.

Phosphopeptide enrichment

Phosphopeptides were enriched from the mixed, labeled peptides on one spin tip from a High-Select™ TiO2 Phosphopeptide Enrichment Kit (capacity of 1–3 mg; Thermo Fisher Scientific, catalog A32993). After preparing spin tips, labeled and mixed peptides were repeatedly applied to the TiO2 spin tip, eluted, and immediately dried as per manufacturer’s instructions. Prior to LC/MS/MS the phosphopeptides were resuspended in 25 µL 0.1% formic acid. The flow through from each tip was saved for global proteomics.

High pH basic fractionation

The phosphopeptide enrichment step flow through and washes were dried down and approximately 120 µg peptides were fractionated using the TMT fractionation protocol of Pierce high pH basic reversed-phase peptide fractionation kit (Thermo Fisher Scientific™ Cat no 84858; with a wash of 5% acetonitrile, 0.1% triethylamine (TEA) followed by elution in 10%, 12.5%, 15%, 17.5%, 20%, 22.5%, 25%, 30%, and 70% acetonitrile, all with 0.1% TEA). Fractions were dried down and resuspended in 50 µL 0.1% FA prior to online LC-MS.

Nano-LC-MS/MS analysis

We performed nano-LC-MS/MS analysis with previously published methods (Grecco et al., 2021a; Grecco et al., 2022a). To summarize our previous methods, nano-LC-MS/MS analyses were performed on an EASY-nLC™ HPLC system (SCR: 014993, Thermo Fisher Scientific) coupled to Orbitrap Fusion™ Lumos™ mass spectrometer (Thermo Fisher Scientific). One-fifth of the phosphopeptides and one-10th of each global peptide fraction was loaded onto a reversed phase EasySpray™ C18 column (2 μm, 100 Å, 75 μm × 25 cm, Thermo Scientific Cat No: ES902A) at 400 nL/min. Peptides were eluted from 4% to 28% with mobile phase B [Mobile phases A: 0.1% FA, water; B: 0.1% FA, 80% Acetonitrile (Fisher Scientific Cat No: LS122500)] over 160 min; 28%–35% B over 5 min; 35%–50% B for 14 min; and dropping from 50% to 10% B over the final 1 min. The mass spectrometer method was operated in positive ion mode with a 4 s cycle time data-dependent acquisition method with advanced peak determination and Easy-IC (internal calibrant). Precursor scans (m/z 400–1750) were done with an orbitrap resolution of 120,000, RF lens% 30, maximum inject time 50 m, standard AGC target, including charges of two to six for fragmentation with 60 s dynamic exclusion. MS2 scans were performed with a fixed first mass of 100 m/z, 34% fixed CE, 50000 resolution, 200% normalized AGC target and dynamic maximum IT. The data were recorded using Thermo Fisher Scientific Xcalibur (4.3) software (Thermo Fisher Scientific Inc.)

Proteome and phosphoproteome data processing

We performed M1, S1, DLS, and DMS data processing similar to our previously published methods with the exception that we analyzed our resulting M1 and previous RAW files in Proteome Discover™ 2.5 (Thermo Fisher Scientific, RRID: SCR_014477) with a more recent Mus musculus UniProt FASTA (both reviewed and unreviewed sequences, downloaded 05/13/2022, 84,754 sequences plus common contaminants (71 sequences)) (Grecco et al., 2021a; Grecco et al., 2022a). SEQUEST HT searches were conducted with a maximum number of three missed cleavages; precursor mass tolerance of 10 ppm; and a fragment mass tolerance of 0.02 Da. Static modifications used for the search were, 1) carbamidomethylation on cysteine C) residues; 2) TMTpro label on lysine K) residues and the N-termini of peptides. Dynamic modifications used for the search were TMTpro label on N-termini of peptides, oxidation of methionines, phosphorylation on serine, threonine or tyrosine, and acetylation, methionine loss or acetylation with methionine loss on protein N-termini. Percolator False Discovery Rate was set to a strict setting of 0.01 and a relaxed setting of 0.05. IMP-ptm-RS node was used for all modification site localization scores. Values from both unique and razor peptides were used for quantification. In the consensus workflows, peptides were normalized by total peptide amount with no scaling. Quantification methods utilized TMTpro isotopic impurity levels available from Thermo Fisher Scientific. Reporter ion quantification was allowed with S/N threshold of seven and co-isolation threshold of 50%. Data shown is for PME/PSE abundance value ratios (AR). Resulting grouped abundance values for each sample type, AR values; and respective p-values (t-test) from Proteome Discover™ were exported to Microsoft Excel. Full datasets are provided in the Supplementary Material. The raw data files can be found at https://github.com/dlhagger/PME-MultiRegion-Proteomics.

Immunohistochemistry

We performed immunohistochemistry and its analysis similar to previously published methods (Grecco et al., 2021a; Grecco et al., 2022a). To summarize our previous methods, animals were anesthetized with isoflurane and perfused with 4% paraformaldehyde prepared in PBS for 10 min at a pump rate of ∼2 mL/min. Fixed brains were sectioned into 100 μm sections in the coronal plane (between bregma −0.1 and −1.94 mm) using a Leica VT-1000 vibrating microtome (Leica Microsystems) and stored in antigen preserved solution (PBS, 50% ethylene glycol and 1% polyvinyl pyrrolidone) at −20°C. For synaptic marker (VGAT, Gephyrin, VGluT1, and PSD95) staining, sections were permeabilized with 2% Triton X100, then incubated with a blocking solution (3% normal goat serum prepared in PBS with 0.3% Triton X-100) and then incubated overnight with primary antibody prepared in blocking solution (See Table 1 for concentration and source). An appropriate secondary antibody conjugated with an Alexa series fluorophore was used to detect the primary antibody. DAPI (100 ng/mL, Thermo Fisher) was included in the secondary antibody solution to stain nuclei.

TABLE 1

TABLE 1. Antibody descriptions.

For imaging synaptic marker staining we followed similar methods which we have previously published (Grecco et al., 2022a). Z-stack confocal images were acquired from both hemispheres with a Nikon A1 confocal microscope with a 60X/NA1.4 objective at 3 times software zoom or Leica SP8 confocal microscope with a 63X/NA1.2 objective at 2.5 times software zoom. The Z-stacks were taken at 0.1 µm intervals (for VGAT + Gephyrin) or 0.2 µm (for VGluT1 + PSD95), and 2–4 µm-total thickness was imaged. Two images from each hemisphere and both hemispheres were imaged per animal. We utilized Imaris (Bitplane, Zurich, Switzerland) to quantify synaptic punctate at the three-dimensional level and establish the data analysis workflow to quantify the synaptic number according to the published literature (Fogarty et al., 2013; Kuljis et al., 2019; Simhal et al., 2019; Zhou et al., 2019). The volume occupied by nuclei and vasculature varied within each image, thus robustly impacting the density of synaptic marker quantification. To accurately estimate the neuropil occupied volume, we first used surface module to create the surface objects of nuclei and vasculature-like structure. Next, the gephyrin- or PSD95-channel was further masked by nuclei and vasculature objects to exclude the volume occupied by nuclei and vasculature. The post-masked gephyrin or PSD95 channel was used to generate a surface object containing the volume (neuropil object) to be analyzed. For spot detection, we followed similar procedures and parameter settings as described before (Fogarty et al., 2013; Kuljis et al., 2019; Simhal et al., 2019; Zhou et al., 2019). Specifically, the pre-synaptic (VGluT1 and VGAT) and postsynaptic (gephyrin and PSD95) punctate were detected by Imaris spot module with 0.5 µm and 0.3 µm diameter according to the published literature. In general, the diameter for synaptic puncta is between 0.25 and 0.8 µm (Kim and Sheng 2009; Dumitriu et al., 2012). To find the optimal detecting threshold for spot detection, we first manually defined the detecting threshold for one image from each animal. The threshold that detected most synaptic punctates without creating artifacts was applied to analyze all images and generated the spot layer for each synaptic marker. Only synaptic punctates inside the neuropil-object were used for subsequence analysis. Next, we determine how many pre-synaptic spots were directly opposed to postsynaptic spots (defined as synapse at anatomical level) with the distance 0.5 µm. The juxtaposed synaptic punctate of VGluT1/PSD95 and VGAT/gephyrin were defined as intracortical excitatory and inhibitory neurochemical inputs. Synaptic density was calculated as the number of synapses detected in a dataset over the volume of the dataset. All image acquisition and data analysis were performed in a blinded manner.

Kinase-substrate enrichment analysis

A kinase-substrate enrichment analysis (KSEA) of the phosphoproteomics data was performed using the KSEA App (https://casecpb.shinyapps.io/ksea/) (Wiredja et al., 2017). All identified phosphopeptides with quantified abundance ratios (PME/PSE) and confirmed phosphosite modifications were utilized for the KSEA. PhosphoSitePlus + NetworKIN (NetworKIN score cutoff of 0 and substrate count cutoff of 0) were used as the kinase-substrate dataset. Results were FDR-corrected (<0.05), and a z-score of enrichment was calculated to determine the normalized magnitude of up- or downregulation of kinases (PME vs. PSE). The kinase scores resulting from the KSEA analysis were exported to Coral and overlayed onto a kinome tree to visualize patterns in kinase regulation where branch color was set to represent the range of FDR significance values, node color was set to represent the range of enrichment z-scores, and the node size was set to represent the absolute value of enrichment z-scores (Metz et al., 2018). The full results of the GO analysis are provided in the Supplementary Material and at https://github.com/dlhagger/PME-MultiRegion-Proteomics.

Gene ontology enrichment analysis

All analyses are presented as PME relative to PSE (log2 abundance ratios of PME/PSE). For protein expression dendrograms, all differentially abundant (p < 0.05) proteins and phosphopeptides were normalized across animal by their associated UniProt Accession numbers and then hierarchically clustered using the Voorhees algorithm. Predicted protein-protein associative networks and functional enrichments were computed by loading the associated GeneIDs and log2 abundance ratios of differentially abundant proteins and phosphopeptides into STRING (Szklarczyk et al., 2021). Protein-protein and phosphopeptide-phosphopeptide interactions were visualized using the highest confidence interaction score (0.9) with disconnected nodes hidden from the visualization. The interaction network was then kmeans clustered (n = 5) and assigned a pseudo-color to represent each subcluster that was post hoc defined by the top one or two strongest KEGG pathway terms returned by each subcluster’s identity. For overrepresentation analysis of Gene Ontology (GO), the STRING analysis output was filtered to return only significant (p < 0.05) terms for the following classifications: GO Process, GO Function, KEGG, and Reactome. The similarity of GeneIDs and Gene Ontology terms for M1 compared to all other brain regions were computed using cosine similarity between the lists of GeneIDs or GO terms to demonstrate how closely related these identifiers were between brain regions. The full results of the GO analysis are provided in the Supplementary Material and at https://github.com/dlhagger/PME-MultiRegion-Proteomics.

Statistics and data presentation

Data are graphically presented as dot plots displaying all individual data points with the inclusion of box plots that display interquartile ranges as whiskers. Individual data points that fall beyond the interquartile range are noted with a diamond next to the individual data point. The level of significance was a priori set at p < 0.05. All experiments were performed using both male and female offspring. To minimize potential litter effects in all completed studies, no more than two males and females per litter were utilized for any study. Immunostaining studies were not sufficiently powered to detect sex differences with sex considered as a factor (determined by power analysis), and therefore the data presented here are collapsed on sex as a factor. Immunostaining statistical analyses were conducted using pingouin (Vallat 2018). T-tests, with Welch’s correction where appropriate, were used for analyzing all immunostaining data.

Results

Differential protein and phosphopeptide expression

To determine the possible unique alterations that underlie dysfunctional development in PME offspring that occur specifically in the motor cortex, we collected brain punches of M1 cortices from adolescent male and female PME and PSE offspring for quantitative proteomic and phosphoproteomic analyses. In M1 alone, we identified 5,974 proteins and 5,794 phosphopeptides in offspring that were altered by PME. For M1, 1,248 proteins and 391 phosphopeptides were differentially abundant (Figures 1A,B). For the M1 proteome, a substantial downregulation of proteins was observed, with 200 proteins being upregulated and 1,048 proteins being downregulated. Whereas for the M1 phosphoproteome, 199 phosphopeptides were upregulated, and 192 phosphopeptides were downregulated. Clustermaps of the differentially abundant proteins show that PME and PSE animals are most dissimilar, yet male and female animals are most dissimilar for differentially abundant phosphopeptides (Figures 1C,D). For a full detailed list of differentially abundant proteins, phosphopeptides, and clustermaps by individual sample, please see Supplementary Figure S1 and the Supplmentary data sheets 1, 2. Together, these data suggest that PME alters more global protein abundances than phosphopeptides within the motor cortex.

FIGURE 1

FIGURE 1. Differential protein and phosphopeptide expression in the primary motor cortex (M1) of prenatal methadone exposed offspring. Volcano plots for the differential proteome (A) and phosphoproteome (B) with blue circles representing individual proteins/phosphopeptides decreased in PME vs. PSE and red circles representing individual proteins/phosphopeptides increased in PME vs. PSE which reach the level of significance. Clusterplots of differentially abundant proteins (C) and phosphopeptides (D) with example GeneIDs (n = 8 (4M:4 F) PME, 8 PSE (4M:4 F)).

Protein-protein and phosphopeptide-phosphopeptide interaction networks

To further understand the differences in the proteome and phosphoproteome of the primary motor cortex, differentially abundant proteins and phosphopeptide interaction networks were computed. By using the highest confidence interactions and showing only connected nodes, the protein interaction network was constructed then k-means clustered (k = 5) to uncover the identifies of subclusters within the network (Figure 2A). The many M1 proteomic interactions show major KEGG defined clusters for alterations in autophagy, MAPK/Rap1/Ras signaling, metabolic pathways alterations, and changes in glutamatergic signaling at synapses. For M1 phosphoproteomics, a substantially smaller interaction network did not return significant KEGG defined subclusters, therefore a smaller k-means clustering (k = 3) of Reactome terms was used to uncover the identities of network subclusters (Figure 2B). Altered mRNA splicing, membrane trafficking, and presynaptic depolarization and calcium channel opening were identified as the strongest clusters for the M1 phosphoproteome. Interestingly, both networks contain strong alterations in synapses, as well as calcium channel alterations associated with neuronal activity.

FIGURE 2

FIGURE 2. Differential protein and phosphopeptide interaction networks in the primary motor cortex (M1) of prenatal methadone exposed offspring. Interaction networks for the differential proteome (A) and phosphoproteome (B). Interaction networks were k-means clustered and the strongest KEGG or Reactome term of the subcluster was used to define the individual subcluster Solid lines represent interactions between proteins or phosphopeptides in the same cluster and dashed lines represent interactions between proteins or phosphopeptides across clusters.

Kinase-substrate enrichment

To further characterize the changes induced in the M1 phosphoproteome, we performed a kinase-substrate enrichment analysis (KSEA) to estimate the changes in kinase pathways based on specific phosphorylation site modifications for the sampled M1 phosphoproteome (Wiredja et al., 2017). The full list of KSEA kinase z-scores and the significant FDR corrected kinase-associated gene IDs were returned (Figure 3A), showing that many of the highly significant altered kinases are known to be calcium dependent and associated with regulating glutamate receptors (Mao et al., 2014). The KSEA output was also overlayed onto a kinome tree using Coral (Metz et al., 2018) to visualize the patterns of enrichment from significantly enriched kinases (Figure 3B). While the CAMK family has some of the strongest enrichments as measured by FDR corrected significance, the STE, CK1, and TKL families of kinases, which are similarly related, in addition to the CGMC and the AGC kinase families have multiple strong z-score alterations in their functions due to PME exposure in M1.

FIGURE 3

FIGURE 3. Dysregulated kinases and kinome mapping in the primary motor cortex (M1). The results of a kinase-substrate enrichment analysis (A) showing kinases that are significantly dysregulated between PME and PSE (blue bars represent decreased kinases and red bars represent increased kinases as measured by z-score) with a table listing all significantly FDR corrected kinases. The results of this kinase-substrate enrichment were mapped onto a kinome tree (B) via Coral to further visualize changes in kinases due to PME. Branch color was set to represent the range of FDR significance values, node color was set to represent the range of enrichment z-scores, and the node size was set to represent the absolute value of enrichment z-scores.

Gene and gene ontology functional enrichments

Although the number of individual proteins and the abundance change magnitudes were larger in the global proteome than the enriched phosphoproteome PME in M1, we sought to further investigate the overlap between PME-induced changes in the proteome and phosphoproteome. By looking at the 89 overlapping associated GeneIDs from differentially abundant proteins and phosphopeptides and returning the five greatest positive and negative changes in log₂ (abundance ratios), there exists an overlap in targets such as Grin2a, Dlgap3, Ank2, Nrgn, and Mobp (Figure 4A). These genes are all known to regulate synaptic organization, activity, and are known targets in disorders associated with aberrant neuronal activity or synapse loss (Bulat et al., 2014; Wan et al., 2014; Hegyi 2017; Rasmussen et al., 2017; Amador et al., 2020; Pereira et al., 2021). Similarly, we looked at the top five positive and negative changes in log₂ (abundance ratios) for proteins that had no accompanied changes in the phosphoproteome (Figure 4A). Interestingly, the most strongly upregulated protein Adra2, the Alpha-2-adrenergic receptor, is a target of a currently used treatment for neonatal opioid exposure, clonidine (Ghazanfarpour et al., 2019). We also examined the top five positive and negative changes in log₂ (abundance ratios) for phosphopeptides that had no accompanied changes in the proteome (Figure 4A). Targets such as Cep170b, Cacan1e, Camk2b, Clasp2, Ppfia4, and Smarca4 are all known to alter synaptic assemblies and plasticity, are associated with calcium channelopathies and modulation of acetylcholine and other receptors, and are associated with disease states defined by aberrant neuronal firing, such as epilepsy (Basu et al., 2014; Nicole and Pacary 2020; Kessi et al., 2021; Li et al., 2021).

FIGURE 4

FIGURE 4. Overlap in the proteome, phosphoproteome, and gene ontology enrichments between prenatal methadone exposed primary motor cortex (M1) (A) Of 1426 differentially abundant proteins and phosphopeptides matched to their corresponding GeneIDs, 89 were both similarly altered in the proteome and phosphoproteome with the top five largest positive and negative log₂ (abundance ratio) listed. Overlap between M1 proteins and phosphopeptides for gene ontology enrichments for GO Process (B), GO Function (C), KEGG (D), and Reactome (E), listed with the top five terms for each enrichment, sorted by the strength of enrichment.

We similarly analyzed overlap in PME-altered expression of M1 proteins and phosphopeptides using gene ontology databases to further characterize the functional changes beyond specific associated gene targets using STRING (Szklarczyk et al., 2021). The overlap between GO Processes showed 279 matching terms (Figure 4B). When ranking these GO Processes by their greatest strength, the five greatest terms all are associated with synapse maintenance, activity, and structure. Similarly, GO Function overlap in M1 shows 58 terms that are shared with the top five strongest terms all relating to structural changes at synapses and glutamate receptor activity (Figure 4C). 31 overlapping KEGG terms show that three of the top five strongest terms are enriched for synaptic plasticity and synaptic events, (Figure 4D), and 31 overlapping Reactome terms show that all of the strongest terms relate to synaptic changes, especially at glutamatergic synapses with NMDA receptors (Figure 4E). Together, the overlap between gene and gene ontology enrichments between protein and phosphopeptides overwhelmingly point to PME induced alterations in synaptic assemblies and activity, especially related to glutamatergic synapses.

Assessment of glutamatergic synaptic markers

Given the large number of differentially abundant proteins and phosphopeptides associated with synaptic functioning and the enrichment in associated GeneIDs and terms associated with the synapse, we investigated both glutamatergic and GABAergic synapse density in layer 2/3 (L2/3) and layer 5 (L5) of M1 using the co-localization of presynaptic vesicular glutamate transporter (VGlut1) and postsynaptic protein PSD95 and co-localization of presynaptic vesicular GABA transporter (VGAT) and postsynaptic protein gephyrin, respectively. The co-localization of these makers was used to infer the presence of a “functional” glutamate or GABA synapse from anatomical data.

For glutamate synapses in L2/3 of M1, the density of VGlut1, a protein conventionally considered as marker of intracortical inputs, was unchanged by PME (Figure 5A). Yet, both PSD95 density (t-test, t₍₅₅₎ = 2.5351, p = 0.0141) and PSD95-VGlut1 co-localization (t-test, t₍₄₉₎ = 4.3503, p = 0.0001) were increased due to PME exposure (Figures 5B,C). A representative confocal stack of PSD95 and VGlut1 double-stained image between PME and PSE visualizes these differences (Figure 5D). Similarly, for glutamate synapses in L5 of M1, the density of VGlut1 was unchanged by PME (Figure 5E). Both, PSD95 density (t-test, t₍₅₈₎ = 3.0883, p = 0.0031) and PSD95-VGlut1 co-localization (t-test, t₍₅₈₎ = 4.0594, p = 0.0001) were also increased due to PME exposure (Figures 5F,G). For both L2/3 and L5, PME substantially increases the density of functional glutamatergic synapses in M1.

FIGURE 5

FIGURE 5. Neurochemical assessment of functional Glutamatergic and GABAergic synapses. In Layer 2/3 (A) vesicular glutamate transporter (VGlut1) density was unaltered, but PME increases (B) post-synaptic density (PSD95) density (t-test, p = 0.0141) and (C) increases PSD95-VGlut1 colocalization (t-test, p = 0.0001). An exemplary (D) confocal stack of PSD95 (green) and VGlut1 (red) double-stained image between PME and PSE, noting the raw image, the output of spot detection using Imaris to identify puncta, and colocalized spots defined by spots within 0.5 µm of each other. In Layer 5 (E) VGlut1 density was unaltered, but PME increases (F) PSD95 density (t-test, p = 0.0031) and (G) increases PSD95-VGlut1 colocalization (t-test, p = 0.0001). n = 8 (3M:5 F) PME, 8 PSE (4M:4 F). In Layer 2/3 both (H) vesicular GABA transporter (VGAT) density and (I) post-synaptic density (gephyrin) density was unaltered, but PME (J) decreases Gephyrin-VGAT colocalization (t-test, p = 0.0005). An exemplary (K) confocal stack of gephyrin (green) and VGAT (red) double-stained image between PME and PSE, noting the raw image, the output of spot detection using Imaris to identify puncta, and colocalized spots defined by spots within 0.5 µm. In Layer 5 (L) VGAT density was unaltered, but PME increases (M) Gephyrin density (t-test, p = 0.0133) and (N) decreases gephyrin-VGAT colocalization (t-test, p = 0.0005). n = 9 (4M:5 F) PME, 8 PSE (4M:4 F).

For GABA synapses in L2/3 of M1, the density of VGAT and gephyrin were both unchanged by PME (Figures 5H,I). Yet, VGAT-gephyrin co-localization (t-test, t₍₅₈₎ = −3.7058, p = 0.0005) was decreased due to PME exposure (Figure 5J). A representative confocal stack of gephyrin and VGAT double-stained image between PME and PSE visualizes these changes (Figure 5K). In L5, VGAT was also unaltered (Figure 5L). Gephyrin density was increased (t-test, t₍₆₀₎ = 2.5499, p = 0.0133) but, gephyrin-VAT co-localization (t-test, t₍₅₈₎ = −3.6675, p = 0.0005) was decreased due to PME exposure (Figures 5M, N). For L2/3 and L5, PME decreased the density of functional GABAergic synapses. Together, the net effect of both glutamate and GABA synapses suggests a net strengthening of glutamatergic synapses in M1, across multiple cortical layers.

Multi-region differential protein and phosphopeptide expression

The combination of protein, phosphopeptide, and anatomical alterations in synapse structure and function in M1 provoked us to ask if there were other brain regions that displayed similar changes due to PME. We updated our analyses of prior proteomic and phosphoproteomic datasets using searches of more up-to-date databases to determine if S1, DLS, and DMS displayed similar changes in protein and phosphopeptide alterations (Grecco et al., 2022a; Grecco et al., 2022b). We preformed the identical quantitative proteomic and phosphoproteomic analyses that we did in M1. In S1, we identified 10,818 proteins and 2,669 phosphopeptides in our samples, but only 1,637 proteins and 263 phosphopeptides were differentially abundant (Figures 6A,B). Similar to the M1 proteome and phosphopeptides, S1 showed a large number of significantly downregulated proteins compared to upregulated proteins, and an even balance of upregulated and downregulated phosphopeptides, which may be most strongly driven by PME males (Figures 6C,D). In DLS, we identified 8,588 proteins and 4,606 phosphopeptides in our samples, but only 839 proteins and 381 phosphopeptides were differentially abundant (Figures 6E,F). The DLS proteome and phosphopeptides were similarly up- and downregulated with the largest effects driven by PME, not sex (Figures 6G,H). Finally, in DMS we identified 9,411proteins and 4,594 phosphopeptides in our samples, but only 1,304 proteins and 282 phosphopeptides were differentially abundant (Figures 6I,J). The DMS proteome, like M1, saw many proteins downregulated (∼80%) and an even distribution in phosphopeptide changes with PME males being most dissimilar from PSE males for the proteome (Figure 6K,L). Altogether, M1 showed similar patterns of changes with S1 and DMS in that there was a bias towards decreased protein expression, but a more balanced change in phosphopeptides, while DLS showed a more even distribution of both proteins and phosphopeptides. In regard to the total number of individual proteins that had their abundance altered by PME, S1 showed the largest PME-induced effect. DMS and M1 were nearly equivalently affected and DLS had the least changes due to PME.

FIGURE 6

FIGURE 6. Differential protein and phosphopeptide expression in the primary somatosensory cortex (S1), dorsolateral striatum (DLS), and dorsomedial striatum (DMS) of prenatal methadone exposed offspring. Volcano plots for the differential S1 proteome (A) and phosphoproteome (B) with blue circles representing individual proteins/phosphopeptides decreased in PME vs. PSE and red circles representing individual proteins/phosphopeptides increased in PME vs. PSE which reach the level of significance. Clusterplots of differentially abundant S1 proteins (C) and phosphopeptides (D) with example GeneIDs. Volcano plots for the differential DLS proteome (E) and phosphoproteome (F) with blue circles representing individual proteins/phosphopeptides decreased in PME vs. PSE and red circles representing individual proteins/phosphopeptides increased in PME vs. PSE which reach the level of significance. Clusterplots of differentially abundant DLS proteins (G) and phosphopeptides (H) with example GeneIDs. Volcano plots for the differential DMS proteome (I) and phosphoproteome (J) with blue circles representing individual proteins/phosphopeptides decreased in PME vs. PSE and red circles representing individual proteins/phosphopeptides increased in PME vs. PSE which reach the level of significance. Clusterplots of differentially abundant DMS proteins (K) and phosphopeptides (L) with example GeneIDs (n = 8 (4M:4 F) PME, 8 PSE (4M:4 F)).

Multi-region gene and gene ontology functional enrichments

We also sought to understand the overlap between the PME-altered proteome and phosphoproteome in other brain regions. Comparing the M1 overlap, which implicated significant changes in synapse structure and activity, to the overlap in these other brain regions would allow us to understand the uniqueness of PME-induced protein and phosphopeptide changes in each region in the context of altered associated genes, and alterations in biological processes. Similarly, we performed gene and gene ontology analyses on S1, DLS, and DMS regions for both the proteome and phosphoproteome. For the proteome of all brain regions, only one associated gene target was shared (Figure 7A), Pdxk which encodes pyridoxal kinase that converts inactive B6 vitamers into the active cofactor (Keller et al., 2020). Further, the overlap of shared GeneIDs between the proteomes M1-SI (2.76%), M1-DMS (2.51%), M1-DLS (3.47%) were extremely low, suggesting that PME alters the proteome of different brain regions in unique ways (Figure 7A). For the phosphoproteome of all brain regions there was a larger overlap of 19 associated GeneIDs (Figure 7B). The GeneIDs of the top five most positively and negatively altered log₂ (abundance ratios) in the phosphoproteome returned items such as Ank2 and Cldn11, which were also shown in the M1 overlap and found to be altered globally (Figure 7B). Also, there were additional associated GeneIDs not found the in M1 overlap that were found in the global phosphoproteome overlap, such as a Add2, Map2, Pclo, Rims1, Bsn, Cacnb4, and Camk2b which are known to alter synapse function and maintenance (Caddick et al., 1999; Fenster and Garner 2002; Mittelstaedt et al., 2010; Lisman et al., 2012; Waites et al., 2013; Kandimalla et al., 2018). These data suggest that PME similarly alters the phosphoproteome globally, but the lack of these GeneIDs being returned in the global proteome suggests that M1, and other brain regions, have unique susceptibility to PME.

FIGURE 7

FIGURE 7. Overlap in the protein, phosphopeptides, and gene ontology enrichments between prenatal methadone exposed primary motor cortex (M1), primary somatosensory cortex (S1), dorsolateral striatum (DLS), and dorsomedial striatum (DMS). Across the global proteome (A), only one associated GeneID was conserved across all brain regions, and across the global phosphoproteome (B), 19 associated GeneIDs were altered, with the top five strongest positive and negative log2 (abundance ratio) GeneIDs listed. For the GO Process gene ontology enrichments, 9.96% of proteins (C) were similar across all brain regions and M1 to non M1-region terms were assessed using cosine similarity and were 35.14% similar. For the phosphoproteome (D), 25.18% of terms were similar and the M1 to non-M1 similarity was 61.98%. For GO Function enrichments, 8.79% of protein terms (E) were similar with a M1 to non-M1 similarity score of 34.71%, and 18.40% of phosphopeptide terms (F) were similar with a M1 to non-Ma score of 49.21%. For KEGG enrichments, there were not any significantly enriched proteins (G) in the S1 region, so the similarity between protein terms was 10.08% and the M1 to non-M1 similarity score was 32.07%. The KEGG enrichments for the phosphoproteome (H) had no overlapping terms and therefore no M1 to non-M1 similarity score. For the Reactome enrichments for the proteome (I) there were 3.17% similarity between protein terms, and a M1 to non-M1 similarity score of 23.46%. For the phosphoprotein Reactome (J) there was 4.08% overlap and a M1 to non-M1 similarity score of 23.25%.

We also looked at gene ontology term overlaps between all brain regions and compared M1 terms to the overlap of all other terms in attempt to quantitatively score how similar changes in M1 were to other regions of interest. For the global proteome, GO Processes were 9.96% similar, with the top five strongest terms focusing on synaptic transmission, protein localization, and neuronal assemblies (Figure 7C). The GO Process terms for M1 were also compared to the set of all non-M1 brain region terms. This list of terms was vectorized and the cosine similarity between the two vectors were computed, which resulted in a similarity of 0.3415, which can be conceptualized as 34.15% similar for GO Processes for the M1 proteome compared to all other brain region proteomes (Figure 7C). For the global phosphoproteome, GO Processes were 25.18% similar, with the strongest relevant terms involving synapse structure and activity (Figure 7D). M1 to non-M1 similarity was also very high, 61.98%, which mirrored the differences seen between the global proteome and global phosphoproteome for associated GeneIDs (Figure 7D). For GO Function terms in the global proteome, there was an overlap of 8.79% with the strongest terms generally referring to protein binding, with a M1 to non-M1 similarity of 34.71% (Figure 7E). In the global phosphoproteome, GO Function was 18.4% similar with terms strongly associated with synapse structure, and glutamate receptor binding, and M1 to non-M1 similarity was 49.21% (Figure 7F). For KEGG terms in the global proteome, there was no significant S1 KEGG enrichments to report, so the global overlap contains M1, DLS, and DMS only. There was a 10.08% overlap between these terms, with axon guidance being the relevant strongly associated term (Figure 7G). For the KEGG global phosphoproteome, there was no overlap between terms, therefore there was also no M1 to non-M1 similarity score, which can be interpreted as zero (Figure 7H). Finally, the global protein Reactome had 3.17% similarity with protein metabolism and membrane tracking being the most relevant strongly enriched terms, with a M1 to non-M1 similarity score of 23.46% (Figure 7I). The global phosphoproteome had 4.08% similarity with synaptic transmission and the neuronal system being significantly enriched, with a M1 to non-M1 similarity score of 23.25% (Figure 7J). For a full detailed list of differentially abundant proteins and phosphopeptides for each brain region, please Supplmentary data sheets 3–8 for S1, DLS, and DMS regions.

Discussion

Together these data suggest that PME-induced changes to the M1 proteome are uniquely different from the S1, DLS, and DMS proteomes. The average cosine similarity between all gene enrichments is 31.35%, and while some significantly enriched terms generally point to differences in neuronal excitability and protein tracking, only M1 terms return narrow, specific enrichments for glutamate synapses. Also, our protein expression assessments identified changes in specific glutamate receptors, such as ionotropic AMPA and NMDA glutamate receptor subunits (GluA2, GluA3, GluA4, GluN1, GluN2a, and GluN2b) as well as metabotropic glutamate receptors (mGluR1, mGluR2, mGluR3, mGluR5, and mGluR8). Consistent with glutamate signaling proteins being especially affected in M1, we also saw an increase in functional glutamate synapse density in M1, which we did not find in our prior study of S1 (Grecco et al., 2022a). We previously found enhanced glutamate input from layer2/3 M1 neurons to layer five M1 pyramidal neurons (Grecco et al., 2021b). The enhanced glutamate synapse density aligns with this. However, the predominant effect of PME on glutamate receptor and glutamate synapse-related proteins was a decrease in expression. Subcellular localization and further functional analyses are needed to determine the interplay of glutamate protein expression and altered synaptic function. Of particular interest, the most differentially expressed protein in M1 as a result of PME is the alpha2-adrenergic receptor, a target of clonidine, which is an adjunct medication used to treat NOWS (Merhar et al., 2021). This receptor also regulates HCN channels which mediate the voltage sag we identified as being enhanced by PME in M1 pyramidal neurons (Sheets et al., 2011; Grecco et al., 2021b). Relatedly, PME decreased HCN1 channel expression as well as filamin A and Trip8b, two proteins that regulate HCN expression and function (Peters et al., 2022). One of the functions of HCN channels is to filter glutamate transmission, which may also contribute to the altered glutamate transmission we discovered previously (Grecco et al., 2021b). Like our findings of enhanced glutamate transmission, but decreased glutamate receptor expression, the paradox of increased voltage sag, but decreased expression of proteins that mediate that phenomenon will require deeper investigation. Why glutamate signaling-related protein expression in M1 is so profoundly affected by PME relative to other brain regions is a mystery that will require much more study.

We found a decrease in the density of VGAT-gephyrin co-localization in both of the layers of M1 we analyzed, suggesting a decrease in functional GABA synapses. An assessment of the distribution of the data in regard to sex seemed to indicate that this may be driven by males, although future work is needed to more rigorously assess this. Gene ontology and pathway enrichment analyses did not overwhelmingly show that PME uniquely affected GABAergic synaptic proteins like they did for glutamatergic synapses. Nevertheless, we did find some changes in GABAergic synapse-related protein expression. Consistent with the reduction of functional GABAergic synapses, gephyrin expression was decreased by PME as was the vesicular GABA transporter GAT3 and various GABAA and GABAB receptor subunits.

Interestingly, PME effects in the phosphoproteome of all overlapping brain regions and in the M1 phosphoproteome compared to other non-M1 brain regions generalize more broadly. PME similarly affects phosphopeptides across many brain regions by altering the expression of synaptic proteins that mediate neuronal signaling. Further, the average similarity between M1 and non-M1 enrichments is 44.81%. Yet, while greater similarities exist between M1 and other brain regions for the phosphoproteome than do similarities between the M1 proteome and other brain regions, the number (164) and percentage (23.91%) of differentially abundant unique phosphopeptides is greater than all three other brain regions, suggesting that the magnitude in which PME can uniquely alter the phosphoproteome is greatest in M1 than any other brain region sampled (Figure 7B). These data suggest that M1 might be a primary site of PME interaction in the brain, and future studies are needed to delve more mechanistically to understand how this occurs. While the present work identifies specific directionality for up- and downregulation of protein targets as well as alterations in kinase activity and specific phosphopeptide changes, the power and resolution needed to definitively point to mechanistic changes remains a goal of future study. Although, the current analyses may provide starting points for determining which specific protein and phosphopeptides may be susceptible to PME, and how these may serve as a starting point for intervention, both experimentally and translationally.

There are some noteworthy limitations to bear in mind when considering this multi-omics analysis. First, although the differential protein/phosphopeptide expression in PME offspring indicated alterations in synaptic signaling were present, these enrichment analyses do not provide “directionality.” For instance, while the protein/phosphopeptide abundances in PME offspring indicate that “Glutamate Receptor Activity” is enriched, a GO enrichment analysis does not tell us that glutamate activity is increased or decreased, but only that this is significantly represented given the list of differentially abundant proteins/phosphopeptides. Our data regarding decreased glutamate receptor expression, but increased glutamate synapse density and glutamate transmission is a perfect example of the limitations in trying to interpret functional changes from proteomics data. Additionally, one-to-one comparisons between proteomics/phosphoproteomics data with synaptic marker findings remain difficult as bulk M1 tissue was taken for the multi-omics analysis. Therefore, the quantified proteins/phosphopeptides may have originated in various M1 layers, glia cells, interneurons of M1, or even presynaptic inputs from other brain regions. Lastly, we determined that we were not fully powered to detect sex differences in all our proteomic and phosphoproteomic analyses which led to our decision to restrict all our comparisons to prenatal exposure alone. Previous studies, including our own, have discovered sex differences when focusing on specific brain regions. This knowledge motivated our decision to denote female and male sexes where possible in our figures even though we lacked the power to return sex specific effects. It is likely that with increased power, these findings and others may indicate sex by prenatal opioid exposure specific findings, yet future work will need to specifically explore these potential differences and what implications these have for neurobehavioral outcomes of PME in males and females.

In summary, our findings indicate PME induces prominent disruptions in M1 as well as other anatomically and functionally connected brain regions. We also investigated thousands of proteins and phosphopeptides across brain regions that display differential abundance in PME offspring with functional enrichment in several relevant pathways including those related to synaptic transmission and assembly. These findings suggest that deficits in motor behavior function that persist throughout the lifespan due to POE may result from persistent neuroadaptations in M1 and connected brain regions induced by opioid exposure during fetal development. Deciphering how POE produces these changes may lead to novel treatments for mitigating the negative impact of POE on motor and other features of development.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The animal study was reviewed and approved by Animal care and research were conducted in accordance with guidelines established by the National Institutes of Health and protocols were approved by the Indiana University School of Medicine Institutional Animal Care and Use Committee.

Author contributions

GG, J-YH, ED, AM, H-CL, and BA designed experiments. GG, generated animals. ED and AM performed proteomics and phosphoproteomic studies. DH completed proteomic and phosphoproteomic analyses. J-YH completed all immunostaining work. All authors discussed the results and contributed to all stages of manuscript preparation.

Funding

The mass spectrometry work performed in this work was done by the Indiana University Proteomics Core. Acquisition of the IUSM Proteomics core instrumentation used for this project was provided by the Indiana University Precision Health Initiative. The proteomics work was supported, in part, with support from the Indiana Clinical and Translational Sciences Institute funded, in part by Award Number UL1TR002529 from the National Institutes of Health, National Center for Advancing Translational Sciences, Clinical and Translational Sciences Award and the Cancer Center Support Grant for the IU Simon Comprehensive Cancer Center (Award Number P30CA082709) from the National Cancer Institute. This work was also supported by grants awarded by the NIH R01AA027214 (BA), NS086794 (HCL) and F30AA028687 (GG), F31AA029297 (DLH), Indiana University (BA and H-CL), Indiana University Health (BA), and the Stark Neurosciences Research Institute (BA and GG).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2023.1124108/full#supplementary-material

Supplementary Figure S1 | Differential protein and phosphopeptide expression in the primary motor cortex (M1) of prenatal methadone exposed offspring. Clusterplots of differentially abundant proteins (A) and phosphopeptides (B) with example GeneIDs. (n = 8 (4M:4F) PME, 8 PSE (4M:4F)).

References

Amador, A., Bostick, C. D., Olson, H., Peters, J., Camp, C. R., Krizay, D., et al. (2020). Modelling and treating GRIN2A developmental and epileptic encephalopathy in mice. Brain 143 (7), 2039–2057. doi:10.1093/brain/awaa147

PubMed Abstract | CrossRef Full Text | Google Scholar

Basu, S., Sladecek, S., Pemble, H., Wittmann, T., Slotman, J. A., van Cappellen, W., et al. (2014). Acetylcholine receptor (AChR) clustering is regulated both by glycogen synthase kinase 3β (GSK3β)-dependent phosphorylation and the level of CLIP-associated protein 2 (CLASP2) mediating the capture of microtubule plus-ends. J. Biol. Chem. 289 (44), 30857–30867. doi:10.1074/jbc.M114.589457

PubMed Abstract | CrossRef Full Text | Google Scholar

Boggess, T., and Risher, W. C. (2022). Clinical and basic research investigations into the long-term effects of prenatal opioid exposure on brain development. J. Neurosci. Res. 100 (1), 396–409. doi:10.1002/jnr.24642

PubMed Abstract | CrossRef Full Text | Google Scholar

Brogly, S. B., Saia, K. E., Werler, M. M., Regan, E., and Hernandez-Diaz, S. (2018). Prenatal treatment and outcomes of women with opioid use disorder. Obstet. Gynecol. 132 (4), 916–922. doi:10.1097/AOG.0000000000002881

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulat, V., Rast, M., and Pielage, J. (2014). Presynaptic CK2 promotes synapse organization and stability by targeting Ankyrin2. J. Cell. Biol. 204 (1), 77–94. doi:10.1083/jcb.201305134

PubMed Abstract | CrossRef Full Text | Google Scholar

Caddick, S. J., Wang, C., Fletcher, C. F., Jenkins, N. A., Copeland, N. G., and Hosford, D. A. (1999). Excitatory but not inhibitory synaptic transmission is reduced in lethargic (Cacnb4(lh)) and tottering (Cacna1atg) mouse thalami. J. Neurophysiol. 81 (5), 2066–2074. doi:10.1152/jn.1999.81.5.2066

PubMed Abstract | CrossRef Full Text | Google Scholar

Committee Opinion (2017). Committee opinion No. 711: Opioid use and opioid use disorder in pregnancy. Obstet. Gynecol. 130 (2), e81–e94. doi:10.1097/AOG.0000000000002235

PubMed Abstract | CrossRef Full Text | Google Scholar

Duffy, C. R., Wright, J. D., Landau, R., Mourad, M. J., Siddiq, Z., Kern-Goldberger, A. R., et al. (2018). Trends and outcomes associated with using long-acting opioids during delivery hospitalizations. Obstet. Gynecol. 132 (4), 937–947. doi:10.1097/AOG.0000000000002861

PubMed Abstract | CrossRef Full Text | Google Scholar

Dumitriu, D., Berger, S. I., Hamo, C., Hara, Y., Bailey, M., Hamo, A., et al. (2012). Vamping: Stereology-based automated quantification of fluorescent puncta size and density. J. Neurosci. Methods 209 (1), 97–105. doi:10.1016/j.jneumeth.2012.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Epstein, R. A., Bobo, W. V., Martin, P. R., Morrow, J. A., Wang, W., Chandrasekhar, R., et al. (2013). Increasing pregnancy-related use of prescribed opioid analgesics. Ann. Epidemiol. 23 (8), 498–503. doi:10.1016/j.annepidem.2013.05.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenster, S. D., and Garner, C. C. (2002). Gene structure and genetic localization of the PCLO gene encoding the presynaptic active zone protein Piccolo. Int. J. Dev. Neurosci. 20 (3-5), 161–171. doi:10.1016/s0736-5748(02)00046-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fogarty, M. J., Hammond, L. A., Kanjhan, R., Bellingham, M. C., and Noakes, P. G. (2013). A method for the three-dimensional reconstruction of Neurobiotin-filled neurons and the location of their synaptic inputs. Front. Neural Circuits 7, 153. doi:10.3389/fncir.2013.00153

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghazanfarpour, M., Najafi, M. N., Roozbeh, N., Mashhadi, M. E., Keramat-Roudi, A., Megarbane, B., et al. (2019). Therapeutic approaches for neonatal abstinence syndrome: A systematic review of randomized clinical trials. Daru 27 (1), 423–431. doi:10.1007/s40199-019-00266-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Grecco, G. G., and Atwood, B. K. (2020). Prenatal opioid exposure enhances responsiveness to future Drug reward and alters sensitivity to pain: A review of preclinical models and contributing mechanisms. eNeuro 7 (6). doi:10.1523/ENEURO.0393-20.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Grecco, G. G., Haggerty, D. L., Doud, E. H., Fritz, B. M., Yin, F., Hoffman, H., et al. (2021a). A multi-omic analysis of the dorsal striatum in an animal model of divergent genetic risk for alcohol use disorder. J. Neurochem. 157 (4), 1013–1031. doi:10.1111/jnc.15226

PubMed Abstract | CrossRef Full Text | Google Scholar

Grecco, G. G., Mork, B. E., Huang, J. Y., Metzger, C. E., Haggerty, D. L., Reeves, K. C., et al. (2021b). Prenatal methadone exposure disrupts behavioral development and alters motor neuron intrinsic properties and local circuitry. Elife 10, e66230. doi:10.7554/eLife.66230

PubMed Abstract | CrossRef Full Text | Google Scholar

Grecco, G. G., Huang, J. Y., Muñoz, B., Doud, E. H., Hines, C. D., Gao, Y., et al. (2022a). Sex-dependent synaptic remodeling of the somatosensory cortex in mice with prenatal methadone exposure. Adv. Drug Alcohol Res. 2. doi:10.3389/adar.2022.10400

CrossRef Full Text | Google Scholar

Grecco, G. G., Munoz, B., Di Prisco, G. V., Doud, E. H., Fritz, B. M., Maulucci, D., et al. (2022b). Prenatal opioid exposure impairs endocannabinoid and glutamate transmission in the dorsal striatum. eNeuro 9 (2). doi:10.1523/ENEURO.0119-22.2022

CrossRef Full Text | Google Scholar

Haight, S. C., Ko, J. Y., Tong, V. T., Bohm, M. K., and Callaghan, W. M. (2018). Opioid use disorder documented at delivery hospitalization - United States, 1999-2014. MMWR Morb. Mortal. Wkly. Rep. 67 (31), 845–849. doi:10.15585/mmwr.mm6731a1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegyi, H. (2017). Connecting myelin-related and synaptic dysfunction in schizophrenia with SNP-rich gene expression hubs. Sci. Rep. 7, 45494. doi:10.1038/srep45494

PubMed Abstract | CrossRef Full Text | Google Scholar

Honein, M. A., Boyle, C., and Redfield, R. R. (2019). Public Health surveillance of prenatal opioid exposure in mothers and infants. Pediatrics 143 (3), e20183801. doi:10.1542/peds.2018-3801

PubMed Abstract | CrossRef Full Text | Google Scholar

Kandimalla, R., Manczak, M., Yin, X., Wang, R., and Reddy, P. H. (2018). Hippocampal phosphorylated tau induced cognitive decline, dendritic spine loss and mitochondrial abnormalities in a mouse model of Alzheimer's disease. Hum. Mol. Genet. 27 (1), 30–40. doi:10.1093/hmg/ddx381

PubMed Abstract | CrossRef Full Text | Google Scholar

Keller, N., Mendoza-Ferreira, N., Maroofian, R., Chelban, V., Khalil, Y., Mills, P. B., et al. (2020). Hereditary polyneuropathy with optic atrophy due to PDXK variant leading to impaired Vitamin B6 metabolism. Neuromuscul. Disord. 30 (7), 583–589. doi:10.1016/j.nmd.2020.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Kessi, M., Chen, B., Peng, J., Yan, F., Yang, L., and Yin, F. (2021). Calcium channelopathies and intellectual disability: A systematic review. Orphanet J. Rare Dis. 16 (1), 219. doi:10.1186/s13023-021-01850-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, E., and Sheng, M. (2009). The postsynaptic density. Curr. Biol. 19 (17), R723–R724. doi:10.1016/j.cub.2009.07.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Ko, J. Y., Patrick, S. W., Tong, V. T., Patel, R., Lind, J. N., and Barfield, W. D. (2016). Incidence of neonatal abstinence syndrome - 28 states, 1999-2013. MMWR Morb. Mortal. Wkly. Rep. 65 (31), 799–802. doi:10.15585/mmwr.mm6531a2

PubMed Abstract | CrossRef Full Text | Google Scholar

Ko, J. Y., Wolicki, S., Barfield, W. D., Patrick, S. W., Broussard, C. S., Yonkers, K. A., et al. (2017). CDC grand rounds: Public Health strategies to prevent neonatal abstinence syndrome. MMWR Morb. Mortal. Wkly. Rep. 66 (9), 242–245. doi:10.15585/mmwr.mm6609a2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuljis, D. A., Park, E., Telmer, C. A., Lee, J., Ackerman, D. S., Bruchez, M. P., et al. (2019). Fluorescence-based quantitative synapse analysis for cell type-specific connectomics. eNeuro 6 (5). doi:10.1523/ENEURO.0193-19.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Labella, M. H., Eiden, R. D., Tabachnick, A. R., Sellers, T., and Dozier, M. (2021). Infant neurodevelopmental outcomes of prenatal opioid exposure and polysubstance use. Neurotoxicol Teratol. 86, 107000. doi:10.1016/j.ntt.2021.107000

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S. J., Woodward, L. J., and Henderson, J. M. T. (2019). Educational achievement at age 9.5 years of children born to mothers maintained on methadone during pregnancy. PLoS One 14 (10), e0223685. doi:10.1371/journal.pone.0223685

PubMed Abstract | CrossRef Full Text | Google Scholar

Levasseur, E. M., Yamada, K., Pineros, A. R., Wu, W., Syed, F., Orr, K. S., et al. (2019). Hypusine biosynthesis in beta cells links polyamine metabolism to facultative cellular proliferation to maintain glucose homeostasis. Sci. Signal 12 (610), eaax0715. doi:10.1126/scisignal.aax0715

PubMed Abstract | CrossRef Full Text | Google Scholar

Leyenaar, J. K., Schaefer, A. P., Wasserman, J. R., Moen, E. L., O'Malley, A. J., and Goodman, D. C. (2021). Infant mortality associated with prenatal opioid exposure. JAMA Pediatr. 175 (7), 706–714. doi:10.1001/jamapediatrics.2020.6364

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Van Vranken, J. G., Pontano Vaites, L., Schweppe, D. K., Huttlin, E. L., Etienne, C., et al. (2020). TMTpro reagents: A set of isobaric labeling mass tags enables simultaneous proteome-wide measurements across 16 samples. Nat. Methods 17 (4), 399–404. doi:10.1038/s41592-020-0781-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Han, Y., Li, D., Yuan, H., Huang, S., Chen, X., et al. (2021). Identification and validation of a dysregulated miRNA-associated mRNA network in temporal lobe epilepsy. Biomed. Res. Int. 2021, 4118216. doi:10.1155/2021/4118216

PubMed Abstract | CrossRef Full Text | Google Scholar

Lisman, J., Yasuda, R., and Raghavachari, S. (2012). Mechanisms of CaMKII action in long-term potentiation. Nat. Rev. Neurosci. 13 (3), 169–182. doi:10.1038/nrn3192

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, L. M., Jin, D. Z., Xue, B., Chu, X. P., and Wang, J. Q. (2014). Phosphorylation and regulation of glutamate receptors by CaMKII. Sheng Li Xue Bao 66 (3), 365–372.

PubMed Abstract | Google Scholar

Merhar, S. L., Ounpraseuth, S., Devlin, L. A., Poindexter, B. B., Young, L. W., Berkey, S. D., et al. (2021). Phenobarbital and clonidine as secondary medications for neonatal opioid withdrawal syndrome. Pediatrics 147 (3), e2020017830. doi:10.1542/peds.2020-017830

PubMed Abstract | CrossRef Full Text | Google Scholar

Metz, K. S., Deoudes, E. M., Berginski, M. E., Jimenez-Ruiz, I., Aksoy, B. A., Hammerbacher, J., et al. (2018). Coral: Clear and customizable visualization of human kinome data. Cell. Syst. 7 (3), 347–350. doi:10.1016/j.cels.2018.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Mittelstaedt, T., Alvarez-Baron, E., and Schoch, S. (2010). RIM proteins and their role in synapse function. Biol. Chem. 391 (6), 599–606. doi:10.1515/BC.2010.064

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicole, O., and Pacary, E. (2020). CaMKIIβ in neuronal development and plasticity: An emerging candidate in brain diseases. Int. J. Mol. Sci. 21 (19), 7272. doi:10.3390/ijms21197272

PubMed Abstract | CrossRef Full Text | Google Scholar

Patrick, S. W., Davis, M. M., Lehmann, C. U., and Cooper, W. O. (2015). Increasing incidence and geographic distribution of neonatal abstinence syndrome: United States 2009 to 2012. J. Perinatol. 35 (8), 650–655. doi:10.1038/jp.2015.36

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereira, J. B., Janelidze, S., Ossenkoppele, R., Kvartsberg, H., Brinkmalm, A., Mattsson-Carlgren, N., et al. (2021). Untangling the association of amyloid-beta and tau with synaptic and axonal loss in Alzheimer's disease. Brain 144 (1), 310–324. doi:10.1093/brain/awaa395

PubMed Abstract | CrossRef Full Text | Google Scholar

Peters, A. J., Liu, H., and Komiyama, T. (2017). Learning in the rodent motor cortex. Annu. Rev. Neurosci. 40, 77–97. doi:10.1146/annurev-neuro-072116-031407

PubMed Abstract | CrossRef Full Text | Google Scholar

Peters, C. H., Singh, R. K., Bankston, J. R., and Proenza, C. (2022). Regulation of HCN channels by protein interactions. Front. Physiol. 13, 928507. doi:10.3389/fphys.2022.928507

PubMed Abstract | CrossRef Full Text | Google Scholar

Rasmussen, A. H., Rasmussen, H. B., and Silahtaroglu, A. (2017). The DLGAP family: Neuronal expression, function and role in brain disorders. Mol. Brain 10 (1), 43. doi:10.1186/s13041-017-0324-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Rose-Jacobs, R., Trevino-Talbot, M., Vibbert, M., Lloyd-Travaglini, C., and Cabral, H. J. (2019). Pregnant women in treatment for opioid use disorder: Material hardships and psychosocial factors. Addict. Behav. 98, 106030. doi:10.1016/j.addbeh.2019.106030

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheets, P. L., Suter, B. A., Kiritani, T., Chan, C. S., Surmeier, D. J., and Shepherd, G. M. (2011). Corticospinal-specific HCN expression in mouse motor cortex: I(h)-dependent synaptic integration as a candidate microcircuit mechanism involved in motor control. J. Neurophysiol. 106 (5), 2216–2231. doi:10.1152/jn.00232.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Simhal, A. K., Zuo, Y., Perez, M. M., Madison, D. V., Sapiro, G., and Micheva, K. D. (2019). Multifaceted changes in synaptic composition and astrocytic involvement in a mouse model of fragile X syndrome. Sci. Rep. 9 (1), 13855. doi:10.1038/s41598-019-50240-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundelin Wahlsten, V., and Sarman, I. (2013). Neurobehavioural development of preschool-age children born to addicted mothers given opiate maintenance treatment with buprenorphine during pregnancy. Acta Paediatr. 102 (5), 544–549. doi:10.1111/apa.12210

PubMed Abstract | CrossRef Full Text | Google Scholar

Svoboda, K., and Li, N. (2018). Neural mechanisms of movement planning: Motor cortex and beyond. Curr. Opin. Neurobiol. 49, 33–41. doi:10.1016/j.conb.2017.10.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., et al. (2021). The STRING database in 2021: Customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 49 (D1), D605–D612. doi:10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, F., Ng, C. M., Bada, H. S., and Leggas, M. (2021). Clinical pharmacology and dosing regimen optimization of neonatal opioid withdrawal syndrome treatments. Clin. Transl. Sci. 14 (4), 1231–1249. doi:10.1111/cts.12994

PubMed Abstract | CrossRef Full Text | Google Scholar

Vallat, R. (2018). Pingouin: Statistics in Python. J. Open Source Softw. 3 (31), 1026. doi:10.21105/joss.01026

CrossRef Full Text | Google Scholar

Waites, C. L., Leal-Ortiz, S. A., Okerlund, N., Dalke, H., Fejtova, A., Altrock, W. D., et al. (2013). Bassoon and Piccolo maintain synapse integrity by regulating protein ubiquitination and degradation. EMBO J. 32 (7), 954–969. doi:10.1038/emboj.2013.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, Y., Ade, K. K., Caffall, Z., Ilcim Ozlu, M., Eroglu, C., Feng, G., et al. (2014). Circuit-selective striatal synaptic dysfunction in the Sapap3 knockout mouse model of obsessive-compulsive disorder. Biol. Psychiatry 75 (8), 623–630. doi:10.1016/j.biopsych.2013.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Winkelman, T. N. A., Villapiano, N., Kozhimannil, K. B., Davis, M. M., and Patrick, S. W. (2018). Incidence and costs of neonatal abstinence syndrome among infants with medicaid: 2004-2014. Pediatrics 141 (4), e20173520. doi:10.1542/peds.2017-3520

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiredja, D. D., Koyuturk, M., and Chance, M. R. (2017). The KSEA App: A web-based tool for kinase activity inference from quantitative phosphoproteomics. Bioinformatics 33 (21), 3489–3491. doi:10.1093/bioinformatics/btx415

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeoh, S. L., Eastwood, J., Wright, I. M., Morton, R., Melhuish, E., Ward, M., et al. (2019). Cognitive and motor outcomes of children with prenatal opioid exposure: A systematic review and meta-analysis. JAMA Netw. Open 2 (7), e197025. doi:10.1001/jamanetworkopen.2019.7025

PubMed Abstract | CrossRef Full Text | Google Scholar

Zankl, A., Martin, J., Davey, J. G., and Osborn, D. A. (2021). Opioid treatment for opioid withdrawal in newborn infants. Cochrane Database Syst. Rev. 7 (7), CD002059. doi:10.1002/14651858.CD002059.pub4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, B., Chen, L., Liao, P., Huang, L., Chen, Z., Liao, D., et al. (2019). Astroglial dysfunctions drive aberrant synaptogenesis and social behavioral deficits in mice with neonatal exposure to lengthy general anesthesia. PLoS Biol. 17 (8), e3000086. doi:10.1371/journal.pbio.3000086

PubMed Abstract | CrossRef Full Text | Google Scholar