Socially Mediated Shift in Neural Circuits Activation Regulated by Synergistic Neuromodulatory Signaling

Animal maintenance

Adult zebrafish (Danio rerio) AB strain of mixed sex were housed communally in 10-liter tanks (∼20 animals) at the core animal facility and kept at a constant temperature at 28°C under a 14/10 h light/dark cycle. Fish were fed twice daily with high protein commercial food (Otohime B2, Reed Mariculture), and once daily with newly hatched artemia (Brine Shrimp Direct). All experiments were performed in accordance with the Institutional Animal Care and Use Committee.

Social isolation and pairing

Only male zebrafish (6–12 months old) were selected from the group-housed (Communal tanks) and initially isolated physically and visually in individual tanks for one week. Following isolation, fish of equal size and age were paired in new tanks over a two-week period to establish a social dominance relationship. Aggressive interactions (number of attacks and retreats) were observed daily to determine and follow stability of social dominance as described elsewhere (Miller et al., 2017). Unless otherwise noted, only communal male zebrafish were selected from the group-housing tanks to serve as controls for the startle/swim, pharmacological and histologic analysis. Communal males were siblings (same genetic background) and same age as the experimental groups (dominants and submissive fish).

Experimental setup

After the pairing phase was complete, fish were temporarily separated, and behavioral testing was performed sequentially. Each fish was placed in a testing chamber (dimensions: 11 × 4 × 3 cm) filled with double distilled water having a resistance of ∼15 MΩ-cm. High-resistance water increases signal-to-noise ratio and was shown to have no adverse effects on the health of the fish (Issa et al., 2011). A pair of conductive electrodes placed on either side of the chamber recorded the electric field potentials. Bare electrodes were 1 mm in thickness with 3- to 5-mm metal exposure. Electrodes were connected to an AC differential amplifier (A-M Systems model 1700), allowing the amplification of signals 1000-fold. Electrical signals were low-pass filtered at 300 Hz and high-pass filtered at 1 kHz. Signals were digitized using a Digidata-1322A digitizer and stored using Axoscope software (Molecular Devices). The animals were acclimatized for 30 min before behavioral testing was initiated. Swimming behavior was recorded following acclimation. Immediately after, startle responses were taken and recorded.

Determination of startle sensitivity

