Data analysis and statistical tests

Fiber-photometry recording data were exported as MATLAB. Mat files from Spike2 software for further analysis. All the raw data were smoothed with a moving average filter (20 ms span) and then segmented and aligned according to the onset of behavioral events within individual trials or bouts. The fluorescence change (ΔF/F) values were calculated as (F−F₀)/F₀, where F₀ is the baseline fluorescence signals averaged over a 1.5 s-long control time window (typically set 0.5 s) prior to a trigger event. To analyze the responses during social interaction, the control time window was set 3.5 s before interaction onset to minimize potential chasing-induced effects. ΔF/F values are presented as heatmaps or as average plots with a shaded area indicating the SEM. Mouse locomotor activity was analyzed using a custom video tracking software developed in house using MATLAB. We plotted the locomotor speed during the 4 s time window following cue onset (0 – 4 s) across individual trials in a behavior session. We then averaged the value per five trials for aversive conditioning and per ten trials for appetitive conditioning.

For the in vivo electrophysiological data, the spikes were sorted off-line with the Spike2 program. Single units were isolated using principal component analysis (PCA) of the spike waveforms that had signal-to-noise ratios of at least 2:1. PETH of spike firing rates (bin width 50 ms) were smoothed with a Gaussian kernel (σ = 50 ms) and then presented either as heatmaps or as average plots. To calculate standard scores, we used the mean firing frequency of a control period. Hierarchical clustering was carried out by reducing the dimensionality of standardized firing activity via principal component analysis (PCA). The first three major principle components (PCs) were then used to calculate a Euclidean distance metric. The complete agglomeration method was applied to build the hierarchy of clusters. Minor adjustments were made by sorting the clusters in the descending order based on the total Z-score values between 5 and 10 s from cue onset.

We applied multivariate permutation tests to analyze the statistical significance of the event-related fluorescence (ERF) change or peri-event time histograms (PETH) of spike firing rates (1000 permutations, α level of 0.05). The null distribution was retrieved from the maximum absolute T-score of all permutations to correct multiple comparisons in two-tailed tests. A series of inferential p values at each time point were generated and the results were superimposed on the average ERF or PETH curves with red and blue lines indicating statistically significant (p<0.05) increases or decreases, respectively.

We plotted receiver operating characteristic (ROC) curves and calculated the area under the curve (AUC) for ERF throughout each trial by comparing the ERF of a 200 ms test window (50 ms advance step) to those in a control time window (200 ms) that occurred 1.8 s preceding the trial onset (−2 to −1.8 s). ROC values >0.5 indicate activation, and values <0.5 indicate inhibition. Differences in the ROC values between the Cue1 and Cue2 in Figure 6E–H were calculated by comparing ERF numbers during the same time windows (200 ms width, 50 ms advance step) throughout the entire trial. ROC values of 1 indicate complete selectivity for the Cue1 stimulus, and ROC values of 0 indicate complete selectivity for the Cue2 stimulus. Permutation tests with 1000 permutations were used to determine the statistical significance of the response strength and selectivity of aversive stimuli or rewards. We performed Kolmogorov–Smirnov tests for the statistical significance of the differences between the cumulative probability distributions. Similarly, we calculated the AUC for event-related fluorescent changes (ΔF/F) as the sum of Ca²⁺ transients.

We performed hierarchical clustering of the reward-related responses in three steps. We first applied principle component analysis (PCA) to reduce the dimensionality of standardized GCaMP signals and firing activity. We then used the first three major principle components (PCs) to define a Euclidean distance metric. Finally, we applied the complete agglomeration method to construct the hierarchy of clusters and plot dendrograms in MATLAB.