Auditory pulses consisting of phasic 1-ms square waves were generated using Audacity open-source audio editor and recorder software (https://www.audacityteam.org/). Sound pulses were amplified with (Realistic, SA-10 solid state stereo amplifier) and connected to a sound speaker (Dell, HK395). Sound intensity was calibrated using a decibel meter (Sinometer, MS6700). Activation of the Mauthner-mediated escape is an all or nothing response with a short latency from stimulus (5–15 ms), and startles were only included if they fell within this range. Non-Mauthner mediated responses >15 ms were excluded (Eaton et al., 2001). Pulse intensity ranged from 70–100 dB re 20 μPa with 5-dB increments and delivery was randomized with a minimum of 2-min intervals between stimuli. Each pulse intensity was delivered three to five times, and response probability for each intensity was tabulated then averaged across animals.

Measurement of swimming activity

Following the 30-min acclimation period and before determination of startle sensitivity, the animal’s swimming behavior was recorded for 1 min. The same methods of data acquisition, amplification, digitization, and storage were used as previously stated. Swimming activity was measured by counting swim bursts with Clampfit software. The “Threshold” function was used for this purpose. A potential was marked as a swim burst if it was at least 8 mV in total amplitude and 30–200 ms in duration. The timing of each swim burst was saved into a Microsoft Excel spreadsheet in reference to the recording start time.

Data analysis and statistics

Startle and swimming data were analyzed using Prism software (GraphPad Inc.) and IBM-SPSS (RRID:SCR_002865). Unless specified otherwise, all comparisons were first subjected to one-way ANOVA (between-subjects factor as group) or mixed design (a mixture of between-subjects factor and repeated measures variables) ANOVA (between-subjects factor as group; within-subjects factor as decibel) or repeated measures ANOVA (within-subjects factors as treatment and decibel) followed by Tukey’s HSD post hoc test or two-sided t test post hoc test for all multiple comparisons. Before using the mixed design or repeated measures ANOVA, sphericity was tested by using Mauchly’s test. When the assumption of sphericity was violated, the degree of freedom in Greenhouse–Geisser correction was used. The one-way ANOVA is considered a relatively robust test against the normality assumption. Startle data were initially curve fitted to multiple equations to determine best fit and highest R² value. The Boltzmann equation produced the best fit with the highest R² values. Thus, nonlinear regressions were performed and fitted using the Boltzmann sigmoidal equation: Y=Bottom+(Top−Bottom)/(1+exp((V50−X)/Slope)) .

Pharmacology

Control experiments were performed on day 14 postpairing during which startle and swimming behaviors were tested and measured. Next day (day 15) fish were injected intramuscularly with the appropriate drug and their startle and swimming behaviors measured again. Immediately after injections the pairs were returned to their pairing tanks but separated with a divider for a period of 1.5 h post-injection of recovery. After the recovery period, the animals were transferred to the behavioral testing chamber and given 30 min of acclimation period.

For injection, glass capillaries were used, having the dimensions 1.0 mm outer diameter (OD) × 0.5 mm inner diameter (ID) × 100 mm in length. These were pulled using Flaming/Brown Micropipette Puller (Model P-87, Sutter Instrument Co). The tip of the micropipette was broken off with a razor blade, before loading with the drug solution. Loaded micropipettes were placed in Pneumatic PicoPump PV 820 for drug administration. A 0.02% tricaine solution was used to anaesthetize the animal before injection. A stock solution of 20 mm each drug was prepared at a concentration of 5 mg/ml. Based on the drugs’ molecular weights, the final concentrations injected were: 25 mm for L-DOPA, 15 mm for SCH23390 (D1 antagonist), 10 mm for bicuculline methobromide (GABA antagonist), and 15 mm for strychnine [glycine receptor (GlyR) antagonist]. These concentrations fall within range of prior reports and our preliminary testing whereby administration of higher concentrations (10 mg/ml) led to increased behavioral deficits (D1 antagonist, bicuculine) and paralysis (strychnine), and prior reports showed that concentration of ≤1.5 mg/ml the drugs had no detrimental effects on motor activity (Roy et al., 2012). All drugs were purchased from Tocris Bioscience except Strychnine, which was purchased from VWR.

Dopamine signaling molecules RNA extraction

RNA extractions were performed using a modified phenol-chloroform protocol using Trizol reagent. Zebrafish were anesthetized with 0.02% MS-222 then euthanized in iced water (10 min). This was followed by dissection and extraction of the brains which were placed in 2-ml microtubes and homogenized by gentle sonication on ice, 400 μl of additional TRIzol was added. The homogenized whole-brain samples were incubated at room temperature for 5 min. A total of 160 μl of chloroform was added, and the samples were mixed. Samples were incubated at 4°C for 20 min with intermittent mixing, followed by centrifugation. The aqueous phase was transferred to a new, sterile tube, and 500-μl ice-cold ethanol was added. The sample was loaded onto a RNeasy spin column (QIAGEN) and spun. Samples were washed with buffer RW1, followed by two washes of Buffer RPE. RNA was eluted with two, 20-μl volumes of RNase-free water into an RNase-free Microfuge tube. RNA extracts were quantified by Nanodrop 2000 (Thermo Scientific), and then stored at −80°C until ready to use.

RT synthesis of cDNA stocks

QIAGEN’s QuantiText Reverse Transcription kit was used for cDNA synthesis of RNA extracts. RNA was brought to a volume of 12 μl with RNase-free water, and 2 μl QIAGEN gDNA wipeout buffer was added. Samples were incubated at 42°C for 2 min to eliminate any DNA contamination. Samples were incubated at 42°C for 30 min for cDNA synthesis, then 3 min at 95°C for enzyme deactivation. Samples were brought up to a final RNA concentration of 10 ng/μl with RNase-free water and stored at −20°C until ready for qPCR.

qPCR primer design and validation

Primers were designed using NCBI Primer-BLAST software and are summarized in Table 1. All primers were designed to span an exon-exon junction, minimize secondary products, produce an amplicon of 75–150 bp, and have an annealing temperature of 60–67°C. Primers were validated to confirm amplicon size, and to determine ideal primer concentration for qPCR reactions. 20 μl PCR reactions were run using GoTaq Flexi DNA Polymerase Reagents from Promega. Primer pairs were assessed at concentrations of 500, 750, and 900 nm final concentrations for amplicon amounts, and the presence of primer-dimer secondary product. Samples were run on a 2% agarose gels containing 0.5× GelRed for visualization under UV light, and pictures were taken with Alpha Innotech’s FluorChem 8900 gel imaging system. Working primer concentrations were determined to be: 500 nm for all primer pairs.

qPCR primer efficiency determination

Standard curves for primer efficiencies were prepared using an initial cDNA amount of 30-ng total cDNA, then 1:10 dilutions were prepared. 20 μl reactions of each cDNA amount and each primer pair were run in triplicate on a 96-well PCR plate from Applied Biosystems. Primer pairs were used in predetermined optimal concentrations, and 10 μl 2× Applied Biosystems Fast SYBR green reaction mix was added for a total volume of 20 μl. Real-time PCR reactions were run on an Applied Biosystems Quantistudio 12k Flex Real-Time PCR System for Ct determination. Primer efficiency standard curves were prepared using the Quantistudio 12k Flex Real-Time PCR System software and are reported in Table 1.

qPCR reaction protocol

qPCR reactions were run using Applied Biosystems Fast SYBR Green Master mix, using Applied Biosystems Quantiostudio 12k Flex Real-Time PCR System. 96 well plates were used with duplicate technical replicates, and 20 μl reactions were prepared. SYBR green detection occurred during each cycle’s amplification stage. Passive reference dye ROX was used for volume standardization. All samples were compared with internal housekeeping gene actb2 for normalization, and isolate fish served as a biological control. Three biological replicates of each social group (dominant, subordinate, and isolate) were analyzed on each 96-well plate. Multiple 96-well plates were run so that each gene had n = 10.

qPCR data analysis

qPCR data were analyzed using the comparative ΔΔCt method. Ct values were determined for each gene of each biological replicate, and experimental genes were compared with internal housekeeping gene actb2 for determination of ΔCt values. ΔΔCt values were determined by normalization of dominant and subordinate social groups to social isolate controls, and percent change was determined from ΔΔCt values using the appropriate primer efficiencies. Percent-change values were transformed to a log base 2 scale for graphing and statistical analysis. Genes percent changes were analyzed for significance by two-sided one sample t test from a theoretical mean of 0, while percent changes between social groups (dominant and subordinate) were analyzed using a two-sided Wilcoxon signed rank test, α = 0.05.

Western blotting

For Western blot analysis four separate Western blot trials were conducted with 10 brains per trial for each social phenotype. Zebrafish were anesthetized with 0.02% MS-222 (1 min) then placed in iced water (10 min). Brains were dissected out, placed in a 1.5-ml microcentrifuge tube, and stored at −20°C until use. Then they were prepared using the total membrane isolation protocol. Brains were homogenized in 1 ml of resuspension buffer: 2.5 ml 2 m sucrose, 2 ml 10 × 10 mm Tris-HCl, 400 μl 0.25 m EDTA, 40 μl 20× protease cocktail inhibitor. The homogenized brains were centrifuged (2000 rpm, 4°C, 10 min) followed by ultracentrifugation (37,000 rpm, 4°C, 60 min). Protein sample concentrations were determined with a Lowery protein assay.

For Western blottings, 10 μg of each protein sample was denatured using 4× buffer containing 10× reducing reagent at 70°C for 10 min and loaded onto a Mini-PROTEAN TGX Precast Protein Gels (Bio-Rad) and run ∼120min at 60V. The proteins were transferred to a 0.2 μm Nitrocellulose membrane using Trans-Blot Turbo RTA Mni Transfer kit (Bio-Rad). The membrane was blocked with PBS containing 5% (w/v) nonfat dry milk and 0.1% (v/v) Tween-20 for 1 h at room temperature and then incubated with primary antibodies overnight at 4°C. After three washes with PBS containing 0.1% (v/v) Tween-20, the membrane was incubated with horseradish peroxidase (HRP)-conjugated goat anti-mouse or goat anti-rabbit IgG (Cell Signaling) secondary antibody for 1 h. After washing, protein detection was performed using HRP detection kit (SuperSignal West Pico, Thermo Fisher Scientific) and visualized using the ChemiDoc Imaging System (Bio-Rad). Bands were quantified with ImageLab software and normalized to β-actin (Bio-Rad). Primary antibodies used were as follows: mouse anti-Drd1b (1:500; Novus Biological), Rabbit anti-DAT (1:500; ProteinTech Group Inc.), β-actin (1:1000; Cell Signaling).

Histology, imaging, and analysis

To determine differences in Drd1b expression in hindbrain glycinergic neurons between social phenotypes, we used the [Tg(GlyT2a:GFP)] transgenic line kindly provided by the Fetcho lab (McLean et al., 2007). Animals’ sex, size, social isolation, pairing, and behavioral analysis were conducted as described above. Animals were then euthanized in a heavy dose of MS-222 (0.25% by volume). Animals were pinned on Sylgard lined Petri dish (Dow Corning) and brains extracted. Brains were fixed overnight in 4% PFA at 4°C, washed 3 × 5 min in 0.1% PBS-T then placed in 30% sucrose for an additional night at 4°C. Brains were then placed in plastic molds containing O.C.T. tissue compound (Tissue-Tek) and frozen in liquid nitrogen and stored at −20°C. Frozen brains were sliced into 60-μm slices using a cryostat (Leica Biosystems) and mounted on Superfrost plus slides (Fisher Brand). Brain slices were incubated in a blocking buffer for 2 h, followed by an incubation in a Rabbit polyclonal Drd1b primary antibody (ThermoFisher, catalog #PA5-33477) 1:500 overnight at 4°C. Primary antibody was washed in PBS 3 × 5 min, then incubated in an Alexa-555 goat anti-rabbit secondary antibody (Invitrogen, catalog #A21428) 1:1000, finally washed in PBS 3 × 5 min. Slices then were mounted with Antifade reagent and cover slip, then sealed 24 h later. All images were acquired with a 40× oil immersion using Carl Zeiss LSM 800 upright confocal microscope whose acquisition parameters were kept similar across experimental conditions.

Confocal image analysis

All acquired images from experimental and control groups were coded and analyzed blindly. Hindbrain glycinergic neurons were identified visually based on their morphologic location in relation to spinal midline and M-cell location. All glycinergic and drd1b expressing neurons were quantified using Imaris software (version 9.3; Oxford Instruments plc) surface and spots functions. This was followed by identification and quantification of all glycinergic neurons that specifically express the Drd1b using the Imaris co-localization function. Ratios of Drd1b expressing cells to total number of glycinergic neurons was tallied and data plotted using GraphPad Prism software. Representative confocal images were imported into Adobe Photoshop, stack projections (60-μm slices) were generated, and dust and scratches noise filter was applied to red channel of all images to eliminate background noise.

drd1b^(−/−) knock-out (KO) line was kindly provided by Teresa Nicolson (Stanford University) and originally created at the Wellcome Sanger Institute (Cambridgeshire, United Kingdom). The mutant line (allele: sa16476) contains a nonsense point mutation (G > A) in the second exon on chromosome 9 that affects amino acid number 314 leading to a truncated nonfunctional protein (Ensembl ID: ENSDARG00000038918); ZFIN ID: ZDB-GENE-070524-2.

Neurocomputational model

In the present study, we consider the behavior of M-cell that receives excitatory and inhibitory inputs, which is modulated by dopamine. The model network is composed of one excitatory neuron, two inhibitory neurons (GABAergic and glycinergic), and one M-cell. The excitatory neuron receives sensory inputs and excites all the other neurons. GABAergic neuron inhibits glycinergic neuron and M-cell. Glycinergic neuron inhibits M-cell only.

All neurons were modeled as a conductance-based modified Morris–Lecar neuronal model (Morris and Lecar, 1981; Izhikevich, 2007; G.B. Ermentrout and Terman, 2010) with additional calcium-dependent potassium current. The membrane potential of each cell obeys the following current balance equation: Cdvdt=−ICa−IK−IL−IKCa−Isyn + Iapp(t), where IK=gKn(v−vK), ICa=gCam∞(v)(v−vCa), IKCa=gKCa{[Ca][Ca] + k1}(v−vK),IL=gL(v−vL) represent the potassium, calcium, calcium-dependent potassium, and leak currents, respectively. Isyn represents the synaptic currents and Iapp(t) is the applied current (see below). [Ca] represents intracellular calcium concentration. For all neurons, gK=8 , vK=−84 , gL=2 , vL=−60 , gCa=4 , vCa=120 , and k1=10 . In M-cell, gKCa=0.3 and C=1 . In other neurons, gKCa=0.25 and C=20 .

m∞ is an instantaneous voltage-dependent gating variable for the calcium current where m∞(v)=0.5(1 + tanh(v−v1v2)), with v1=−1.2 and v2=18.

The concentration of intracellular Ca2+ is governed by the calcium balance equation d[Ca]dt=ε(−μ ICa−kCa[Ca]), where ε=0.005 , μ=0.19 . In M-cell, kCa=0.9 while kCa=1 in all other neurons.

n is a gating variable for the potassium current obeying dndt=ϕ(n∞(v)−n)τn(v), n∞(v)=0.5(1 + tanh(v−v3v4)), τn(v)=1/cosh(v−v32v4), where ϕ=0.23 , v3=12 , and v4=17 .

In an excitatory cell and inhibitory cells, synaptic variable, s, is modeled by an equation for the fraction of activated channel dsdt=αs∞(v)(1−s)−βs, where s∞(v)=1/(1 + exp(−v + θsσs)) with θs=0 and σs=4 . Here, α=15 and β=0.1 in an excitatory cell, α=4 and β=0.08 in a GABAergic cell, and α=8 and β=0.08 in a glycinergic cell.

The synaptic current (Isyn ), the sum of synaptic inputs from other cells, is given by Isyn=gsyn(v−vsyn)∑jsj, where the summation is over s variables from all neurons projecting to a given neuron.

In the current neuronal network, an excitatory cell does not receive any synaptic input but receives an external stimulus to simulate the effect from the hair cell or VIIIth nerve fiber. Thus, in an excitatory cell, Isyn=0 and Iapp(t)=IE0+WE*I(τ) , where IE0 is a baseline external input, WE is the stimulus strength and I(τ) is the stimulus at time τ. IE0=43.9 and WE∈{10,15,20,...,55,60} in the model. Here, I(τ) resembles the unit square pulse with height 1.

GABAergic cell does not receive any external stimulus but receives a synaptic input from the excitatory cell. Thus, in a GABAergic cell, Isyn=gE→GA(v−vE→GA)sE , where sE is the synaptic variable from the excitatory cell and Iapp(t)=IGA0=36 , a baseline external input to the GABAergic cell. We use vE→GA=40andgE→GA=0.3 in the model.

Glycinergic cell receives synaptic inputs from the excitatory cell and the GABAergic cell. Thus, in glycinergic cell, Isyn=gE→GL(v−vE→GL)sE+gGA→GL(v−vGA→GL)sGA , where sE is a synaptic variable from the excitatory cell and sGA is a synaptic variable from GABAergic cell. Iapp(t)=IGL0=36 , a baseline external input to the glycinergic cell. To reflect dopaminergic modulation over excitatory connection in the model, the excitatory maximal synaptic conductance (gE→GL ) was modified as explained below.

In M-cell, which receives synaptic inputs from all the other cells in the network, Isyn=gE→M(v−vE→M)sE + (gGA→MsGA + gGL→MsGL)(v−vG→M) + gM→M(v−vM→M)sM, where sM is the synaptic variable from another M-cell, which is assumed to be a small constant (sM=0.029) in the current study. sGL is a synaptic variable from the glycinergic cell. Iapp(t)=IM0=19.5 , a baseline external input to the M-cell. Now, calcium is known to modulate inhibitory presynaptic neurotransmitter release via retrograde signaling (Diana and Bregestovski, 2005). In the model M-cell, we assumed that intracellular calcium level reciprocally modulates the presynaptic inputs to the M-cell and Isyn is updated as follows: Isyn=gE→MgI(v−vE→M)sE + gI(gGA→MsGA + gGL→MsGL)(v−vG→M) + gM→MgI(v−vM→M)sM, where gI obeys the following equation dgIdt=gImax[Ca] + k2−gIρ, where gImax is the maximal gI value, ρ is the time constant of gI , and [Ca] is the intracellular calcium concentration of the M-cell. Here, parameter values are gE→M=0.24,vE→ M=40, vG→M=−50,gM→M=0.1,vM→M=−50,gImax=15,k2=10,ρ=10000,gGA→M=0.4,gGL→M=0.2 . Now we implemented the dopamine modulation of synaptic inputs in the glycinergic cell, and the M-cell. In the current network, we assume that D₁ enhances the input that a cell receives.

In the glycinergic cell, the excitatory input is enhanced by D₁, which is reflected in Isyn as follows: Isyn=gE→GL(1 + D1GL)(v−vE→GL)sE + gGA→GL(v−vGA→GL)sGA, where gE→GL=0.3,vE→GL=40 , and vGA→GL=−50. We also used D1GL=0.65 and gGA→GL=0.2 for dominant-like model, and D1GL=0.25 and gGA→GL=0.4 for subordinate-like model.

Lastly, in the M-cell, D₁ enhances the amplitudes of excitatory input, which is reflected in Isyn as follows: Isyn=gE→MgI(1 + D1M)(v−vE→M)sE + gI(gGA→MsGA + gGL→MsGL)(v−vG→M) + gM→MgI(v−vM→M)sM where D1M=0.015 .

In summary, we varied three parameters (D1GL , gGA→GL , β) between dominant-like model and subordinate-like model: D1GL=0.65 , gGA→GL=0.2 for dominant-like model and D1GL=0.25, gGA→GL=0.4 for subordinate-like model. We also varied the decay time constant β in the synaptic variable for the glycinergic input to simulate the blockage of the GABAergic receptor: β=0.024 for dominant-like model and β=0.0072 for subordinate-like model. These parameters were chosen based on the model hypothesis that dominant-like model has stronger dopaminergic→glycinergic→M, while subordinate-like model has stronger GABAergic→glycinergic pathways.

Simulations were performed on a personal computer using the software XPP (B. Ermentrout, 2002). The numerical method used was an adaptive-step fourth order Runge–Kutta method with a step size 0.01 ms.

Code accessibility

The code described is deposited and freely available online at https://modeldb.science/2015414 and the code is available as Extended Data 1.